diff --git a/src/algorithms/calorimetry/ImagingTopoCluster.cc b/src/algorithms/calorimetry/ImagingTopoCluster.cc index 8031c7c6a6..3125590768 100644 --- a/src/algorithms/calorimetry/ImagingTopoCluster.cc +++ b/src/algorithms/calorimetry/ImagingTopoCluster.cc @@ -63,6 +63,12 @@ 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; + sameLayerDistXYZ[0] = m_cfg.sameLayerDistXYZ[0] / dd4hep::mm; + sameLayerDistXYZ[1] = m_cfg.sameLayerDistXYZ[1] / dd4hep::mm; + sameLayerDistXYZ[2] = m_cfg.sameLayerDistXYZ[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; sameLayerDistEtaPhi[0] = m_cfg.sameLayerDistEtaPhi[0]; sameLayerDistEtaPhi[1] = m_cfg.sameLayerDistEtaPhi[1] / dd4hep::rad; diffLayerDistEtaPhi[0] = m_cfg.diffLayerDistEtaPhi[0]; @@ -89,6 +95,16 @@ void ImagingTopoCluster::init() { "Local [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].", sameLayerDistXY[0], sameLayerDistXY[1]); break; + case ImagingTopoClusterConfig::ELayerMode::xyz: + if (m_cfg.sameLayerDistXYZ.size() != 3) { + const std::string msg = "Expected 3 values (x_dist, y_dist, z_dist) for sameLayerDistXYZ"; + error(msg); + throw std::runtime_error(msg); + } + info("Same-layer clustering (same sector and same layer): " + "Local [x, y, z] distance between hits <= [{:.4f} mm, {:.4f} mm, {:.4f} mm].", + sameLayerDistXYZ[0], sameLayerDistXYZ[1], sameLayerDistXYZ[2]); + break; case ImagingTopoClusterConfig::ELayerMode::etaphi: if (m_cfg.sameLayerDistEtaPhi.size() != 2) { const std::string msg = "Expected 2 values (eta_dist, phi_dist) for sameLayerDistEtaPhi"; @@ -135,6 +151,16 @@ 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, y_dist) for diffLayerDistXYZ"; + error(msg); + throw std::runtime_error(msg); + } + info("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, diffLayerDistXYZ[0], diffLayerDistXYZ[1], diffLayerDistXYZ[2]); + break; case ImagingTopoClusterConfig::ELayerMode::tz: if (m_cfg.diffLayerDistTZ.size() != 2) { const std::string msg = "Expected 2 values (t_dist, z_dist) for diffLayerDistTZ"; @@ -260,6 +286,11 @@ bool ImagingTopoCluster::is_neighbour(const edm4eic::CalorimeterHit& h1, 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]) && @@ -296,6 +327,11 @@ bool ImagingTopoCluster::is_neighbour(const edm4eic::CalorimeterHit& h1, 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())); diff --git a/src/algorithms/calorimetry/ImagingTopoCluster.h b/src/algorithms/calorimetry/ImagingTopoCluster.h index 1d7d30aae9..1d7654b09b 100644 --- a/src/algorithms/calorimetry/ImagingTopoCluster.h +++ b/src/algorithms/calorimetry/ImagingTopoCluster.h @@ -55,6 +55,8 @@ class ImagingTopoCluster : public ImagingTopoClusterAlgorithm, // unitless counterparts of the input parameters std::array sameLayerDistXY{0, 0}; std::array diffLayerDistXY{0, 0}; + std::array sameLayerDistXYZ{0, 0, 0}; + std::array diffLayerDistXYZ{0, 0, 0}; std::array sameLayerDistEtaPhi{0, 0}; std::array diffLayerDistEtaPhi{0, 0}; std::array sameLayerDistTZ{0, 0}; diff --git a/src/algorithms/calorimetry/ImagingTopoClusterConfig.h b/src/algorithms/calorimetry/ImagingTopoClusterConfig.h index 01bc1f1845..639a5d0dc0 100644 --- a/src/algorithms/calorimetry/ImagingTopoClusterConfig.h +++ b/src/algorithms/calorimetry/ImagingTopoClusterConfig.h @@ -15,22 +15,26 @@ struct ImagingTopoClusterConfig { // 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 layerMode==xy) + // 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 global (eta, phi) to be considered as neighbors at same layers (if layerMode==etaphi) + // 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 layerMode==tz) + // 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 layerMode==xy) + // maximum distance of global (x, y) to be considered as neighbors at different layers (if samelayerMode==xy) std::vector> diffLayerDistXY = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; - // maximum distance of global (eta, phi) to be considered as neighbors at different layers (if layerMode==etaphi) + // 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 samelayerMode==etaphi) std::vector diffLayerDistEtaPhi = {0.01, 0.01}; - // maximum distance of global (t, z) to be considered as neighbors at different layers (if layerMode==tz) + // maximum distance of global (t, z) to be considered as neighbors at different layers (if samelayerMode==tz) std::vector diffLayerDistTZ = {1.0 * dd4hep::mm, 1.0 * dd4hep::mm}; // Layermodes - enum class ELayerMode { etaphi = 0, xy = 1, tz = 2 }; + enum class ELayerMode { etaphi = 0, xy = 1, xyz = 2, tz = 3 }; // 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 // determines how neighbors are determined for hits in different layers (using either eta and phi, or x and y) @@ -57,7 +61,9 @@ 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") { + } else if (s == "xyz" or s == "2") { + layerMode = ImagingTopoClusterConfig::ELayerMode::xyz; + } else if (s == "tz" or s == "3") { layerMode = ImagingTopoClusterConfig::ELayerMode::tz; } else { in.setstate(std::ios::failbit); // Set the fail bit if the input is not valid @@ -73,6 +79,9 @@ std::ostream& operator<<(std::ostream& out, const ImagingTopoClusterConfig::ELay case ImagingTopoClusterConfig::ELayerMode::xy: out << "xy"; break; + case ImagingTopoClusterConfig::ELayerMode::xyz: + out << "xyz"; + break; case ImagingTopoClusterConfig::ELayerMode::tz: out << "tz"; break; diff --git a/src/detectors/BEMC/BEMC.cc b/src/detectors/BEMC/BEMC.cc index 6e26e58adf..ed006fbbfa 100644 --- a/src/detectors/BEMC/BEMC.cc +++ b/src/detectors/BEMC/BEMC.cc @@ -203,7 +203,7 @@ void InitPlugin(JApplication* app) { .readout = "EcalBarrelScFiHits", .layerField = "layer", .sectorField = "sector", - .localDetFields = {"system"}, + .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", @@ -248,6 +248,43 @@ void InitPlugin(JApplication* app) { {"EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations"}, {.longitudinalShowerInfoAvailable = true, .energyWeight = "log", .logWeightBase = 6.2}, app)); + // Imaging TopoClustering on ScFi + app->Add(new JOmniFactoryGeneratorT( + "EcalBarrelScFiProtoClusters_Topo", {"EcalBarrelScFiRecHits"}, + {"EcalBarrelScFiProtoClusters_Topo"}, + { + .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}, + .sameLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xyz, + .diffLayerMode = eicrecon::ImagingTopoClusterConfig::ELayerMode::xyz, + .sectorDist = 5.0 * dd4hep::cm, + .minClusterHitEdep = 0, + .minClusterCenterEdep = 0, + .minClusterEdep = 100 * dd4hep::MeV, + .minClusterNhits = 10, + + }, + app // TODO: Remove me once fixed + )); + + app->Add(new JOmniFactoryGeneratorT( + "EcalBarrelScFiTopoClustersWithoutShapes", + {"EcalBarrelScFiProtoClusters_Topo", // edm4eic::ProtoClusterCollection + "EcalBarrelScFiRawHitAssociations"}, // edm4eic::MCRecoCalorimeterHitAssociation + {"EcalBarrelScFiTopoClustersWithoutShapes", // edm4eic::Cluster + "EcalBarrelScFiTopoClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation + {.energyWeight = "log", .sampFrac = 1.0, .logWeightBase = 6.2, .enableEtaBounds = false}, + app // TODO: Remove me once fixed + )); + app->Add(new JOmniFactoryGeneratorT( + "EcalBarrelScFiTopoClusters", + {"EcalBarrelScFiTopoClustersWithoutShapes", + "EcalBarrelScFiTopoClusterAssociationsWithoutShapes"}, + {"EcalBarrelScFiTopoClusters", "EcalBarrelScFiTopoClusterAssociations"}, + {.longitudinalShowerInfoAvailable = true, .energyWeight = "log", .logWeightBase = 6.2}, app)); + // Make sure digi and reco use the same value decltype(CalorimeterHitDigiConfig::capADC) EcalBarrelImaging_capADC = 8192; //8192, 13bit ADC decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalBarrelImaging_dyRangeADC = 3 * dd4hep::MeV; diff --git a/src/factories/calorimetry/ImagingTopoCluster_factory.h b/src/factories/calorimetry/ImagingTopoCluster_factory.h index c7deb78279..a6e9bb3a4a 100644 --- a/src/factories/calorimetry/ImagingTopoCluster_factory.h +++ b/src/factories/calorimetry/ImagingTopoCluster_factory.h @@ -31,6 +31,10 @@ class ImagingTopoCluster_factory config().diffLayerDistEtaPhi}; ParameterRef> m_ldtz_same{this, "sameLayerDistTZ", config().sameLayerDistTZ}; ParameterRef> m_ldtz_diff{this, "diffLayerDistTZ", config().diffLayerDistTZ}; + ParameterRef> m_ldxyz_same{this, "sameLayerDistXYZ", + config().sameLayerDistXYZ}; + ParameterRef> m_ldxyz_diff{this, "diffLayerDistXYZ", + config().diffLayerDistXYZ}; ParameterRef m_sameLayerMode{ this, "sameLayerMode", config().sameLayerMode}; ParameterRef m_diffLayerMode{ diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 3e36aeb0f9..84388bac49 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -307,6 +307,8 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "EcalBarrelScFiRecHits", "EcalBarrelScFiClusters", "EcalBarrelScFiClusterAssociations", + "EcalBarrelScFiTopoClusters", // ScFi clusters based in Topological Clustering + "EcalBarrelScFiTopoClusterAssociations", "EcalLumiSpecRawHits", "EcalLumiSpecRecHits", "EcalLumiSpecTruthClusters",