Imaging topological clustering for Scfi

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()));

src/algorithms/calorimetry/ImagingTopoCluster.h

@@ -55,6 +55,8 @@ class ImagingTopoCluster : public ImagingTopoClusterAlgorithm,
   // unitless counterparts of the input parameters
   std::array<double, 2> sameLayerDistXY{0, 0};
   std::array<double, 2> diffLayerDistXY{0, 0};
+  std::array<double, 3> sameLayerDistXYZ{0, 0, 0};
+  std::array<double, 3> diffLayerDistXYZ{0, 0, 0};
   std::array<double, 2> sameLayerDistEtaPhi{0, 0};
   std::array<double, 2> diffLayerDistEtaPhi{0, 0};
   std::array<double, 2> sameLayerDistTZ{0, 0};

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<std::variant<std::string, double>> 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<double> 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<double> 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<double> 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<std::variant<std::string, double>> 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<double> 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<double> 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<double> 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;

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<ImagingTopoCluster_factory>(
+      "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<CalorimeterClusterRecoCoG_factory>(
+      "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<CalorimeterClusterShape_factory>(
+      "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;

src/factories/calorimetry/ImagingTopoCluster_factory.h

@@ -31,6 +31,10 @@ class ImagingTopoCluster_factory
                                                 config().diffLayerDistEtaPhi};
   ParameterRef<std::vector<double>> m_ldtz_same{this, "sameLayerDistTZ", config().sameLayerDistTZ};
   ParameterRef<std::vector<double>> m_ldtz_diff{this, "diffLayerDistTZ", config().diffLayerDistTZ};
+  ParameterRef<std::vector<double>> m_ldxyz_same{this, "sameLayerDistXYZ",
+                                                 config().sameLayerDistXYZ};
+  ParameterRef<std::vector<double>> m_ldxyz_diff{this, "diffLayerDistXYZ",
+                                                 config().diffLayerDistXYZ};
   ParameterRef<eicrecon::ImagingTopoClusterConfig::ELayerMode> m_sameLayerMode{
       this, "sameLayerMode", config().sameLayerMode};
   ParameterRef<eicrecon::ImagingTopoClusterConfig::ELayerMode> m_diffLayerMode{

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",