eic · mhkim-anl · Apr 1, 2025 · Apr 1, 2025 · Apr 1, 2025 · Apr 1, 2025
diff --git a/src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc b/src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc
@@ -0,0 +1,224 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Minho Kim, Sylvester Joosten, Derek Anderson, Wouter Deconinck
+
+#include "SimCalorimeterHitProcessor.h"
+
+#include <DD4hep/Detector.h>
+#include <DD4hep/IDDescriptor.h>
+#include <DD4hep/Readout.h>
+#include <DD4hep/config.h>
+#include <DDSegmentation/BitFieldCoder.h>
+#include <Evaluator/DD4hepUnits.h>
+#include <edm4hep/Vector3f.h>
+#include <fmt/core.h>
+#include <podio/ObjectID.h>
+#include <podio/RelationRange.h>
+#include <edm4eic/unit_system.h>
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <gsl/pointers>
+#include <limits>
+#include <map>
+#include <numeric>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <vector>
+#include <ranges>
+
+#include "algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h"
+
+using namespace dd4hep;
+
+namespace eicrecon {
+
+void SimCalorimeterHitProcessor::init() {
+
+  // readout checks
+  if (m_cfg.readout.empty()) {
+    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 {
+    m_id_spec = m_geo.detector()->readout(m_cfg.readout).idSpec();
+  } catch (...) {
+    debug("Failed to load ID decoder for {}", m_cfg.readout);
+    throw std::runtime_error(fmt::format("Failed to load ID decoder for {}", m_cfg.readout));
+  }
+
+  // get m_hit_id_mask for adding up hits with the same dimensions that are merged over
+  if (!m_cfg.hitMergeFields.empty()) {
+    uint64_t id_inverse_mask = 0;
+    for (auto& field : m_cfg.hitMergeFields) {
+      id_inverse_mask |= m_id_spec.field(field)->mask();
+    }
+    m_hit_id_mask = ~id_inverse_mask;
+    if (m_hit_id_mask)
+      debug("ID mask in {:s}: {:#064b}", m_cfg.readout, m_hit_id_mask.value());
+  }
+
+  // get m_contribution_id_mask for adding up contributions with the same dimensions that are merged over
+  if (!m_cfg.contributionMergeFields.empty()) {
+    uint64_t id_inverse_mask = 0;
+    for (auto& field : m_cfg.contributionMergeFields) {
+      id_inverse_mask |= m_id_spec.field(field)->mask();
+    }
+    m_contribution_id_mask = ~id_inverse_mask;
+    if (m_contribution_id_mask)
+      debug("ID mask in {:s}: {:#064b}", m_cfg.readout, m_contribution_id_mask.value());
+  }
+
+  // get reference position for attenuating hits and contributions
+  if (!m_cfg.attenuationReferencePositionName.empty()) {
+    m_attenuationReferencePosition_mm =
+        m_geo.detector()->constant<double>(m_cfg.attenuationReferencePositionName) *
+        edm4eic::unit::mm / dd4hep::mm;
+  }
+}
+
+void SimCalorimeterHitProcessor::process(const SimCalorimeterHitProcessor::Input& input,
+                                         const SimCalorimeterHitProcessor::Output& output) const {
+
+  const auto [in_hits]              = input;
+  auto [out_hits, out_hit_contribs] = output;
+
+  // map for regrouping
+  std::map<edm4hep::MCParticle, std::vector<edm4hep::MutableSimCalorimeterHit>> mapMCParToSimCalHit;
+
+  // regroup the sim hits by mc particle
+  // keep an eye on whether these for-loops grow large enough to affect memory limits
+  // on the running of reconstruction, especially forward and backward calorimeters
+  for (const auto& ih : *in_hits) {
+    for (const auto& contrib : ih.getContributions()) {
+      edm4hep::MCParticle primary = get_primary(contrib);
+
+      auto& simhit = mapMCParToSimCalHit[primary].emplace_back(ih.getCellID(), contrib.getEnergy(),
+                                                               ih.getPosition());
+      simhit.addToContributions(contrib);
+
+      trace("Identified primary: id = {}, pid = {}, total energy = {}, contributed = {}",
+            primary.getObjectID().index, primary.getPDG(), primary.getEnergy(),
+            mapMCParToSimCalHit[primary].back().getEnergy());
+    }
+  }
+
+  // Attenuate energies of the sim hits
+  // 1. sum the hits if they have the same dimension
+  // 2. attenuate the summed hits
+  for (const auto& [par, hits] : mapMCParToSimCalHit) {<
8000
/span>
+    double attFactor = 1.;
+
+    // when merging hits is necessary
+    if (m_hit_id_mask) {
+      // map
+      std::unordered_map<uint64_t, std::vector<std::size_t>> hit_merge_map;
+
+      // map for merging hits
+      std::size_t ix = 0;
+      for (const auto& ahit : hits) {
+        uint64_t hid = ahit.getCellID() & m_hit_id_mask.value();
+        hit_merge_map[hid].push_back(ix);
+        ix++;
+      }
+
+      for (const auto& [id_hit, ixs_hit] : hit_merge_map) {
+        std::unordered_map<uint64_t, std::vector<std::size_t>> contribution_merge_map;
+
+        // map for merging contributions
+        for (const auto& ix : ixs_hit) {
+          uint64_t hid = hits[ix].getCellID() & m_contribution_id_mask.value();
+          contribution_merge_map[hid].push_back(ix);
+        }
+
+        auto out_hit   = out_hits->create();
+        double edepSum = 0;
+
+        for (const auto& [id_contrib, ixs_contrib] : contribution_merge_map) {
+          auto leading_hit     = hits[ixs_contrib[0]];
+          auto leading_contrib = hits[ixs_contrib[0]].getContributions(0);
+
+          // add up the energy deposits considering the attenuation
+          double edepSumLocal = std::accumulate(
+              ixs_contrib.begin(), ixs_contrib.end(), 0.0, [&hits](double sum, size_t ix) {
+                return sum + hits[ix].getContributions(0).getEnergy();
+              });
+
+          // attenuation
+          if (m_attenuationReferencePosition_mm) {
+            attFactor = get_attenuation(leading_hit.getPosition().z);
+            edepSumLocal *= attFactor;
+          }
+
+          edepSum += edepSumLocal;
+
+          // find the earliest time
+          float timeEarLocal =
+              std::ranges::min(ixs_contrib | std::views::transform([&hits](size_t ix) {
+                                 return hits[ix].getContributions(0).getTime();
+                               }));
+
+          auto out_hit_contrib = out_hit_contribs->create();
+          out_hit_contrib.setPDG(leading_contrib.getPDG());
+          out_hit_contrib.setEnergy(static_cast<float>(edepSumLocal));
+          out_hit_contrib.setTime(timeEarLocal);
+          out_hit_contrib.setStepPosition(leading_contrib.getStepPosition());
+          out_hit_contrib.setParticle(par);
+          out_hit.addToContributions(out_hit_contrib);
+        }
+
+        auto leading_hit = hits[ixs_hit[0]];
+
+        out_hit.setCellID(leading_hit.getCellID());
+        out_hit.setEnergy(static_cast<float>(edepSum));
+        out_hit.setPosition(leading_hit.getPosition());
+      }
+      // when merging hits are not necessary
+    } else {
+      for (const auto& hit : hits) {
+        auto contrib = hit.getContributions(0);
+
+        // attenuation
+        if (m_attenuationReferencePosition_mm) {
+          attFactor = get_attenuation(hit.getPosition().z);
+        }
+
+        auto out_hit_contrib = out_hit_contribs->create();
+        out_hit_contrib.setPDG(contrib.getPDG());
+        out_hit_contrib.setEnergy(attFactor * contrib.getEnergy());
+        out_hit_contrib.setTime(contrib.getTime());
+        out_hit_contrib.setStepPosition(contrib.getStepPosition());
+        out_hit_contrib.setParticle(par);
+
+        auto out_hit = out_hits->create();
+        out_hit.setCellID(hit.getCellID());
+        out_hit.setEnergy(attFactor * hit.getEnergy());
+        out_hit.setPosition(hit.getPosition());
+        out_hit.addToContributions(out_hit_contrib);
+      }
+    }
+  }
+}
+
+edm4hep::MCParticle
+SimCalorimeterHitProcessor::get_primary(const edm4hep::CaloHitContribution& contrib) const {
+  const auto contributor = contrib.getParticle();
+
+  edm4hep::MCParticle primary = contributor;
+  while (primary.parents_size() > 0) {
+    if (primary.getGeneratorStatus() != 0)
+      break;
+    primary = primary.getParents(0);
+  }
+  return primary;
+}
+
+double SimCalorimeterHitProcessor::get_attenuation(double zpos) const {
+  double length = std::abs(m_attenuationReferencePosition_mm.value() - zpos);
+  double factor = m_cfg.attPars[0] * std::exp(-length / m_cfg.attPars[1]) +
+                  (1 - m_cfg.attPars[0]) * std::exp(-length / m_cfg.attPars[2]);
+  return factor;
+}
+} // namespace eicrecon
diff --git a/src/algorithms/calorimetry/SimCalorimeterHitProcessor.h b/src/algorithms/calorimetry/SimCalorimeterHitProcessor.h
@@ -0,0 +1,60 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Minho Kim, Sylvester Joosten, Derek Anderson, Wouter Deconinck
+
+#pragma once
+
+#include <DD4hep/IDDescriptor.h>
+#include <algorithms/algorithm.h>
+#include <algorithms/geo.h>
+#include <edm4hep/CaloHitContributionCollection.h>
+#include <edm4hep/MCParticleCollection.h>
+#include <edm4hep/SimCalorimeterHitCollection.h>
+#include <stdint.h>
+#include <optional>
+#include <string>
+#include <string_view>
+
+#include "SimCalorimeterHitProcessorConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using SimCalorimeterHitProcessorAlgorithm =
+    algorithms::Algorithm<algorithms::Input<edm4hep::SimCalorimeterHitCollection>,
+                          algorithms::Output<edm4hep::SimCalorimeterHitCollection,
+                                             edm4hep::CaloHitContributionCollection>>;
+
+class SimCalorimeterHitProcessor : public SimCalorimeterHitProcessorAlgorithm,
+                                   public WithPodConfig<SimCalorimeterHitProcessorConfig> {
+
+public:
+  SimCalorimeterHitProcessor(std::string_view name)
+      : SimCalorimeterHitProcessorAlgorithm{
+            name,
+            {"inputHitCollection"},
+            {"outputHitCollection", "outputHitContributionCollection"},
+            "Regroup the hits by particle, add up the hits if"
+            "they have the same z-segmentation, and attenuate."} {}
+
+  void init() final;
+  void process(const Input&, const Output&) const final;
+
+private:
+  std::optional<uint64_t> m_hit_id_mask;
+  std::optional<uint64_t> m_contribution_id_mask;
+
+  dd4hep::IDDescriptor m_id_spec;
+
+  
F438
const algorithms::GeoSvc& m_geo = algorithms::GeoSvc::instance();
+
+  // a reference value for attenuation
+  std::optional<double> m_attenuationReferencePosition_mm;
+
+private:
+  edm4hep::MCParticle get_primary(const edm4hep::CaloHitContribution& contrib) const;
+
+  // attenuation function
+  double get_attenuation(double zpos) const;
+};
+
+} // namespace eicrecon
diff --git a/src/algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h b/src/algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h
@@ -0,0 +1,24 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Minho Kim, Sylvester Joosten, Derek Anderson, Wouter Deconinck
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+namespace eicrecon {
+
+struct SimCalorimeterHitProcessorConfig {
+
+  // parameters for attenuation function
+  std::vector<double> attPars;
+
+  std::string readout{""};
+  std::string attenuationReferencePositionName{""};
+  // fields for merging hits
+  std::vector < std::string > hitMergeFields{};
+  // fildes for merging contributions
+  std::vector < std::string > contributionMergeFields{};
+};
+
+} // namespace eicrecon
diff --git a/src/detectors/BEMC/BEMC.cc b/src/detectors/BEMC/BEMC.cc
@@ -4,14 +4,17 @@
 //
 
 #include <edm4eic/EDM4eicVersion.h>
+#include <edm4eic/unit_system.h>
 #include <Evaluator/DD4hepUnits.h>
 #include <JANA/JApplication.h>
 #include <math.h>
 #include <string>
 
 #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h"
+#include "algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h"
 #include "extensions/jana/JOmniFactoryGeneratorT.h"
 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
+#include "factories/calorimetry/SimCalorimeterHitProcessor_factory.h"
 #include "factories/calorimetry/CalorimeterHitDigi_factory.h"
 #include "factories/calorimetry/CalorimeterHitReco_factory.h"
 #include "factories/calorimetry/CalorimeterIslandCluster_factory.h"
@@ -28,13 +31,43 @@ void InitPlugin(JApplication* app) {
 
   InitJANAPlugin(app);
 
+  // Make sure left and right use the same value
+      0.416212, 747.39875 * edm4eic::unit::mm, 7521.88383 * edm4eic::unit::mm};
+  decltype(SimCalorimeterHitProcessorConfig::hitMergeFields) EcalBarrelScFi_hitMergeFields = {"fiber", "z"};
+  decltype(SimCalorimeterHitProcessorConfig::contributionMergeFields) EcalBarrelScFi_contributionMergeFields = {"fiber"};
+
   // Make sure digi and reco use the same value
   decltype(CalorimeterHitDigiConfig::capADC) EcalBarrelScFi_capADC = 16384; //16384,  14bit ADC
   decltype(CalorimeterHitDigiConfig::dyRangeADC) EcalBarrelScFi_dyRangeADC   = 1500 * dd4hep::MeV;
   decltype(CalorimeterHitDigiConfig::pedMeanADC) EcalBarrelScFi_pedMeanADC   = 100;
   decltype(CalorimeterHitDigiConfig::pedSigmaADC) EcalBarrelScFi_pedSigmaADC = 1;
   decltype(CalorimeterHitDigiConfig::resolutionTDC) EcalBarrelScFi_resolutionTDC =
       10 * dd4hep::picosecond;
+  app->Add(new JOmniFactoryGeneratorT<SimCalorimeterHitProcessor_factory>(
+      "EcalBarrelScFiPAttenuatedHits", {"EcalBarrelScFiHits"},
+      {"EcalBarrelScFiPAttenuatedHits", "EcalBarrelScFiPAttenuatedHitContributions"},
+      {
+          .attPars                          = EcalBarrelScFi_attPars,
+          .readout                          = "EcalBarrelScFiHits",
+          .attenuationReferencePositionName = "EcalBarrel_Readout_zmax",
+          .hitMergeFields                   = EcalBarrelScFi_hitMergeFields,
+	  .contributionMergeFields          = EcalBarrelScFi_contributionMergeFields,
+      },
+      app // TODO: Remove me once fixed
+      ));
+  app->Add(new JOmniFactoryGeneratorT<SimCalorimeterHitProcessor_factory>(
+      "EcalBarrelScFiNAttenuatedHits", {"EcalBarrelScFiHits"},
+      {"EcalBarrelScFiNAttenuatedHits", "EcalBarrelScFiNAttenuatedHitContributions"},
+      {
+          .attPars                          = EcalBarrelScFi_attPars,
+          .readout                          = "EcalBarrelScFiHits",
+          .attenuationReferencePositionName = "EcalBarrel_Readout_zmin",
+	  .hitMergeFields                   = EcalBarrelScFi_hitMergeFields,
+          .contributionMergeFields          = EcalBarrelScFi_contributionMergeFields,
+      },
+      app // TODO: Remove me once fixed
+      ));
   app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
       "EcalBarrelScFiRawHits", {"EcalBarrelScFiHits"},
 #if EDM4EIC_VERSION_MAJOR >= 7