diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index b95476fae3..1f768b4e88 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -8,6 +8,9 @@ * Author: Chao Peng (ANL), 09/27/2020 */ +#include +#include +#include #include #include #include @@ -27,10 +30,11 @@ #include #include #include +#include #include #include -#include "CalorimeterClusterRecoCoG.h" +#include "algorithms/calorimetry/CalorimeterClusterRecoCoG.h" #include "algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h" namespace eicrecon { @@ -38,6 +42,19 @@ namespace eicrecon { using namespace dd4hep; void CalorimeterClusterRecoCoG::init() { + + // (Re)acquire detector and, if configured, IDDescriptor for system-aware behavior + m_detector = algorithms::GeoSvc::instance().detector(); + if (!m_detector) { + error("Failed to get detector from GeoSvc"); + throw std::runtime_error("Detector not available"); + } + + if (!m_cfg.readout.empty()) { + m_idSpec = m_detector->readout(m_cfg.readout).idSpec(); + } else { + warning("No readout configured; system-aware ScFi/Imaging weighting is disabled."); + } // select weighting method std::string ew = m_cfg.energyWeight; // make it case-insensitive @@ -92,6 +109,24 @@ void CalorimeterClusterRecoCoG::process(const CalorimeterClusterRecoCoG::Input& } } +// System ID from cell ID + +static const int SYSID_SCIFI = 105; +static const int SYSID_IMAGING = 101; + +inline int getSystemID(const edm4eic::CalorimeterHit& hit, const dd4hep::IDDescriptor& m_idSpec) { + static thread_local const auto* sys_field = m_idSpec.field("system"); + return sys_field->value(hit.getCellID()); +} + +inline bool isSciFiHit(const edm4eic::CalorimeterHit& hit, const dd4hep::IDDescriptor& m_idSpec) { + return getSystemID(hit, m_idSpec) == SYSID_SCIFI; +} + +inline bool isImagingHit(const edm4eic::CalorimeterHit& hit, const dd4hep::IDDescriptor& m_idSpec) { + return getSystemID(hit, m_idSpec) == SYSID_IMAGING; +} + std::optional CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoCluster& pcl) const { edm4eic::MutableCluster cl; @@ -111,19 +146,54 @@ CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoCluster& pcl) const { float maxHitEta = std::numeric_limits::min(); auto time = 0; auto timeError = 0; + + bool hasSciFi = false; + bool hasImaging = false; + + for (const auto& hit : pcl.getHits()) { + if (isSciFiHit(hit, m_idSpec)) { + hasSciFi = true; + } + if (isImagingHit(hit, m_idSpec)) { + hasImaging = true; + } + } + + bool specialMode = hasSciFi && hasImaging; + + // ----------------------------------------------------------------------------------- + // ScFi hit alone contributes to energy weight & for Img hit energy weight = 0 + // ----------------------------------------------------------------------------------- + for (unsigned i = 0; i < pcl.getHits().size(); ++i) { + const auto& hit = pcl.getHits()[i]; const auto weight = pcl.getWeights()[i]; - debug("hit energy = {} hit weight: {}", hit.getEnergy(), weight); - auto energy = hit.getEnergy() * weight; + + float energy = 0.0f; + + if (specialMode) { + + if (isSciFiHit(hit, m_idSpec)) { + energy = hit.getEnergy() * weight; + } else if (isImagingHit(hit, m_idSpec)) { + energy = 0.0f; // Imaging has no energy weight contribution + } + } else { + energy = hit.getEnergy() * weight; + } + totalE += energy; time += (hit.getTime() - time) * energy / totalE; + cl.addToHits(hit); cl.addToHitContributions(energy); - const float eta = edm4hep::utils::eta(hit.getPosition()); - minHitEta = std::min(eta, minHitEta); - maxHitEta = std::max(eta, maxHitEta); + + float eta = edm4hep::utils::eta(hit.getPosition()); + minHitEta = std::min(minHitEta, eta); + maxHitEta = std::max(maxHitEta, eta); } + cl.setEnergy(totalE / m_cfg.sampFrac); cl.setEnergyError(0.); cl.setTime(time); @@ -142,13 +212,30 @@ CalorimeterClusterRecoCoG::reconstruct(const edm4eic::ProtoCluster& pcl) const { } } + // -------------------------------------------------------------------------------------- + // Imaging hit alone contributes to Position weight & for ScFi hit position weight = 0 + // -------------------------------------------------------------------------------------- for (unsigned i = 0; i < pcl.getHits().size(); ++i) { - const auto& hit = pcl.getHits()[i]; - const auto weight = pcl.getWeights()[i]; - // _DBG_<<" -- weight = " << weight << " E=" << hit.getEnergy() << " totalE=" < +#include #include +#include #include #include #include @@ -23,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -81,6 +85,10 @@ class CalorimeterClusterRecoCoG : public CalorimeterClusterRecoCoGAlgorithm, private: std::function weightFunc; + // Pointer to the geometry service + dd4hep::IDDescriptor m_idSpec; + const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()}; + private: std::optional reconstruct(const edm4eic::ProtoCluster& pcl) const; void associate(const edm4eic::Cluster& cl, diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h index df38aca48c..a18f1d834d 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h @@ -10,7 +10,7 @@ namespace eicrecon { struct CalorimeterClusterRecoCoGConfig { - + std::string readout = ""; std::string energyWeight; double sampFrac = 1.; diff --git a/src/algorithms/calorimetry/ImagingTopoCluster.cc b/src/algorithms/calorimetry/ImagingTopoCluster.cc index 4022ed2523..6988e92c0a 100644 --- a/src/algorithms/calorimetry/ImagingTopoCluster.cc +++ b/src/algorithms/calorimetry/ImagingTopoCluster.cc @@ -22,11 +22,16 @@ #include "algorithms/calorimetry/ImagingTopoCluster.h" #include +#include +#include #include #include #include +#include +#include #include #include +#include #include #include #include @@ -43,6 +48,14 @@ template struct multilambda : L... { void ImagingTopoCluster::init() { + m_detector = algorithms::GeoSvc::instance().detector(); + if (!m_detector) { + error("Failed to get detector from GeometrySvc for readout '{}'", m_cfg.readout); + throw std::runtime_error("Detector not available"); + } + + m_idSpec = m_detector->readout(m_cfg.readout).idSpec(); + multilambda _toDouble = { [](const std::string& v) { return dd4hep::_toDouble(v); }, [](const double& v) { return v; }, @@ -62,16 +75,47 @@ void ImagingTopoCluster::init() { sameLayerDistXY[1] = std::visit(_toDouble, m_cfg.sameLayerDistXY[1]) / dd4hep::mm; diffLayerDistXY[0] = std::visit(_toDouble, m_cfg.diffLayerDistXY[0]) / dd4hep::mm; diffLayerDistXY[1] = std::visit(_toDouble, m_cfg.diffLayerDistXY[1]) / dd4hep::mm; + + ScFi_sameLayerDistXY[0] = std::visit(_toDouble, m_cfg.ScFi_sameLayerDistXY[0]) / dd4hep::mm; + ScFi_sameLayerDistXY[1] = std::visit(_toDouble, m_cfg.ScFi_sameLayerDistXY[1]) / dd4hep::mm; + ScFi_diffLayerDistXY[0] = std::visit(_toDouble, m_cfg.ScFi_diffLayerDistXY[0]) / dd4hep::mm; + ScFi_diffLayerDistXY[1] = std::visit(_toDouble, m_cfg.ScFi_diffLayerDistXY[1]) / dd4hep::mm; + + ScFi_sameLayerDistXZ[0] = m_cfg.ScFi_sameLayerDistXZ[0] / dd4hep::mm; + ScFi_sameLayerDistXZ[1] = m_cfg.ScFi_sameLayerDistXZ[1] / dd4hep::mm; + ScFi_diffLayerDistXZ[0] = m_cfg.ScFi_diffLayerDistXZ[0] / dd4hep::mm; + ScFi_diffLayerDistXZ[1] = m_cfg.ScFi_diffLayerDistXZ[1] / dd4hep::mm; + + Img_sameLayerDistXY[0] = std::visit(_toDouble, m_cfg.Img_sameLayerDistXY[0]) / dd4hep::mm; + Img_sameLayerDistXY[1] = std::visit(_toDouble, m_cfg.Img_sameLayerDistXY[1]) / dd4hep::mm; + Img_diffLayerDistXY[0] = std::visit(_toDouble, m_cfg.Img_diffLayerDistXY[0]) / dd4hep::mm; + Img_diffLayerDistXY[1] = std::visit(_toDouble, m_cfg.Img_diffLayerDistXY[1]) / dd4hep::mm; + if (m_cfg.sameLayerDistXYZ.size() >= 3) { - sameLayerDistXYZ[0] = std::visit(_toDouble, m_cfg.sameLayerDistXYZ[0]) / dd4hep::mm; - sameLayerDistXYZ[1] = std::visit(_toDouble, m_cfg.sameLayerDistXYZ[1]) / dd4hep::mm; - sameLayerDistXYZ[2] = std::visit(_toDouble, m_cfg.sameLayerDistXYZ[2]) / dd4hep::mm; + sameLayerDistXYZ[0] = m_cfg.sameLayerDistXYZ[0] / dd4hep::mm; + sameLayerDistXYZ[1] = m_cfg.sameLayerDistXYZ[1] / dd4hep::mm; + sameLayerDistXYZ[2] = m_cfg.sameLayerDistXYZ[2] / dd4hep::mm; } if (m_cfg.diffLayerDistXYZ.size() >= 3) { - diffLayerDistXYZ[0] = std::visit(_toDouble, m_cfg.diffLayerDistXYZ[0]) / dd4hep::mm; - diffLayerDistXYZ[1] = std::visit(_toDouble, m_cfg.diffLayerDistXYZ[1]) / dd4hep::mm; - diffLayerDistXYZ[2] = std::visit(_toDouble, m_cfg.diffLayerDistXYZ[2]) / dd4hep::mm; + diffLayerDistXYZ[0] = m_cfg.diffLayerDistXYZ[0] / dd4hep::mm; + diffLayerDistXYZ[1] = m_cfg.diffLayerDistXYZ[1] / dd4hep::mm; + diffLayerDistXYZ[2] = m_cfg.diffLayerDistXYZ[2] / dd4hep::mm; } + + ScFi_sameLayerDistXYZ[0] = m_cfg.ScFi_sameLayerDistXYZ[0] / dd4hep::mm; + ScFi_sameLayerDistXYZ[1] = m_cfg.ScFi_sameLayerDistXYZ[1] / dd4hep::mm; + ScFi_sameLayerDistXYZ[2] = m_cfg.ScFi_sameLayerDistXYZ[2] / dd4hep::mm; + ScFi_diffLayerDistXYZ[0] = m_cfg.ScFi_diffLayerDistXYZ[0] / dd4hep::mm; + ScFi_diffLayerDistXYZ[1] = m_cfg.ScFi_diffLayerDistXYZ[1] / dd4hep::mm; + ScFi_diffLayerDistXYZ[2] = m_cfg.ScFi_diffLayerDistXYZ[2] / dd4hep::mm; + + Img_sameLayerDistXYZ[0] = m_cfg.Img_sameLayerDistXYZ[0] / dd4hep::mm; + Img_sameLayerDistXYZ[1] = m_cfg.Img_sameLayerDistXYZ[1] / dd4hep::mm; + Img_sameLayerDistXYZ[2] = m_cfg.Img_sameLayerDistXYZ[2] / dd4hep::mm; + Img_diffLayerDistXYZ[0] = m_cfg.Img_diffLayerDistXYZ[0] / dd4hep::mm; + Img_diffLayerDistXYZ[1] = m_cfg.Img_diffLayerDistXYZ[1] / dd4hep::mm; + Img_diffLayerDistXYZ[2] = m_cfg.Img_diffLayerDistXYZ[2] / dd4hep::mm; + sameLayerDistEtaPhi[0] = m_cfg.sameLayerDistEtaPhi[0]; sameLayerDistEtaPhi[1] = m_cfg.sameLayerDistEtaPhi[1] / dd4hep::rad; diffLayerDistEtaPhi[0] = m_cfg.diffLayerDistEtaPhi[0]; @@ -81,10 +125,26 @@ void ImagingTopoCluster::init() { diffLayerDistTZ[0] = m_cfg.diffLayerDistTZ[0] / dd4hep::mm; diffLayerDistTZ[1] = m_cfg.diffLayerDistTZ[1] / dd4hep::mm; - sectorDist = m_cfg.sectorDist / dd4hep::mm; - minClusterHitEdep = m_cfg.minClusterHitEdep / dd4hep::GeV; - minClusterCenterEdep = m_cfg.minClusterCenterEdep / dd4hep::GeV; - minClusterEdep = m_cfg.minClusterEdep / dd4hep::GeV; + Img_sameLayerDistEtaPhi[0] = m_cfg.Img_sameLayerDistEtaPhi[0]; + Img_sameLayerDistEtaPhi[1] = m_cfg.Img_sameLayerDistEtaPhi[1] / dd4hep::rad; + Img_diffLayerDistEtaPhi[0] = m_cfg.Img_diffLayerDistEtaPhi[0]; + Img_diffLayerDistEtaPhi[1] = m_cfg.Img_diffLayerDistEtaPhi[1] / dd4hep::rad; + Img_sameLayerDistTZ[0] = m_cfg.Img_sameLayerDistTZ[0] / dd4hep::mm; + Img_sameLayerDistTZ[1] = m_cfg.Img_sameLayerDistTZ[1] / dd4hep::mm; + Img_diffLayerDistTZ[0] = m_cfg.Img_diffLayerDistTZ[0] / dd4hep::mm; + Img_diffLayerDistTZ[1] = m_cfg.Img_diffLayerDistTZ[1] / dd4hep::mm; + + cross_system_DistXYZ[0] = m_cfg.cross_system_DistXYZ[0] / dd4hep::mm; + cross_system_DistXYZ[1] = m_cfg.cross_system_DistXYZ[1] / dd4hep::mm; + cross_system_DistXYZ[2] = m_cfg.cross_system_DistXYZ[2] / dd4hep::mm; + + sectorDist = m_cfg.sectorDist / dd4hep::mm; + cross_system_sectorDist = m_cfg.cross_system_sectorDist / dd4hep::mm; + ScFi_sectorDist = m_cfg.ScFi_sectorDist / dd4hep::mm; + Img_sectorDist = m_cfg.Img_sectorDist / dd4hep::mm; + minClusterHitEdep = m_cfg.minClusterHitEdep / dd4hep::GeV; + minClusterCenterEdep = m_cfg.minClusterCenterEdep / dd4hep::GeV; + minClusterEdep = m_cfg.minClusterEdep / dd4hep::GeV; // same layer clustering parameters switch (m_cfg.sameLayerMode) { @@ -154,9 +214,10 @@ void ImagingTopoCluster::init() { "Global [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", m_cfg.neighbourLayersRange, diffLayerDistXY[0], diffLayerDistXY[1]); break; + case ImagingTopoClusterConfig::ELayerMode::xyz: if (m_cfg.diffLayerDistXYZ.size() != 3) { - const std::string msg = "Expected 3 values (x_dist, y_dist, z_dist) for diffLayerDistXYZ"; + const std::string msg = "Expected 3 values (x_dist, y_dist, y_dist) for diffLayerDistXYZ"; error(msg); throw std::runtime_error(msg); } @@ -181,6 +242,221 @@ void ImagingTopoCluster::init() { info("Neighbour sectors clustering (different sector): " "Global distance between hits <= {:.4f} mm.", sectorDist); + + //ScFi_layer mode + switch (m_cfg.ScFi_sameLayerMode) { + + case ImagingTopoClusterConfig::ELayerMode::xy: + if (m_cfg.ScFi_sameLayerDistXY.size() != 2) { + const std::string msg = "Expected 2 values (x_dist, y_dist) for sameLayerDistXY"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Same-layer clustering (same sector and same layer): " + "Local [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", + ScFi_sameLayerDistXY[0], ScFi_sameLayerDistXY[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xz: + if (m_cfg.ScFi_sameLayerDistXZ.size() != 2) { + const std::string msg = "Expected 2 values (x_dist, z_dist) for sameLayerDistXZ"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Same-layer clustering (same sector and same layer): " + "Local [x, z] distance between hits <= [{:.4f} mm, {:.4f} mm].", + ScFi_sameLayerDistXZ[0], ScFi_sameLayerDistXZ[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xyz: + if (m_cfg.ScFi_sameLayerDistXYZ.size() != 3) { + const std::string msg = + "Expected 3 values (x_dist, y_dist, z_dist) for ScFi_sameLayerDistXYZ"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Same-layer clustering (same sector and same layer): " + "Local [x, y, z] distance between hits <= [{:.4f} mm, {:.4f} mm, {:.4f} mm].", + ScFi_sameLayerDistXYZ[0], ScFi_sameLayerDistXYZ[1], ScFi_sameLayerDistXYZ[2]); + break; + case ImagingTopoClusterConfig::ELayerMode::etaphi: + if (m_cfg.ScFi_sameLayerDistEtaPhi.size() != 2) { + const std::string msg = "Expected 2 values (eta_dist, phi_dist) for ScFi_sameLayerDistEtaPhi"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Same-layer clustering (same sector and same layer): " + "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].", + ScFi_sameLayerDistEtaPhi[0], ScFi_sameLayerDistEtaPhi[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::tz: + if (m_cfg.ScFi_sameLayerDistTZ.size() != 2) { + const std::string msg = "Expected 2 values (t_dist, z_dist) for ScFi_sameLayerDistTZ"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Same-layer clustering (same sector and same layer): " + "Global [t, z] distance between hits <= [{:.4f} mm, {:.4f} mm].", + ScFi_sameLayerDistTZ[0], ScFi_sameLayerDistTZ[1]); + break; + default: + throw std::runtime_error("Unknown same-layer mode."); + } + + switch (m_cfg.ScFi_diffLayerMode) { + case ImagingTopoClusterConfig::ELayerMode::etaphi: + if (m_cfg.ScFi_diffLayerDistEtaPhi.size() != 2) { + const std::string msg = "Expected 2 values (eta_dist, phi_dist) for ScFi_diffLayerDistEtaPhi"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].", + m_cfg.ScFi_neighbourLayersRange, ScFi_diffLayerDistEtaPhi[0], ScFi_diffLayerDistEtaPhi[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xy: + if (m_cfg.ScFi_diffLayerDistXY.size() != 2) { + const std::string msg = "Expected 2 values (x_dist, y_dist) for ScFi_diffLayerDistXY"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", + m_cfg.ScFi_neighbourLayersRange, ScFi_diffLayerDistXY[0], ScFi_diffLayerDistXY[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xz: + if (m_cfg.ScFi_diffLayerDistXZ.size() != 2) { + const std::string msg = "Expected 2 values (x_dist, z_dist) for ScFi_diffLayerDistXZ"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [x, z] distance between hits <= [{:.4f} mm, {:.4f} mm].", + m_cfg.ScFi_neighbourLayersRange, ScFi_diffLayerDistXZ[0], ScFi_diffLayerDistXZ[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xyz: + if (m_cfg.ScFi_diffLayerDistXYZ.size() != 3) { + const std::string msg = + "Expected 3 values (x_dist, y_dist, y_dist) for ScFi_diffLayerDistXYZ"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [x, y, z] distance between hits <= [{:.4f} mm, {:.4f} mm, {:.4f} mm].", + m_cfg.ScFi_neighbourLayersRange, ScFi_diffLayerDistXYZ[0], ScFi_diffLayerDistXYZ[1], + ScFi_diffLayerDistXYZ[2]); + break; + case ImagingTopoClusterConfig::ELayerMode::tz: + if (m_cfg.ScFi_diffLayerDistTZ.size() != 2) { + const std::string msg = "Expected 2 values (t_dist, z_dist) for ScFi_diffLayerDistTZ"; + error(msg); + throw std::runtime_error(msg); + } + info("ScFi_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [t, z] distance between hits <= [{:.4f} mm, {:.4f} mm].", + m_cfg.ScFi_neighbourLayersRange, ScFi_diffLayerDistTZ[0], ScFi_diffLayerDistTZ[1]); + break; + default: + error("Unknown different-layer mode."); + throw std::runtime_error("Unknown different-layer mode."); + } + info("ScFi_Neighbour sectors clustering (different sector): " + "Global distance between hits <= {:.4f} mm.", + ScFi_sectorDist); + + //Imaging Layer Mode + switch (m_cfg.Img_sameLayerMode) { + case ImagingTopoClusterConfig::ELayerMode::xy: + if (m_cfg.Img_sameLayerDistXY.size() != 2) { + const std::string msg = "Expected 2 values (x_dist, y_dist) for Img_sameLayerDistXY"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Same-layer clustering (same sector and same layer): " + "Local [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", + Img_sameLayerDistXY[0], Img_sameLayerDistXY[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xyz: + if (m_cfg.Img_sameLayerDistXYZ.size() != 3) { + const std::string msg = "Expected 3 values (x_dist, y_dist, z_dist) for Img_sameLayerDistXYZ"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Same-layer clustering (same sector and same layer): " + "Local [x, y, z] distance between hits <= [{:.4f} mm, {:.4f} mm, {:.4f} mm].", + Img_sameLayerDistXYZ[0], Img_sameLayerDistXYZ[1], Img_sameLayerDistXYZ[2]); + break; + case ImagingTopoClusterConfig::ELayerMode::etaphi: + if (m_cfg.Img_sameLayerDistEtaPhi.size() != 2) { + const std::string msg = "Expected 2 values (eta_dist, phi_dist) for Img_sameLayerDistEtaPhi"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Same-layer clustering (same sector and same layer): " + "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].", + Img_sameLayerDistEtaPhi[0], Img_sameLayerDistEtaPhi[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::tz: + if (m_cfg.Img_sameLayerDistTZ.size() != 2) { + const std::string msg = "Expected 2 values (t_dist, z_dist) for Img_sameLayerDistTZ"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Same-layer clustering (same sector and same layer): " + "Global [t, z] distance between hits <= [{:.4f} mm, {:.4f} mm].", + Img_sameLayerDistTZ[0], Img_sameLayerDistTZ[1]); + break; + default: + throw std::runtime_error("Unknown same-layer mode."); + } + + switch (m_cfg.Img_diffLayerMode) { + case ImagingTopoClusterConfig::ELayerMode::etaphi: + if (m_cfg.Img_diffLayerDistEtaPhi.size() != 2) { + const std::string msg = "Expected 2 values (eta_dist, phi_dist) for Img_diffLayerDistEtaPhi"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].", + m_cfg.neighbourLayersRange, Img_diffLayerDistEtaPhi[0], Img_diffLayerDistEtaPhi[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xy: + if (m_cfg.Img_diffLayerDistXY.size() != 2) { + const std::string msg = "Expected 2 values (x_dist, y_dist) for Img_diffLayerDistXY"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", + m_cfg.neighbourLayersRange, Img_diffLayerDistXY[0], Img_diffLayerDistXY[1]); + break; + case ImagingTopoClusterConfig::ELayerMode::xyz: + if (m_cfg.Img_diffLayerDistXYZ.size() != 3) { + const std::string msg = "Expected 3 values (x_dist, y_dist, y_dist) for Img_diffLayerDistXYZ"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [x, y, z] distance between hits <= [{:.4f} mm, {:.4f} mm, {:.4f} mm].", + m_cfg.neighbourLayersRange, Img_diffLayerDistXYZ[0], Img_diffLayerDistXYZ[1], + Img_diffLayerDistXYZ[2]); + break; + case ImagingTopoClusterConfig::ELayerMode::tz: + if (m_cfg.Img_diffLayerDistTZ.size() != 2) { + const std::string msg = "Expected 2 values (t_dist, z_dist) for Img_diffLayerDistTZ"; + error(msg); + throw std::runtime_error(msg); + } + info("Img_Neighbour layers clustering (same sector and layer id within +- {:d}): " + "Global [t, z] distance between hits <= [{:.4f} mm, {:.4f} mm].", + m_cfg.neighbourLayersRange, Img_diffLayerDistTZ[0], Img_diffLayerDistTZ[1]); + break; + default: + error("Unknown different-layer mode."); + throw std::runtime_error("Unknown different-layer mode."); + } + info("Img_Neighbour sectors clustering (different sector): " + "Global distance between hits <= {:.4f} mm.", + Img_sectorDist); } void ImagingTopoCluster::process(const Input& input, const Output& output) const { @@ -188,172 +464,546 @@ void ImagingTopoCluster::process(const Input& input, const Output& output) const const auto [hits] = input; auto [proto] = output; - // Sort hit indices (podio collections do not support std::sort) - auto compare = [&hits](const auto& a, const auto& b) { - // if !(a < b) and !(b < a), then a and b are equivalent - // and only one of them will be allowed in a set - if ((*hits)[a].getLayer() == (*hits)[b].getLayer()) { - return (*hits)[a].getObjectID().index < (*hits)[b].getObjectID().index; - } - return (*hits)[a].getLayer() < (*hits)[b].getLayer(); - }; - // indices contains the remaining hit indices that have not - // been assigned to a group yet - std::set indices(compare); - // set does not have a size yet, so cannot fill with iota - for (std::size_t i = 0; i < hits->size(); ++i) { - indices.insert(i); - } - // ensure no hits were dropped due to equivalency in set - if (hits->size() != indices.size()) { - error("equivalent hits were dropped: #hits {:d}, #indices {:d}", hits->size(), indices.size()); + // Group hit indices by system type + std::map> hits_by_system; + + auto* sys_field = m_idSpec.field("system"); + if (!sys_field) { + error("Field 'system' not found in IDSpec for readout {}", m_cfg.readout); + return; } - // Group neighbouring hits - std::vector> groups; - // because indices changes, the loop over indices requires some care: - // - we must use iterators instead of range-for - // - erase returns an incremented iterator and therefore acts as idx++ - // - when the set becomes empty on erase, idx is invalid and idx++ will be too - // (also applies to loop in bfs_group below) - for (auto idx = indices.begin(); idx != indices.end(); - indices.empty() ? idx = indices.end() : idx) { - - trace("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {}), energy = {}", - *idx, (*hits)[*idx].getLocal().x, (*hits)[*idx].getLocal().y, (*hits)[*idx].getLocal().z, - (*hits)[*idx].getPosition().x, (*hits)[*idx].getPosition().y, - (*hits)[*idx].getPosition().z, (*hits)[*idx].getEnergy()); - - // not energetic enough for cluster center, but could still be cluster hit - if ((*hits)[*idx].getEnergy() < minClusterCenterEdep) { - idx++; - continue; - } - - // create a new group, and group all the neighbouring hits - groups.emplace_back(std::list{*idx}); - bfs_group(*hits, indices, groups.back(), *idx); - - // wait with erasing until after bfs_group to ensure iterator is not invalidated in bfs_group - idx = indices.erase(idx); // takes role of idx++ + for (size_t i = 0; i < hits->size(); ++i) { + int sys = sys_field->value((*hits)[i].getCellID()); + hits_by_system[sys].push_back(i); } - debug("found {} potential clusters (groups of hits)", groups.size()); - for (std::size_t i = 0; i < groups.size(); ++i) { - debug("group {}: {} hits", i, groups[i].size()); - for (auto idx : groups[i]) { - const auto& hit = (*hits)[idx]; - trace(" hit {} -> energy = {:.6f}, layer = {}, sector = {}, local = ({:.2f}, {:.2f}, " - "{:.2f}), global = ({:.2f}, {:.2f}, {:.2f})", - idx, hit.getEnergy(), hit.getLayer(), hit.getSector(), hit.getLocal().x, - hit.getLocal().y, hit.getLocal().z, hit.getPosition().x, hit.getPosition().y, - hit.getPosition().z); + + // for cluster info + std::vector> all_clusters; + std::map> clusters_by_system; + + // Loop over systems + for (auto& [sys, hit_indices] : hits_by_system) { + + debug("Processing system {} with {} hits", sys, hit_indices.size()); + + // Sort hit indices (podio collections do not support std::sort) + auto compare = [&hits](const auto& a, const auto& b) { + // if !(a < b) and !(b < a), then a and b are equivalent + // and only one of them will be allowed in a set + const auto& ha = (*hits)[a]; + const auto& hb = (*hits)[b]; + + if (ha.getObjectID().collectionID != hb.getObjectID().collectionID) + return ha.getObjectID().collectionID < hb.getObjectID().collectionID; + if (ha.getLayer() != hb.getLayer()) + return ha.getLayer() < hb.getLayer(); + return ha.getObjectID().index < hb.getObjectID().index; + }; + + // indices contains the remaining hit indices that have not + // been assigned to a group yet + std::set indices(compare); + + // set does not have a size yet, so cannot fill with iota + for (auto i : hit_indices) { + indices.insert(i); + } + + if (hit_indices.size() != indices.size()) { + error("System {}: equivalent hits were dropped: #hits {:d}, #indices {:d}", sys, + hit_indices.size(), indices.size()); + for (auto i : hit_indices) { + const auto& h = (*hits)[i]; + error(" dropped hit {}: cellID=0x{:x}, layer={}, index={}", i, h.getCellID(), h.getLayer(), + h.getObjectID().index); + } + } + + // Group neighboring hits + std::vector> groups; + std::vector>> group_edges; //for is_neighbour + + // because indices changes, the loop over indices requires some care: + // - we must use iterators instead of range-for + // - erase returns an incremented iterator and therefore acts as idx++ + // - when the set becomes empty on erase, idx is invalid and idx++ will be too + // (also applies to loop in bfs_group below) + + for (auto idx = indices.begin(); idx != indices.end(); + indices.empty() ? idx = indices.end() : idx) { + + trace("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {}), energy = {}", + *idx, (*hits)[*idx].getLocal().x, (*hits)[*idx].getLocal().y, + (*hits)[*idx].getLocal().z, (*hits)[*idx].getPosition().x, + (*hits)[*idx].getPosition().y, (*hits)[*idx].getPosition().z, + (*hits)[*idx].getEnergy()); + + // not energetic enough for cluster center, but could still be cluster hit + if ((*hits)[*idx].getEnergy() < minClusterCenterEdep) { + idx++; + continue; + } + + // create a new group, and group all the neighbouring hits + groups.emplace_back(std::list{*idx}); + group_edges.emplace_back(); // create matching edges vector + bfs_group(*hits, indices, groups.back(), group_edges.back(), *idx); + + // wait with erasing until after bfs_group to ensure iterator is not invalidated in bfs_group + idx = indices.erase(idx); // takes role of idx++ + } + + debug("found {} potential clusters (groups of hits)", groups.size()); + for (std::size_t i = 0; i < groups.size(); ++i) { + debug("group {}: {} hits", i, groups[i].size()); + for (auto idx : groups[i]) { + const auto& hit = (*hits)[idx]; + double eta = edm4hep::utils::eta(hit.getPosition()); + double phi = edm4hep::utils::angleAzimuthal(hit.getPosition()); + + debug("hit {} , energy = {:.6f}|layer = {}, sector = {} local = ({:.2f}, {:.2f}, " + "{:.2f}),global = ({:.2f}, {:.2f}, {:.2f}) , eta = {:.2f}, phi = {:.2f} \n", + idx, hit.getEnergy(), hit.getLayer(), hit.getSector(), hit.getLocal().x, + hit.getLocal().y, hit.getLocal().z, hit.getPosition().x, hit.getPosition().y, + hit.getPosition().z, eta, phi); + } + } + + // form clusters + for (const auto& group : groups) { + if (group.size() < m_cfg.minClusterNhits) + continue; + + double energy = 0.; + for (auto idx : group) + energy += (*hits)[idx].getEnergy(); + if (energy < minClusterEdep) + continue; + + clusters_by_system[sys].push_back(all_clusters.size()); + all_clusters.emplace_back(group.begin(), group.end()); } } - // form clusters - for (const auto& group : groups) { - if (group.size() < m_cfg.minClusterNhits) { - continue; + // mergeCrossSystemClusters + std::vector> final_clusters; + std::set used_Img; + std::set used_ScFi; + + // Track cross-system neighbor pairs for visualization + std::set> cross_system_neighbor_pairs; + + bool have_Img = clusters_by_system.count(101); + bool have_ScFi = clusters_by_system.count(105); + + if (have_Img && have_ScFi) { + + const auto& Img = clusters_by_system.at(101); + const auto& ScFi = clusters_by_system.at(105); + + debug("Performing simple cross-system merging: {} Img clusters, {} ScFi clusters", Img.size(), + ScFi.size()); + + // Loop over Imaging clusters + for (size_t idx_Img : Img) { + + if (used_Img.count(idx_Img)) + continue; + + const auto& cl_Img = all_clusters[idx_Img]; + std::vector merged = cl_Img; + used_Img.insert(idx_Img); + + // Check against all ScFi clusters + for (size_t idx_ScFi : ScFi) { + + const auto& cl_ScFi = all_clusters[idx_ScFi]; + bool neighbor_found = false; + + for (auto h1 : cl_Img) { + for (auto h2 : cl_ScFi) { + if (cross_system_is_neighbour((*hits)[h1], (*hits)[h2])) { + debug(" Cross system neighbour: Img cluster {} hit {} <-> ScFi cluster {} hit {}", + idx_Img, h1, idx_ScFi, h2); + // Track this cross-system neighbor pair + size_t min_h = std::min(h1, h2); + size_t max_h = std::max(h1, h2); + cross_system_neighbor_pairs.insert({min_h, max_h}); + neighbor_found = true; + break; + } + } + if (neighbor_found) + break; + } + + if (neighbor_found) { + merged.insert(merged.end(), cl_ScFi.begin(), cl_ScFi.end()); + used_ScFi.insert(idx_ScFi); + + // merge any other Imaging clusters neighbour to ScFi + for (size_t idx2_Img : Img) { + if (used_Img.count(idx2_Img)) + continue; + + const auto& cl2_Img = all_clusters[idx2_Img]; + bool second_neighbour = false; + + for (auto h1 : cl2_Img) { + for (auto h2 : cl_ScFi) { + if (cross_system_is_neighbour((*hits)[h1], (*hits)[h2])) { + debug(" second Neighbour: Img cluster {} hit {} is neighbour to same ScFi " + "cluster {} via hit {}", + idx2_Img, h1, idx_ScFi, h2); + // Track this cross-system neighbor pair + size_t min_h = std::min(h1, h2); + size_t max_h = std::max(h1, h2); + cross_system_neighbor_pairs.insert({min_h, max_h}); + second_neighbour = true; + break; + } + } + if (second_neighbour) + break; + } + + if (second_neighbour) { + merged.insert(merged.end(), cl2_Img.begin(), cl2_Img.end()); + used_Img.insert(idx2_Img); + } + } + } + } + + final_clusters.push_back(std::move(merged)); } - double energy = 0.; - for (std::size_t idx : group) { - energy += (*hits)[idx].getEnergy(); + + // Add ScFi clusters that are not a part of any other cluster + for (size_t idx_ScFi : ScFi) { + if (!used_ScFi.count(idx_ScFi)) { + final_clusters.push_back(all_clusters[idx_ScFi]); + } } - if (energy < minClusterEdep) { - continue; + } else { + // If only one subsystem exists : no merging + final_clusters = all_clusters; + } + + debug("Outputting {} final proto-clusters", final_clusters.size()); + + for (std::size_t i = 0; i < final_clusters.size(); ++i) { + const auto& cluster = final_clusters[i]; + + debug("FINAL cluster {}: {} hits", i, cluster.size()); + + double total_energy = 0.0; + + for (auto idx : cluster) { + const auto& hit = (*hits)[idx]; + + int sys = sys_field->value(hit.getCellID()); + + double eta = edm4hep::utils::eta(hit.getPosition()); + double phi = edm4hep::utils::angleAzimuthal(hit.getPosition()); + + total_energy += hit.getEnergy(); + + debug(" hit {} , system = {} , energy = {:.6f} | layer = {}, sector = {} " + "local = ({:.2f}, {:.2f}, {:.2f}), " + "global = ({:.2f}, {:.2f}, {:.2f}), " + "eta = {:.2f}, phi = {:.2f}", + idx, sys, hit.getEnergy(), hit.getLayer(), hit.getSector(), hit.getLocal().x, + hit.getLocal().y, hit.getLocal().z, hit.getPosition().x, hit.getPosition().y, + hit.getPosition().z, eta, phi); } + + debug(" --> total cluster energy = {:.6f}", total_energy); + } + + // Write output + for (auto& cl : final_clusters) { auto pcl = proto->create(); - for (std::size_t idx : group) { + for (auto idx : cl) { pcl.addToHits((*hits)[idx]); pcl.addToWeights(1); } } } +bool ImagingTopoCluster::cross_system_is_neighbour(const edm4eic::CalorimeterHit& h1, + const edm4eic::CalorimeterHit& h2) const { + + // Get the "system" field from the ID specification + auto* sys_field = m_idSpec.field("system"); + if (!sys_field) { + error("Field 'system' not found in IDSpec for readout {}", m_cfg.readout); + return false; + } + // Extract system IDs for both hits + int sys1 = sys_field->value(h1.getCellID()); + int sys2 = sys_field->value(h2.getCellID()); + + // If different systems, allow Imaging(101) <-> ScFi(105) cross-linking + if (sys1 != sys2) { + if (h1.getSector() != h2.getSector()) { + return std::hypot((h1.getPosition().x - h2.getPosition().x), + (h1.getPosition().y - h2.getPosition().y), + (h1.getPosition().z - h2.getPosition().z)) <= cross_system_sectorDist; + } else { + return (std::abs(h1.getPosition().x - h2.getPosition().x) <= cross_system_DistXYZ[0]) && + (std::abs(h1.getPosition().y - h2.getPosition().y) <= cross_system_DistXYZ[1]) && + (std::abs(h1.getPosition().z - h2.getPosition().z) <= cross_system_DistXYZ[2]); + } + } + return false; +} + // helper function to group hits bool ImagingTopoCluster::is_neighbour(const edm4eic::CalorimeterHit& h1, const edm4eic::CalorimeterHit& h2) const { - // different sectors, simple distance check - if (h1.getSector() != h2.getSector()) { - return std::hypot((h1.getPosition().x - h2.getPosition().x), - (h1.getPosition().y - h2.getPosition().y), - (h1.getPosition().z - h2.getPosition().z)) <= sectorDist; + + // Get the "system" field from the ID specification + auto* sys_field = m_idSpec.field("system"); + if (!sys_field) { + error("Field 'system' not found in IDSpec for readout {}", m_cfg.readout); + return false; } + // Extract system IDs for both hits + int sys1 = sys_field->value(h1.getCellID()); + int sys2 = sys_field->value(h2.getCellID()); - // layer check - int ldiff = std::abs(h1.getLayer() - h2.getLayer()); - // same layer, check local positions - if (ldiff == 0) { - switch (m_cfg.sameLayerMode) { - case ImagingTopoClusterConfig::ELayerMode::xy: - return (std::abs(h1.getLocal().x - h2.getLocal().x) <= sameLayerDistXY[0]) && - (std::abs(h1.getLocal().y - h2.getLocal().y) <= sameLayerDistXY[1]); - - case ImagingTopoClusterConfig::ELayerMode::xyz: - return (std::abs(h1.getLocal().x - h2.getLocal().x) <= sameLayerDistXYZ[0]) && - (std::abs(h1.getLocal().y - h2.getLocal().y) <= sameLayerDistXYZ[1]) && - (std::abs(h1.getLocal().z - h2.getLocal().z) <= sameLayerDistXYZ[2]); - - case ImagingTopoClusterConfig::ELayerMode::etaphi: - return (std::abs(edm4hep::utils::eta(h1.getPosition()) - - edm4hep::utils::eta(h2.getPosition())) <= sameLayerDistEtaPhi[0]) && - (std::abs(edm4hep::utils::angleAzimuthal(h1.getPosition()) - - edm4hep::utils::angleAzimuthal(h2.getPosition())) <= sameLayerDistEtaPhi[1]); - - case ImagingTopoClusterConfig::ELayerMode::tz: { - // Layer mode 'tz' uses the average phi of the hits to define a rotated direction. The coordinate is a distance, not an angle. - auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + - edm4hep::utils::angleAzimuthal(h2.getPosition())); - auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); - auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); - auto h1_z = h1.getPosition().z; - auto h2_z = h2.getPosition().z; - - return (std::abs(h1_t - h2_t) <= sameLayerDistTZ[0]) && - (std::abs(h1_z - h2_z) <= sameLayerDistTZ[1]); - } - - default: - error("Unknown layer mode for same-layer clustering."); - return false; - } - } else if (ldiff <= m_cfg.neighbourLayersRange) { - switch (m_cfg.diffLayerMode) { - case eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi: - return (std::abs(edm4hep::utils::eta(h1.getPosition()) - - edm4hep::utils::eta(h2.getPosition())) <= diffLayerDistEtaPhi[0]) && - (std::abs(edm4hep::utils::angleAzimuthal(h1.getPosition()) - - edm4hep::utils::angleAzimuthal(h2.getPosition())) <= diffLayerDistEtaPhi[1]); - - case eicrecon::ImagingTopoClusterConfig::ELayerMode::xy: - // Here, the xy layer mode is based on global XY positions rather than local XY positions, and thus it only works for endcap detectors. - return (std::abs(h1.getPosition().x - h2.getPosition().x) <= diffLayerDistXY[0]) && - (std::abs(h1.getPosition().y - h2.getPosition().y) <= diffLayerDistXY[1]); - - case ImagingTopoClusterConfig::ELayerMode::xyz: - return (std::abs(h1.getPosition().x - h2.getPosition().x) <= diffLayerDistXYZ[0]) && - (std::abs(h1.getPosition().y - h2.getPosition().y) <= diffLayerDistXYZ[1]) && - (std::abs(h1.getPosition().z - h2.getPosition().z) <= diffLayerDistXYZ[2]); - - case eicrecon::ImagingTopoClusterConfig::ELayerMode::tz: { - auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + - edm4hep::utils::angleAzimuthal(h2.getPosition())); - auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); - auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); - auto h1_z = h1.getPosition().z; - auto h2_z = h2.getPosition().z; - - return (std::abs(h1_t - h2_t) <= diffLayerDistTZ[0]) && - (std::abs(h1_z - h2_z) <= diffLayerDistTZ[1]); - } - - default: - error("Unknown layer mode for different-layer clustering."); - return false; - } + if (sys1 != sys2) { + return false; } - // not in adjacent layers + // ScFi + if (sys1 == 105) { + // different sectors,distance check + if (h1.getSector() != h2.getSector()) { + return std::hypot((h1.getPosition().x - h2.getPosition().x), + (h1.getPosition().y - h2.getPosition().y), + (h1.getPosition().z - h2.getPosition().z)) <= ScFi_sectorDist; + } + + int ldiff = std::abs(h1.getLayer() - h2.getLayer()); + + double eta1 = edm4hep::utils::eta(h1.getPosition()); + double phi1 = edm4hep::utils::angleAzimuthal(h1.getPosition()); + double eta2 = edm4hep::utils::eta(h2.getPosition()); + double phi2 = edm4hep::utils::angleAzimuthal(h2.getPosition()); + + if (ldiff == 0) { + switch (m_cfg.ScFi_sameLayerMode) { + case ImagingTopoClusterConfig::ELayerMode::xy: + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= ScFi_sameLayerDistXY[0]) && + (std::abs(h1.getLocal().y - h2.getLocal().y) <= ScFi_sameLayerDistXY[1]); + case ImagingTopoClusterConfig::ELayerMode::xz: + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= ScFi_sameLayerDistXZ[0]) && + (std::abs(h1.getLocal().z - h2.getLocal().z) <= ScFi_sameLayerDistXZ[1]); + case ImagingTopoClusterConfig::ELayerMode::xyz: + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= ScFi_sameLayerDistXYZ[0]) && + (std::abs(h1.getLocal().y - h2.getLocal().y) <= ScFi_sameLayerDistXYZ[1]) && + (std::abs(h1.getLocal().z - h2.getLocal().z) <= ScFi_sameLayerDistXYZ[2]); + case ImagingTopoClusterConfig::ELayerMode::etaphi: + return (std::abs(eta1 - eta2) <= ScFi_sameLayerDistEtaPhi[0]) && + (std::abs(phi1 - phi2) <= ScFi_sameLayerDistEtaPhi[1]); + case ImagingTopoClusterConfig::ELayerMode::tz: { + auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + + edm4hep::utils::angleAzimuthal(h2.getPosition())); + auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); + auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); + auto h1_z = h1.getPosition().z; + auto h2_z = h2.getPosition().z; + return (std::abs(h1_t - h2_t) <= ScFi_sameLayerDistTZ[0]) && + (std::abs(h1_z - h2_z) <= ScFi_sameLayerDistTZ[1]); + } + default: + error("Unknown layer mode for same-layer clustering."); + return false; + } + } else if (ldiff <= m_cfg.ScFi_neighbourLayersRange) { + switch (m_cfg.ScFi_diffLayerMode) { + case eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi: + return (std::abs(eta1 - eta2) <= ScFi_diffLayerDistEtaPhi[0]) && + (std::abs(phi1 - phi2) <= ScFi_diffLayerDistEtaPhi[1]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::xy: + return (std::abs(h1.getPosition().x - h2.getPosition().x) <= ScFi_diffLayerDistXY[0]) && + (std::abs(h1.getPosition().y - h2.getPosition().y) <= ScFi_diffLayerDistXY[1]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::xz: + // return (std::abs(h1.getPosition().x - h2.getPosition().x) <= ScFi_diffLayerDistXZ[0]) && + // (std::abs(h1.getPosition().z - h2.getPosition().z) <= ScFi_diffLayerDistXZ[1]); + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= ScFi_diffLayerDistXZ[0]) && + (std::abs(h1.getLocal().z - h2.getLocal().z) <= ScFi_diffLayerDistXZ[1]); + case ImagingTopoClusterConfig::ELayerMode::xyz: + // return (std::abs(h1.getPosition().x - h2.getPosition().x) <= ScFi_diffLayerDistXYZ[0]) && + // (std::abs(h1.getPosition().y - h2.getPosition().y) <= ScFi_diffLayerDistXYZ[1]) && + // (std::abs(h1.getPosition().z - h2.getPosition().z) <= ScFi_diffLayerDistXYZ[2]); + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= ScFi_diffLayerDistXYZ[0]) && + (std::abs(h1.getLocal().y - h2.getLocal().y) <= ScFi_diffLayerDistXYZ[1]) && + (std::abs(h1.getLocal().z - h2.getLocal().z) <= ScFi_diffLayerDistXYZ[2]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::tz: { + auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + + edm4hep::utils::angleAzimuthal(h2.getPosition())); + auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); + auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); + auto h1_z = h1.getPosition().z; + auto h2_z = h2.getPosition().z; + return (std::abs(h1_t - h2_t) <= ScFi_diffLayerDistTZ[0]) && + (std::abs(h1_z - h2_z) <= ScFi_diffLayerDistTZ[1]); + } + default: + error("Unknown layer mode for same-layer clustering."); + return false; + } + } + } // sys1 == 105 + + // Imaging + else if (sys1 == 101) { + // different sectors, simple distance check + if (h1.getSector() != h2.getSector()) { + return std::hypot((h1.getPosition().x - h2.getPosition().x), + (h1.getPosition().y - h2.getPosition().y), + (h1.getPosition().z - h2.getPosition().z)) <= Img_sectorDist; + } + + int ldiff = std::abs(h1.getLayer() - h2.getLayer()); + + double eta1 = edm4hep::utils::eta(h1.getPosition()); + double phi1 = edm4hep::utils::angleAzimuthal(h1.getPosition()); + double eta2 = edm4hep::utils::eta(h2.getPosition()); + double phi2 = edm4hep::utils::angleAzimuthal(h2.getPosition()); + + if (ldiff == 0) { + switch (m_cfg.Img_sameLayerMode) { + case ImagingTopoClusterConfig::ELayerMode::xy: + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= Img_sameLayerDistXY[0]) && + (std::abs(h1.getLocal().y - h2.getLocal().y) <= Img_sameLayerDistXY[1]); + case ImagingTopoClusterConfig::ELayerMode::xyz: + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= Img_sameLayerDistXYZ[0]) && + (std::abs(h1.getLocal().y - h2.getLocal().y) <= Img_sameLayerDistXYZ[1]) && + (std::abs(h1.getLocal().z - h2.getLocal().z) <= Img_sameLayerDistXYZ[2]); + case ImagingTopoClusterConfig::ELayerMode::etaphi: + return (std::abs(eta1 - eta2) <= Img_sameLayerDistEtaPhi[0]) && + (std::abs(phi1 - phi2) <= Img_sameLayerDistEtaPhi[1]); + case ImagingTopoClusterConfig::ELayerMode::tz: { + auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + + edm4hep::utils::angleAzimuthal(h2.getPosition())); + auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); + auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); + auto h1_z = h1.getPosition().z; + auto h2_z = h2.getPosition().z; + return (std::abs(h1_t - h2_t) <= Img_sameLayerDistTZ[0]) && + (std::abs(h1_z - h2_z) <= Img_sameLayerDistTZ[1]); + } + default: + error("Unknown layer mode for same-layer clustering."); + return false; + } + } else if (ldiff <= m_cfg.neighbourLayersRange) { + switch (m_cfg.Img_diffLayerMode) { + case eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi: + return (std::abs(eta1 - eta2) <= Img_diffLayerDistEtaPhi[0]) && + (std::abs(phi1 - phi2) <= Img_diffLayerDistEtaPhi[1]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::xy: + return (std::abs(h1.getPosition().x - h2.getPosition().x) <= Img_diffLayerDistXY[0]) && + (std::abs(h1.getPosition().y - h2.getPosition().y) <= Img_diffLayerDistXY[1]); + case ImagingTopoClusterConfig::ELayerMode::xyz: + return (std::abs(h1.getPosition().x - h2.getPosition().x) <= Img_diffLayerDistXYZ[0]) && + (std::abs(h1.getPosition().y - h2.getPosition().y) <= Img_diffLayerDistXYZ[1]) && + (std::abs(h1.getPosition().z - h2.getPosition().z) <= Img_diffLayerDistXYZ[2]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::tz: { + auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + + edm4hep::utils::angleAzimuthal(h2.getPosition())); + auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); + auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); + auto h1_z = h1.getPosition().z; + auto h2_z = h2.getPosition().z; + return (std::abs(h1_t - h2_t) <= Img_diffLayerDistTZ[0]) && + (std::abs(h1_z - h2_z) <= Img_diffLayerDistTZ[1]); + } + default: + error("Hits are not neighbors: idx1 = {}, idx2 = {}, layer1 = {}, layer2 = {}, ldiff = {}", + h1.getObjectID().index, h2.getObjectID().index, h1.getLayer(), h2.getLayer(), ldiff); + return false; + } + } + } // sys1 == 101 + + // any other system other than ScFi and Img + else { + + // different sectors, simple distance check + if (h1.getSector() != h2.getSector()) { + return std::hypot((h1.getPosition().x - h2.getPosition().x), + (h1.getPosition().y - h2.getPosition().y), + (h1.getPosition().z - h2.getPosition().z)) <= sectorDist; + } + + // layer check + int ldiff = std::abs(h1.getLayer() - h2.getLayer()); + if (ldiff == 0) { + switch (m_cfg.sameLayerMode) { + case ImagingTopoClusterConfig::ELayerMode::xy: + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= sameLayerDistXY[0]) && + (std::abs(h1.getLocal().y - h2.getLocal().y) <= sameLayerDistXY[1]); + case ImagingTopoClusterConfig::ELayerMode::xyz: + return (std::abs(h1.getLocal().x - h2.getLocal().x) <= sameLayerDistXYZ[0]) && + (std::abs(h1.getLocal().y - h2.getLocal().y) <= sameLayerDistXYZ[1]) && + (std::abs(h1.getLocal().z - h2.getLocal().z) <= sameLayerDistXYZ[2]); + case ImagingTopoClusterConfig::ELayerMode::etaphi: + return (std::abs(edm4hep::utils::eta(h1.getPosition()) - + edm4hep::utils::eta(h2.getPosition())) <= sameLayerDistEtaPhi[0]) && + (std::abs(edm4hep::utils::angleAzimuthal(h1.getPosition()) - + edm4hep::utils::angleAzimuthal(h2.getPosition())) <= + sameLayerDistEtaPhi[1]); + case ImagingTopoClusterConfig::ELayerMode::tz: { + auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + + edm4hep::utils::angleAzimuthal(h2.getPosition())); + auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); + auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); + auto h1_z = h1.getPosition().z; + auto h2_z = h2.getPosition().z; + return (std::abs(h1_t - h2_t) <= sameLayerDistTZ[0]) && + (std::abs(h1_z - h2_z) <= sameLayerDistTZ[1]); + } + default: + error("Unknown layer mode for same-layer clustering."); + return false; + } + } else if (ldiff <= m_cfg.neighbourLayersRange) { + switch (m_cfg.diffLayerMode) { + case eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi: + return (std::abs(edm4hep::utils::eta(h1.getPosition()) - + edm4hep::utils::eta(h2.getPosition())) <= diffLayerDistEtaPhi[0]) && + (std::abs(edm4hep::utils::angleAzimuthal(h1.getPosition()) - + edm4hep::utils::angleAzimuthal(h2.getPosition())) <= + diffLayerDistEtaPhi[1]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::xy: + return (std::abs(h1.getPosition().x - h2.getPosition().x) <= diffLayerDistXY[0]) && + (std::abs(h1.getPosition().y - h2.getPosition().y) <= diffLayerDistXY[1]); + case ImagingTopoClusterConfig::ELayerMode::xyz: + return (std::abs(h1.getPosition().x - h2.getPosition().x) <= diffLayerDistXYZ[0]) && + (std::abs(h1.getPosition().y - h2.getPosition().y) <= diffLayerDistXYZ[1]) && + (std::abs(h1.getPosition().z - h2.getPosition().z) <= diffLayerDistXYZ[2]); + case eicrecon::ImagingTopoClusterConfig::ELayerMode::tz: { + auto phi = 0.5 * (edm4hep::utils::angleAzimuthal(h1.getPosition()) + + edm4hep::utils::angleAzimuthal(h2.getPosition())); + auto h1_t = (h1.getPosition().x * sin(phi)) - (h1.getPosition().y * cos(phi)); + auto h2_t = (h2.getPosition().x * sin(phi)) - (h2.getPosition().y * cos(phi)); + auto h1_z = h1.getPosition().z; + auto h2_z = h2.getPosition().z; + return (std::abs(h1_t - h2_t) <= diffLayerDistTZ[0]) && + (std::abs(h1_z - h2_z) <= diffLayerDistTZ[1]); + } + default: + error("Unknown layer mode for different-layer clustering."); + return false; + } + } + } // any other system + return false; } diff --git a/src/algorithms/calorimetry/ImagingTopoCluster.h b/src/algorithms/calorimetry/ImagingTopoCluster.h index 1d7654b09b..78da4d33ec 100644 --- a/src/algorithms/calorimetry/ImagingTopoCluster.h +++ b/src/algorithms/calorimetry/ImagingTopoCluster.h @@ -18,18 +18,26 @@ * hits while keeping iterators valid * */ + #pragma once +#include +#include +#include #include +#include // Event Model related classes #include #include #include #include +#include +#include #include #include #include #include +#include #include "ImagingTopoClusterConfig.h" #include "algorithms/interfaces/WithPodConfig.h" @@ -55,13 +63,37 @@ class ImagingTopoCluster : public ImagingTopoClusterAlgorithm, // unitless counterparts of the input parameters std::array sameLayerDistXY{0, 0}; std::array diffLayerDistXY{0, 0}; + std::array ScFi_sameLayerDistXY{0, 0}; + std::array ScFi_diffLayerDistXY{0, 0}; + std::array ScFi_sameLayerDistXZ{0, 0}; + std::array ScFi_diffLayerDistXZ{0, 0}; + std::array Img_sameLayerDistXY{0, 0}; + std::array Img_diffLayerDistXY{0, 0}; std::array sameLayerDistXYZ{0, 0, 0}; std::array diffLayerDistXYZ{0, 0, 0}; + std::array ScFi_sameLayerDistXYZ{0, 0, 0}; + std::array ScFi_diffLayerDistXYZ{0, 0, 0}; + std::array Img_sameLayerDistXYZ{0, 0, 0}; + std::array Img_diffLayerDistXYZ{0, 0, 0}; std::array sameLayerDistEtaPhi{0, 0}; std::array diffLayerDistEtaPhi{0, 0}; std::array sameLayerDistTZ{0, 0}; std::array diffLayerDistTZ{0, 0}; + std::array ScFi_sameLayerDistEtaPhi{0, 0}; + std::array ScFi_diffLayerDistEtaPhi{0, 0}; + std::array ScFi_sameLayerDistTZ{0, 0}; + std::array ScFi_diffLayerDistTZ{0, 0}; + std::array Img_sameLayerDistEtaPhi{0, 0}; + std::array Img_diffLayerDistEtaPhi{0, 0}; + std::array Img_sameLayerDistTZ{0, 0}; + std::array Img_diffLayerDistTZ{0, 0}; + + std::array cross_system_DistXYZ{0, 0, 0}; + double sectorDist{0}; + double cross_system_sectorDist{0}; + double ScFi_sectorDist{0}; + double Img_sectorDist{0}; double minClusterHitEdep{0}; double minClusterCenterEdep{0}; double minClusterEdep{0}; @@ -70,16 +102,40 @@ class ImagingTopoCluster : public ImagingTopoClusterAlgorithm, void init(); void process(const Input& input, const Output& output) const final; + // based on system Id and across the system neighbouring + private: // helper function to group hits + + // std::vector> mergeCrossSystemClusters(const std::map>& clusters_by_system, + // const std::vector>& all_clusters, + // const edm4eic::CalorimeterHitCollection& hits) const; + + bool cross_system_is_neighbour(const edm4eic::CalorimeterHit& h1, + const edm4eic::CalorimeterHit& h2) const; bool is_neighbour(const edm4eic::CalorimeterHit& h1, const edm4eic::CalorimeterHit& h2) const; + // Pointer to the geometry service + dd4hep::IDDescriptor m_idSpec; + + const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()}; + // grouping function with Breadth-First Search // note: template to allow Compare only known in local scope of caller template void bfs_group(const edm4eic::CalorimeterHitCollection& hits, std::set& indices, std::list& group, - const std::size_t idx) const { + std::vector>& edges, const std::size_t idx) const { + + auto* sys_field = m_idSpec.field("system"); + if (!sys_field) { + error("Field 'system' not found in IDSpec for BFS grouping"); + return; + } + + int sys = sys_field->value(hits[idx].getCellID()); + // debug("Starting BFS for hit {} in system {}", idx, sys); // loop over group as it grows, until the end is stable and we reach it for (auto idx1 = group.begin(); idx1 != group.end(); ++idx1) { @@ -94,6 +150,15 @@ class ImagingTopoCluster : public ImagingTopoClusterAlgorithm, continue; } + // debug("Checking neighbor for hit {} in system {}", *idx2, sys_field->value(hits[*idx2].getCellID())); + + // skip hits from other system + if (sys_field->value(hits[*idx2].getCellID()) != sys) { + debug(" Skipping hit {}: different system", *idx2); + ++idx2; + continue; + } + // skip rest of list of hits when we're past relevant layers //if (hits[*idx2].getLayer() - hits[*idx1].getLayer() > m_cfg.neighbourLayersRange) { // break; @@ -106,6 +171,7 @@ class ImagingTopoCluster : public ImagingTopoClusterAlgorithm, } if (is_neighbour(hits[*idx1], hits[*idx2])) { + edges.emplace_back(*idx1, *idx2); group.push_back(*idx2); idx2 = indices.erase(idx2); // takes role of idx2++ } else { diff --git a/src/algorithms/calorimetry/ImagingTopoClusterConfig.h b/src/algorithms/calorimetry/ImagingTopoClusterConfig.h index 3a8fb62269..5468e7e407 100644 --- a/src/algorithms/calorimetry/ImagingTopoClusterConfig.h +++ b/src/algorithms/calorimetry/ImagingTopoClusterConfig.h @@ -13,37 +13,84 @@ namespace eicrecon { struct ImagingTopoClusterConfig { + std::string readout = ""; + // maximum difference in layer numbers that can be considered as neighbours - int neighbourLayersRange = 1; - // maximum distance of global (x, y) to be considered as neighbors at same layers (if sameLayerMode==xy) - std::vector> sameLayerDistXY = {1.0 * dd4hep::mm, - 1.0 * dd4hep::mm}; - // maximum distance of local (x, y,z) to be considered as neighbors at same layers (if sameLayerMode==xyz) - std::vector> sameLayerDistXYZ = { - 1.0 * dd4hep::mm, 1.0 * dd4hep::mm, 20.0 * dd4hep::mm}; - // maximum distance of global (eta, phi) to be considered as neighbors at same layers (if sameLayerMode==etaphi) - std::vector sameLayerDistEtaPhi = {0.01, 0.01}; - // maximum distance of global (t, z) to be considered as neighbors at same layers (if sameLayerMode==tz) - std::vector sameLayerDistTZ = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; - // maximum distance of global (x, y) to be considered as neighbors at different layers (if diffLayerMode==xy) - std::vector> diffLayerDistXY = {1.0 * dd4hep::mm, - 1.0 * dd4hep::mm}; - // maximum distance of global (x, y, z) to be considered as neighbors at different layers (if diffLayerMode==xyz) - std::vector> diffLayerDistXYZ = { - 1.0 * dd4hep::mm, 1.0 * dd4hep::mm, 20.0 * dd4hep::mm}; - // maximum distance of global (eta, phi) to be considered as neighbors at different layers (if diffLayerMode==etaphi) - std::vector diffLayerDistEtaPhi = {0.01, 0.01}; - // maximum distance of global (t, z) to be considered as neighbors at different layers (if diffLayerMode==tz) - std::vector diffLayerDistTZ = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; + int neighbourLayersRange = 1; + int ScFi_neighbourLayersRange = 3; + // maximum distance of local (x, y) to be considered as neighbors at same layers (if samelayerMode==xy) + std::vector> sameLayerDistXY = {1.0 * dd4hep::mm, + 1.0 * dd4hep::mm}; + std::vector> ScFi_sameLayerDistXY = {1.0 * dd4hep::mm, + 1.0 * dd4hep::mm}; + std::vector> Img_sameLayerDistXY = {1.0 * dd4hep::mm, + 1.0 * dd4hep::mm}; + std::vector ScFi_sameLayerDistXZ = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; + std::vector ScFi_diffLayerDistXZ = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; + + // maximum distance of local (x, y,z) to be considered as neighbors at same layers (if samelayerMode==xyz) + std::vector sameLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, 40.0 * dd4hep::mm}; + std::vector ScFi_sameLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, + 40.0 * dd4hep::mm}; + std::vector Img_sameLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, + 40.0 * dd4hep::mm}; + + // maximum distance of global (eta, phi) to be considered as neighbors at same layers (if samelayerMode==etaphi) + std::vector sameLayerDistEtaPhi = {0.01, 0.01}; + std::vector ScFi_sameLayerDistEtaPhi = {0.01, 0.01}; + std::vector Img_sameLayerDistEtaPhi = {0.01, 0.01}; + + // maximum distance of global (t, z) to be considered as neighbors at same layers (if samelayerMode==tz) + std::vector sameLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}; + std::vector ScFi_sameLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}; + std::vector Img_sameLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}; + + // maximum distance of global (x, y) to be considered as neighbors at different layers (if difflayerMode==xy) + std::vector> diffLayerDistXY = {1.0 * dd4hep::mm, + 1.0 * dd4hep::mm}; + std::vector> ScFi_diffLayerDistXY = {1.0 * dd4hep::mm, + 1.0 * dd4hep::mm}; + std::vector> Img_diffLayerDistXY = {1.0 * dd4hep::mm, + 1.0 * dd4hep::mm}; + + // maximum distance of global (x, y,z) to be considered as neighbors at different layers (if difflayerMode==xyz) + std::vector diffLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, 40.0 * dd4hep::mm}; + std::vector ScFi_diffLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, + 40.0 * dd4hep::mm}; + std::vector Img_diffLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, + 40.0 * dd4hep::mm}; + + // maximum distance of global (eta, phi) to be considered as neighbors at different layers (if difflayerMode==etaphi) + std::vector diffLayerDistEtaPhi = {0.01, 0.01}; + std::vector ScFi_diffLayerDistEtaPhi = {0.01, 0.01}; + std::vector Img_diffLayerDistEtaPhi = {0.01, 0.01}; + + // maximum distance of global (t, z) to be considered as neighbors at different layers (if difflayerMode==tz) + std::vector diffLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}; + std::vector ScFi_diffLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}; + std::vector Img_diffLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}; + + std::vector cross_system_DistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, + 40.0 * dd4hep::mm}; + // Layermodes - enum class ELayerMode { etaphi = 0, xy = 1, tz = 2, xyz = 3 }; + enum class ELayerMode { etaphi = 0, xy = 1, xz = 2, xyz = 3, tz = 4 }; + // determines how neighbors are determined for hits in same layers (using either eta and phi, or x and y) - ELayerMode sameLayerMode = ELayerMode::xy; // for ldiff =0 + ELayerMode sameLayerMode = ELayerMode::xy; // for ldiff =0 + ELayerMode ScFi_sameLayerMode = ELayerMode::xyz; + ELayerMode Img_sameLayerMode = ELayerMode::tz; + // determines how neighbors are determined for hits in different layers (using either eta and phi, or x and y) - ELayerMode diffLayerMode = ELayerMode::xy; // for ldiff <= neighbourLayersRange + ELayerMode diffLayerMode = ELayerMode::etaphi; // for ldiff <= neighbourLayersRange + ELayerMode ScFi_diffLayerMode = ELayerMode::xyz; + ELayerMode Img_diffLayerMode = ELayerMode::etaphi; // maximum global distance to be considered as neighbors in different sectors - double sectorDist = 1.0 * dd4hep::cm; + double sectorDist = 3.0 * dd4hep::cm; + double cross_system_sectorDist = 5.0 * dd4hep::cm; + double ScFi_sectorDist = 5.0 * dd4hep::cm; + double Img_sectorDist = 3.0 * dd4hep::cm; // minimum hit energy to participate clustering double minClusterHitEdep = 0.; @@ -63,10 +110,12 @@ std::istream& operator>>(std::istream& in, ImagingTopoClusterConfig::ELayerMode& layerMode = ImagingTopoClusterConfig::ELayerMode::etaphi; } else if (s == "xy" or s == "1") { layerMode = ImagingTopoClusterConfig::ELayerMode::xy; - } else if (s == "tz" or s == "2") { - layerMode = ImagingTopoClusterConfig::ELayerMode::tz; + } else if (s == "xz" or s == "2") { + layerMode = ImagingTopoClusterConfig::ELayerMode::xz; } else if (s == "xyz" or s == "3") { layerMode = ImagingTopoClusterConfig::ELayerMode::xyz; + } else if (s == "tz" or s == "4") { + layerMode = ImagingTopoClusterConfig::ELayerMode::tz; } else { in.setstate(std::ios::failbit); // Set the fail bit if the input is not valid } @@ -81,6 +130,9 @@ std::ostream& operator<<(std::ostream& out, const ImagingTopoClusterConfig::ELay case ImagingTopoClusterConfig::ELayerMode::xy: out << "xy"; break; + case ImagingTopoClusterConfig::ELayerMode::xz: + out << "xz"; + break; case ImagingTopoClusterConfig::ELayerMode::xyz: out << "xyz"; break; diff --git a/src/detectors/BEMC/BEMC.cc b/src/detectors/BEMC/BEMC.cc index 0628a86917..a79d21f1bd 100644 --- a/src/detectors/BEMC/BEMC.cc +++ b/src/detectors/BEMC/BEMC.cc @@ -5,9 +5,14 @@ #include #include #include +#include +#include +#include #include #include +#include #include +#include #include #include #include @@ -38,6 +43,7 @@ #include "factories/digi/PulseCombiner_factory.h" #include "factories/digi/PulseGeneration_factory.h" #include "factories/digi/PulseNoise_factory.h" +#include "factories/meta/CollectionCollector_factory.h" extern "C" { void InitPlugin(JApplication* app) { @@ -267,6 +273,7 @@ void InitPlugin(JApplication* app) { .layerField = "layer", .sectorField = "sector", .localDetFields = {"system", "sector"}, + // here we want to use grid center position (XY) but keeps the z information from fiber-segment // TODO: a more realistic way to get z is to reconstruct it from timing .maskPos = "xy", @@ -306,7 +313,11 @@ void InitPlugin(JApplication* app) { {"EcalBarrelScFiClustersWithoutShapes", // edm4eic::Cluster "EcalBarrelScFiClusterLinksWithoutShapes", "EcalBarrelScFiClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation - {.energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, .enableEtaBounds = false}, + {.readout = "EcalBarrelScFiHits", + .energyWeight = "log", + .sampFrac = 1.0, + .logWeightBase = 6.2, + .enableEtaBounds = false}, app // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( @@ -315,15 +326,20 @@ void InitPlugin(JApplication* app) { {"EcalBarrelScFiClusters", "EcalBarrelScFiClusterLinks", "EcalBarrelScFiClusterAssociations"}, {.longitudinalShowerInfoAvailable = true, .energyWeight = "log", .logWeightBase = 6.2}, app)); - // Imaging TopoClustering on ScFi + //-------------------------------------------------------------------------------------------------------- + // TopoClustering on ScFi + //-------------------------------------------------------------------------------------------------------- + app->Add(new JOmniFactoryGeneratorT( "EcalBarrelScFiProtoTopoClusters", {"EcalBarrelScFiRecHits"}, {"EcalBarrelScFiProtoTopoClusters"}, { + .readout = "EcalBarrelScFiHits", .neighbourLayersRange = 2, // # id diff for adjacent layer .sameLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, - 40.0 * dd4hep::mm}, // # same layer - .diffLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, 40.0 * dd4hep::mm}, + 80.0 * dd4hep::mm}, // # same layer + .diffLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, + 80.0 * dd4hep::mm}, // # diff layer .sameLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xyz, .diffLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xyz, .sectorDist = 5.0 * dd4hep::cm, @@ -338,13 +354,17 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "EcalBarrelScFiTopoClustersWithoutShapes", - {"EcalBarrelScFiProtoTopoClusters", // edm4eic::ProtoClusterCollection - "EcalBarrelScFiRawHitLinks", // edm4eic::MCRecoCalorimeterHitLink + {"EcalBarrelScFiProtoTopoClusters", // edm4eic::ProtoClusterCollection + "EcalBarrelScFiRawHitLinks", "EcalBarrelScFiRawHitAssociations"}, // edm4eic::MCRecoCalorimeterHitAssociation {"EcalBarrelScFiTopoClustersWithoutShapes", // edm4eic::Cluster "EcalBarrelScFiTopoClusterLinksWithoutShapes", "EcalBarrelScFiTopoClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation - {.energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, .enableEtaBounds = false}, + {.readout = "EcalBarrelScFiHits", + .energyWeight = "log", + .sampFrac = 1.0, + .logWeightBase = 6.2, + .enableEtaBounds = false}, app // TODO: Remove me once fixed )); app->Add(new JOmniFactoryGeneratorT( @@ -414,12 +434,14 @@ void InitPlugin(JApplication* app) { "EcalBarrelImagingProtoClusters", {"EcalBarrelImagingRecHits"}, {"EcalBarrelImagingProtoClusters"}, { + .readout = "EcalBarrelImagingHits", .neighbourLayersRange = 2, // # id diff for adjacent layer - .sameLayerDistTZ = {2.0 * dd4hep::mm, 2 * dd4hep::mm}, // # same layer - .diffLayerDistEtaPhi = {10 * dd4hep::mrad, 10 * dd4hep::mrad}, // # adjacent layer + .sameLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}, // # same layer + .diffLayerDistEtaPhi = {10.0 * dd4hep::mrad, 10.0 * dd4hep::mrad}, // # adjacent layer .sameLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::tz, .diffLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi, .sectorDist = 3.0 * dd4hep::cm, + .minClusterHitEdep = 0, .minClusterCenterEdep = 0, .minClusterEdep = 100 * dd4hep::MeV, @@ -439,6 +461,7 @@ void InitPlugin(JApplication* app) { }, app // TODO: Remove me once fixed )); + app->Add(new JOmniFactoryGeneratorT( "EcalBarrelImagingClusters", {"EcalBarrelImagingClustersWithoutShapes", @@ -447,6 +470,82 @@ void InitPlugin(JApplication* app) { "EcalBarrelImagingClusterAssociations"}, {.longitudinalShowerInfoAvailable = false, .energyWeight = "log", .logWeightBase = 6.2}, app)); + + //-------------------------------------------------------------------------------------------------------- + // Combined Collection of ScFi and Imaging Hits + //-------------------------------------------------------------------------------------------------------- + + app->Add(new JOmniFactoryGeneratorT>( + "EcalBarrelTopoRecHits", {"EcalBarrelScFiRecHits", "EcalBarrelImagingRecHits"}, + {"EcalBarrelTopoRecHits"}, app)); + + app->Add(new JOmniFactoryGeneratorT< + CollectionCollector_factory>( + "EcalBarrelTopoRawHitAssociations", + {"EcalBarrelImagingRawHitAssociations", "EcalBarrelScFiRawHitAssociations"}, + {"EcalBarrelTopoRawHitAssociations"}, app)); + + app->Add( + new JOmniFactoryGeneratorT>( + "EcalBarrelTopoRawHitLinks", + {"EcalBarrelImagingRawHitLinks", "EcalBarrelScFiRawHitLinks"}, + {"EcalBarrelTopoRawHitLinks"}, app)); + + //-------------------------------------------------------------------------------------------------------- + // TopoClustering on Combined clustering ScFi and Imaging Hits + // (SystemID based Clustering) + //-------------------------------------------------------------------------------------------------------- + + app->Add(new JOmniFactoryGeneratorT( + "EcalBarrelTopoProtoClusters", {"EcalBarrelTopoRecHits"}, {"EcalBarrelTopoProtoClusters"}, + { + .readout = + "EcalBarrelScFiHits", // we need a readout name to connect the cell ID to a system ID only, so it does not matter that this is more specific than the hits we feed here + .neighbourLayersRange = 2, // # id diff for adjacent layer + .ScFi_sameLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, + 40.0 * dd4hep::mm}, // # same layer + .Img_sameLayerDistTZ = {2.0 * dd4hep::mm, 2.0 * dd4hep::mm}, + .ScFi_diffLayerDistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, 80.0 * dd4hep::mm}, + .Img_diffLayerDistEtaPhi = {10.0 * dd4hep::mrad, 10.0 * dd4hep::mrad}, + .cross_system_DistXYZ = {80.0 * dd4hep::mm, 80.0 * dd4hep::mm, 80.0 * dd4hep::mm}, + .ScFi_sameLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xyz, + .Img_sameLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::tz, + .ScFi_diffLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xyz, + .Img_diffLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::etaphi, + .cross_system_sectorDist = 5.0 * dd4hep::cm, + .ScFi_sectorDist = 5.0 * dd4hep::cm, + .Img_sectorDist = 3.0 * dd4hep::cm, + .minClusterHitEdep = 0, + .minClusterCenterEdep = 0, + .minClusterEdep = 100 * dd4hep::MeV, + .minClusterNhits = 10, + }, + app // TODO: Remove me once fixed + )); + + app->Add(new JOmniFactoryGeneratorT( + "EcalBarrelTopoClustersWithoutShapes", + {"EcalBarrelTopoProtoClusters", // edm4eic::ProtoClusterCollection + "EcalBarrelTopoRawHitLinks", + "EcalBarrelTopoRawHitAssociations"}, // edm4eic::MCRecoCalorimeterHitAssociation + {"EcalBarrelTopoClustersWithoutShapes", // edm4eic::Cluster + "EcalBarrelTopoClusterLinksWithoutShapes", + "EcalBarrelTopoClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation + { + .readout = "EcalBarrelScFiHits", + .energyWeight = "log", + .sampFrac = 1.0, + .logWeightBase = 6.2, + .enableEtaBounds = false, + }, + app // TODO: Remove me once fixed + )); + app->Add(new JOmniFactoryGeneratorT( + "EcalBarrelTopoClusters", + {"EcalBarrelTopoClustersWithoutShapes", "EcalBarrelTopoClusterAssociationsWithoutShapes"}, + {"EcalBarrelTopoClusters", "EcalBarrelTopoClusterLinks", "EcalBarrelTopoClusterAssociations"}, + {.longitudinalShowerInfoAvailable = true, .energyWeight = "log", .logWeightBase = 6.2}, app)); + app->Add(new JOmniFactoryGeneratorT( "EcalBarrelClustersWithoutShapes", {"EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations", "EcalBarrelImagingClusters", diff --git a/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc index f9c6dea44e..8499c22879 100644 --- a/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc +++ b/src/tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc @@ -45,6 +45,7 @@ TEST_CASE("the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]") { logger->set_level(spdlog::level::trace); CalorimeterClusterRecoCoGConfig cfg; + cfg.readout = "MockCalorimeterHits"; cfg.energyWeight = "log"; cfg.sampFrac = 0.0203; cfg.logWeightBaseCoeffs = {5.0, 0.65, 0.31}; diff --git a/src/tests/algorithms_test/calorimetry_ImagingTopoCluster.cc b/src/tests/algorithms_test/calorimetry_ImagingTopoCluster.cc index 75ea739908..03ac9375ed 100644 --- a/src/tests/algorithms_test/calorimetry_ImagingTopoCluster.cc +++ b/src/tests/algorithms_test/calorimetry_ImagingTopoCluster.cc @@ -33,11 +33,14 @@ TEST_CASE("the clustering algorithm runs", "[ImagingTopoCluster]") { logger->set_level(spdlog::level::trace); ImagingTopoClusterConfig cfg; + cfg.readout = "MockCalorimeterHits"; cfg.sameLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xy; cfg.minClusterHitEdep = 0. * dd4hep::GeV; cfg.minClusterCenterEdep = 0. * dd4hep::GeV; - cfg.sameLayerDistXY = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; //mm - cfg.diffLayerDistXY = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; //mm + cfg.sameLayerMode = ImagingTopoClusterConfig::ELayerMode::xy; + cfg.sameLayerDistXY = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; + cfg.diffLayerMode = ImagingTopoClusterConfig::ELayerMode::xy; + cfg.diffLayerDistXY = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; cfg.minClusterEdep = 9 * dd4hep::MeV; // minimum number of hits (to save this cluster) cfg.minClusterNhits = 1; diff --git a/src/tests/algorithms_test/particle_flow_TrackProtoClusterMatchPromoter.cc b/src/tests/algorithms_test/particle_flow_TrackProtoClusterMatchPromoter.cc index 0c4660fd8c..d4a548a1b8 100644 --- a/src/tests/algorithms_test/particle_flow_TrackProtoClusterMatchPromoter.cc +++ b/src/tests/algorithms_test/particle_flow_TrackProtoClusterMatchPromoter.cc @@ -127,6 +127,7 @@ TEST_CASE("the TrackProtoClusterMatchPromoter algorithm runs", "[TrackProtoClust // configure reco algorithm to match EEEMCAL eicrecon::CalorimeterClusterRecoCoGConfig cfg_reco; + cfg_reco.readout = "MockCalorimeterHits"; cfg_reco.energyWeight = "log"; cfg_reco.sampFrac = 1.0; cfg_reco.logWeightBase = 3.6;