Scfi Img CrossSystem Topo Clustering (#2861)

src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc

@@ -8,6 +8,9 @@
  *  Author: Chao Peng (ANL), 09/27/2020
  */
 
+#include <DD4hep/IDDescriptor.h>
+#include <DD4hep/Readout.h>
+#include <DDSegmentation/BitFieldCoder.h>
 #include <Evaluator/DD4hepUnits.h>
 #include <boost/algorithm/string/join.hpp>
 #include <boost/range/adaptor/map.hpp>
@@ -27,17 +30,31 @@
 #include <limits>
 #include <map>
 #include <optional>
+#include <stdexcept>
 #include <tuple>
 #include <vector>
 
-#include "CalorimeterClusterRecoCoG.h"
+#include "algorithms/calorimetry/CalorimeterClusterRecoCoG.h"
 #include "algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h"
 
 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<edm4eic::MutableCluster>
 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<float>::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=" <<totalE << " log(E/totalE)=" << std::log(hit.getEnergy()/totalE) << std::endl;
-    float w = weightFunc(hit.getEnergy() * weight, totalE, logWeightBase, 0);
-    tw += w;
-    v = v + (hit.getPosition() * w);
+    const auto& hit = pcl.getHits()[i];
+    float w         = 0.0f;
+
+    if (specialMode) {
+
+      if (isImagingHit(hit, m_idSpec)) {
+        const auto weight = pcl.getWeights()[i];
+        w                 = weightFunc(hit.getEnergy() * weight, totalE, logWeightBase, 0);
+      } else if (isSciFiHit(hit, m_idSpec)) {
+        w = 0.0f; // ScFi has no position weight
+      }
+
+      tw += w;
+      v = v + hit.getPosition() * w;
+    } else {
+      const auto weight = pcl.getWeights()[i];
+      float w           = weightFunc(hit.getEnergy() * weight, totalE, logWeightBase, 0);
+      tw += w;
+      v = v + (hit.getPosition() * w);
+    }
   }
   if (tw == 0.) {
     warning("zero total weights encountered, you may want to adjust your weighting parameter.");

src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h

@@ -10,7 +10,10 @@
 
 #pragma once
 
+#include <DD4hep/Detector.h>
+#include <DD4hep/IDDescriptor.h>
 #include <algorithms/algorithm.h>
+#include <algorithms/geo.h>
 #include <edm4eic/ClusterCollection.h>
 #include <edm4eic/MCRecoCalorimeterHitAssociationCollection.h>
 #include <edm4eic/MCRecoCalorimeterHitLinkCollection.h>
@@ -23,6 +26,7 @@
 #include <algorithm>
 #include <cmath>
 #include <functional>
+#include <gsl/pointers>
 #include <map>
 #include <memory>
 #include <optional>
@@ -81,6 +85,10 @@ class CalorimeterClusterRecoCoG : public CalorimeterClusterRecoCoGAlgorithm,
 private:
   std::function<double(double, double, double, int)> weightFunc;
 
+  // Pointer to the geometry service
+  dd4hep::IDDescriptor m_idSpec;
+  const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
+
 private:
   std::optional<edm4eic::MutableCluster> reconstruct(const edm4eic::ProtoCluster& pcl) const;
   void associate(const edm4eic::Cluster& cl,

src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h

@@ -10,7 +10,7 @@
 namespace eicrecon {
 
 struct CalorimeterClusterRecoCoGConfig {
-
+  std::string readout = "";
   std::string energyWeight;
 
   double sampFrac      = 1.;