feat: Calorimeter hit digi simplification

src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h

@@ -56,9 +56,6 @@ namespace eicrecon {
     const dd4hep::Detector* m_detector;
     std::shared_ptr<spdlog::logger> m_log;
 
-    std::string m_input_simhit_tag;
-    std::string m_input_protoclust_tag;
-
     std::function<double(double, double, double, int)> weightFunc;
 
   private:

src/algorithms/calorimetry/CalorimeterHitDigi.cc

@@ -18,10 +18,8 @@
 #include <fmt/format.h>
 using namespace dd4hep;
 
-//
-// This algorithm converted from:
-//
-//  https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugDigi/src/components/CalorimeterHitDigi.cpp
+namespace eicrecon {
+
 //
 // TODO:
 // - Array type configuration parameters are not yet supported in JANA (needs to be added)
@@ -31,13 +29,9 @@ using namespace dd4hep;
 //   values here. This needs to be confirmed.
 
 
-
-//------------------------
-// AlgorithmInit
-//------------------------
-void CalorimeterHitDigi::AlgorithmInit(std::shared_ptr<spdlog::logger>& logger) {
-
-    // Assume all configuration parameter data members have been filled in already.
+void CalorimeterHitDigi::init(const dd4hep::Detector* detector, std::shared_ptr<spdlog::logger>& logger) {
+    m_detector = detector;
+    m_log = logger;
 
     // Gaudi implements a random number generator service. It is not clear to me how this
     // can work. There are multiple race conditions that occur in parallel event processing:
@@ -54,68 +48,53 @@ void CalorimeterHitDigi::AlgorithmInit(std::shared_ptr<spdlog::logger>& logger)
     // now, just use default values defined in header file.
 
     // set energy resolution numbers
-    m_log=logger;
-
-    if (u_eRes.empty()) {
-      u_eRes.resize(3);
-    } else if (u_eRes.size() != 3) {
-      m_log->error("Invalid u_eRes.size()");
-      throw std::runtime_error("Invalid u_eRes.size()");
+    if (m_cfg.eRes.empty()) {
+      m_cfg.eRes.resize(3);
+    } else if (m_cfg.eRes.size() != 3) {
+      m_log->error("Invalid m_cfg.eRes.size()");
+      throw std::runtime_error("Invalid m_cfg.eRes.size()");
     }
 
     // using juggler internal units (GeV, mm, radian, ns)
-    tRes       = m_tRes / dd4hep::ns;
-    stepTDC    = dd4hep::ns / m_resolutionTDC;
+    tRes       = m_cfg.tRes / dd4hep::ns;
+    stepTDC    = dd4hep::ns / m_cfg.resolutionTDC;
 
     // all these are for signal sum at digitization level
     merge_hits = false;
-    if (!u_fields.empty()) {
+    if (!m_cfg.fields.empty()) {
         // sanity checks
-        if (!m_geoSvc) {
-            m_log->error("Unable to locate Geometry Service.");
+        if (!m_detector) {
+            m_log->error("Unable to locate geometry.");
             throw std::runtime_error("Unable to locate Geometry Service.");
         }
-        if (m_readout.empty()) {
+        if (m_cfg.readout.empty()) {
             m_log->error("readoutClass is not provided, it is needed to know the fields in readout ids.");
             throw std::runtime_error("readoutClass is not provided.");
         }
 
         // get decoders
         try {
-            auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec();
+            auto id_desc = m_detector->readout(m_cfg.readout).idSpec();
             id_mask = 0;
-            for (size_t i = 0; i < u_fields.size(); ++i) {
-                id_mask |= id_desc.field(u_fields[i])->mask();
+            for (size_t i = 0; i < m_cfg.fields.size(); ++i) {
+                id_mask |= id_desc.field(m_cfg.fields[i])->mask();
             }
         } catch (...) {
             // a workaround to avoid breaking the whole analysis if a field is not in some configurations
             // TODO: it should be a fatal error to not cause unexpected analysis results
-            m_log->warn("Failed to load ID decoder for {}, hits will not be merged.", m_readout);
-            // throw::runtime_error(fmt::format("Failed to load ID decoder for {}", m_readout));
+            m_log->warn("Failed to load ID decoder for {}, hits will not be merged.", m_cfg.readout);
+            // throw::runtime_error(fmt::format("Failed to load ID decoder for {}", m_cfg.readout));
             return;
         }
         id_mask = ~id_mask;
-        m_log->info("ID mask in {:s}: {:#064b}", m_readout, id_mask);
+        m_log->info("ID mask in {:s}: {:#064b}", m_cfg.readout, id_mask);
         // all checks passed
         merge_hits = true;
     }
 }
 
 
-
-//------------------------
-// AlgorithmChangeRun
-//------------------------
-void CalorimeterHitDigi::AlgorithmChangeRun() {
-    /// This is automatically run before Process, when a new run number is seen
-    /// Usually we update our calibration constants by asking a JService
-    /// to give us the latest data for this run number
-}
-
-//------------------------
-// AlgorithmProcess
-//------------------------
-std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::AlgorithmProcess(const edm4hep::SimCalorimeterHitCollection &simhits)  {
+std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::process(const edm4hep::SimCalorimeterHitCollection &simhits)  {
     if (merge_hits) {
         return std::move(signal_sum_digi(simhits));
     } else {
@@ -132,35 +111,31 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::single
         const double eDep    = ahit.getEnergy();
 
         // apply additional calorimeter noise to corrected energy deposit
-        const double eResRel = (eDep > m_threshold)
+        const double eResRel = (eDep > m_cfg.threshold)
                                ? m_normDist(generator) * std::sqrt(
-                                    std::pow(u_eRes[0] / std::sqrt(eDep), 2) +
-                                    std::pow(u_eRes[1], 2) +
-                                    std::pow(u_eRes[2] / (eDep), 2)
-                )
+                                    std::pow(m_cfg.eRes[0] / std::sqrt(eDep), 2) +
+                                    std::pow(m_cfg.eRes[1], 2) +
+                                    std::pow(m_cfg.eRes[2] / (eDep), 2)
+                                 )
                                : 0;
-//       const double eResRel = (eDep > 1e-6)
-//                               ? m_normDist(generator) * std::sqrt(std::pow(u_eRes[0] / std::sqrt(eDep), 2) +
-//                                                          std::pow(u_eRes[1], 2) + std::pow(u_eRes[2] / (eDep), 2))
-//                               : 0;
 
-        const double ped    = m_pedMeanADC + m_normDist(generator) * m_pedSigmaADC;
-        const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / m_dyRangeADC * m_capADC);
+        const double ped    = m_cfg.pedMeanADC + m_normDist(generator) * m_cfg.pedSigmaADC;
+        const long long adc = std::llround(ped + eDep * (m_cfg.corrMeanScale + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC);
 
         double time = std::numeric_limits<double>::max();
         for (const auto& c : ahit.getContributions()) {
             if (c.getTime() <= time) {
                 time = c.getTime();
             }
         }
-        if (time > m_capTime) continue;
+        if (time > m_cfg.capTime) continue;
 
         const long long tdc = std::llround((time + m_normDist(generator) * tRes) * stepTDC);
 
-        if (eDep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t cell ID {}", eDep, adc, time, m_capTime, tdc, ahit.getCellID());
+        if (eDep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t cell ID {}", eDep, adc, time, m_cfg.capTime, tdc, ahit.getCellID());
         rawhits->create(
                 ahit.getCellID(),
-                (adc > m_capADC ? m_capADC : adc),
+                (adc > m_cfg.capADC ? m_cfg.capADC : adc),
                 tdc
         );
     }
@@ -177,8 +152,8 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::signal
     for (const auto &ahit : simhits) {
         uint64_t hid = ahit.getCellID() & id_mask;
 
-        m_log->trace("org cell ID in {:s}: {:#064b}", m_readout, ahit.getCellID());
-        m_log->trace("new cell ID in {:s}: {:#064b}", m_readout, hid);
+        m_log->trace("org cell ID in {:s}: {:#064b}", m_cfg.readout, ahit.getCellID());
+        m_log->trace("new cell ID in {:s}: {:#064b}", m_cfg.readout, hid);
 
         merge_map[hid].push_back(ix);
 
@@ -202,7 +177,7 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::signal
                     timeC = c.getTime();
                 }
             }
-            if (timeC > m_capTime) continue;
+            if (timeC > m_cfg.capTime) continue;
             edep += hit.getEnergy();
             m_log->trace("adding {} \t total: {}", hit.getEnergy(), edep);
 
@@ -221,29 +196,25 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::signal
             }
         }
 
-//        double eResRel = 0.;
         // safety check
-        const double eResRel = (edep > m_threshold)
-                ? m_normDist(generator) * u_eRes[0] / std::sqrt(edep) +
-                  m_normDist(generator) * u_eRes[1] +
-                  m_normDist(generator) * u_eRes[2] / edep
+        const double eResRel = (edep > m_cfg.threshold)
+                ? m_normDist(generator) * m_cfg.eRes[0] / std::sqrt(edep) +
+                  m_normDist(generator) * m_cfg.eRes[1] +
+                  m_normDist(generator) * m_cfg.eRes[2] / edep
                   : 0;
-//        if (edep > 1e-6) {
-//            eResRel = m_normDist(generator) * u_eRes[0] / std::sqrt(edep) +
-//                      m_normDist(generator) * u_eRes[1] +
-//                      m_normDist(generator) * u_eRes[2] / edep;
-//        }
-        double    ped     = m_pedMeanADC + m_normDist(generator) * m_pedSigmaADC;
-        unsigned long long adc     = std::llround(ped + edep * (m_corrMeanScale + eResRel) / m_dyRangeADC * m_capADC);
+        double    ped     = m_cfg.pedMeanADC + m_normDist(generator) * m_cfg.pedSigmaADC;
+        unsigned long long adc     = std::llround(ped + edep * (m_cfg.corrMeanScale + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC);
         unsigned long long tdc     = std::llround((time + m_normDist(generator) * tRes) * stepTDC);
 
-        if (edep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {}", edep, adc, time, m_capTime, tdc);
+        if (edep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {}", edep, adc, time, m_cfg.capTime, tdc);
         rawhits->create(
                 mid,
-                (adc > m_capADC ? m_capADC : adc),
+                (adc > m_cfg.capADC ? m_cfg.capADC : adc),
                 tdc
         );
     }
 
     return std::move(rawhits);
 }
+
+} // namespace eicrecon

src/algorithms/calorimetry/CalorimeterHitDigi.h

@@ -16,72 +16,42 @@
 #include <memory>
 #include <random>
 
-#include "services/geometry/dd4hep/JDD4hep_service.h"
+#include <DD4hep/Detector.h>
 
 #include <edm4hep/SimCalorimeterHitCollection.h>
 #include <edm4hep/RawCalorimeterHitCollection.h>
 #include <spdlog/spdlog.h>
 
-class CalorimeterHitDigi {
+#include "algorithms/interfaces/WithPodConfig.h"
+#include "CalorimeterHitDigiConfig.h"
 
-    // Insert any member variables here
+namespace eicrecon {
 
-public:
-    CalorimeterHitDigi() = default;
-    void AlgorithmInit(std::shared_ptr<spdlog::logger>& logger);
-    void AlgorithmChangeRun() ;
-    std::unique_ptr<edm4hep::RawCalorimeterHitCollection> AlgorithmProcess(const edm4hep::SimCalorimeterHitCollection &simhits) ;
+  class CalorimeterHitDigi : public WithPodConfig<CalorimeterHitDigiConfig> {
 
-    //-------- Configuration Parameters ------------
-    //instantiate new spdlog logger
-    std::shared_ptr<spdlog::logger> m_log;
+  public:
+    void init(const dd4hep::Detector* detector, std::shared_ptr<spdlog::logger>& logger);
+    std::unique_ptr<edm4hep::RawCalorimeterHitCollection> process(const edm4hep::SimCalorimeterHitCollection &simhits) ;
 
-    // Name of input data type (collection)
-    std::string              m_input_tag;
-
-    // additional smearing resolutions
-    std::vector<double>      u_eRes;
-    double                   m_tRes;
-
-    // single hit energy deposition threshold
-    double                   m_threshold=1.0*dd4hep::keV; // {this, "threshold", 1. * keV};
-
-    // digitization settings
-    unsigned int             m_capADC;
-    double                   m_capTime = 1000.;
-    double                   m_dyRangeADC;
-    unsigned int             m_pedMeanADC;
-    double                   m_pedSigmaADC;
-    double                   m_resolutionTDC;
-    double                   m_corrMeanScale;
-
-    // signal sums
-    std::vector<std::string> u_fields;
-    std::string              m_geoSvcName;
-    std::string              m_readout;
-
-    // This may be used to declare the data members as JANA configuration parameters.
-    // This should compile OK even without JANA so long as you don't try using it.
-    // To use it, do something like the following:
-    //
-    //    mycalohitdigi->SetJANAConfigParameters( japp, "BEMC");
-    //
-    // The above will register config. parameters like: "BEMC:tag".
-    // The configuration parameter members of this class should be set to their
-    // defaults *before* calling this.
-    //-----------------------------------------------
+  private:
 
     // unitless counterparts of inputs
     double           dyRangeADC{0}, stepTDC{0}, tRes{0};
     // variables for merging at digitization step
     bool             merge_hits = false;
-    std::shared_ptr<JDD4hep_service> m_geoSvc;
+
     uint64_t         id_mask{0};
 
-private:
+  private:
+    const dd4hep::Detector* m_detector;
+    std::shared_ptr<spdlog::logger> m_log;
+
     std::default_random_engine generator; // TODO: need something more appropriate here
     std::normal_distribution<double> m_normDist; // defaults to mean=0, sigma=1
 
     std::unique_ptr<edm4hep::RawCalorimeterHitCollection> single_hits_digi(const edm4hep::SimCalorimeterHitCollection &simhits);
     std::unique_ptr<edm4hep::RawCalorimeterHitCollection> signal_sum_digi(const edm4hep::SimCalorimeterHitCollection &simhits);
-};
+
+  };
+
+} // namespace eicrecon

src/algorithms/calorimetry/CalorimeterHitDigiConfig.h

@@ -0,0 +1,34 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2023 Wouter Deconinck
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+namespace eicrecon {
+
+  struct CalorimeterHitDigiConfig {
+
+    std::vector<double>      eRes;
+    double                   tRes;
+
+    // single hit energy deposition threshold
+    double                   threshold{1.0*dd4hep::keV};
+
+    // digitization settings
+    unsigned int             capADC{1};
+    double                   capTime{1000}; // dynamic range in ns
+    double                   dyRangeADC{1};
+    unsigned int             pedMeanADC{0};
+    double                   pedSigmaADC{0};
+    double                   resolutionTDC{1};
+    double                   corrMeanScale{1};
+
+    // signal sums
+    std::string              readout{""};
+    std::vector<std::string> fields{};
+
+  };
+
+} // eicrecon