Pulse generation from time-clustered Simcalorimeterhits using a shared template

src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc

@@ -13,7 +13,6 @@
 #include <edm4hep/Vector3f.h>
 #include <edm4hep/utils/vector_utils.h>
 #include <fmt/core.h>
-#include <podio/ObjectID.h>
 #include <podio/RelationRange.h>
 #include <podio/podioVersion.h>
 #include <cmath>
@@ -57,12 +56,13 @@ template <> struct hash<edm4hep::MCParticle> {
 
 // Hash for tuple<edm4hep::MCParticle, uint64_t>
 // --> not yet supported by any compiler at the moment
-template <> struct hash<std::tuple<edm4hep::MCParticle, uint64_t>> {
-  size_t operator()(const std::tuple<edm4hep::MCParticle, uint64_t>& key) const noexcept {
-    const auto& [particle, cellID] = key;
-    size_t h1                      = hash<edm4hep::MCParticle>{}(particle);
-    size_t h2                      = hash<uint64_t>{}(cellID);
-    return h1 ^ (h2 << 1);
+template <> struct hash<std::tuple<edm4hep::MCParticle, uint64_t, int>> {
+  size_t operator()(const std::tuple<edm4hep::MCParticle, uint64_t, int>& key) const noexcept {
+    const auto& [particle, cellID, timeID] = key;
+    size_t h1                              = hash<edm4hep::MCParticle>{}(particle);
+    size_t h2                              = hash<uint64_t>{}(cellID);
+    size_t h3                              = hash<int>{}(timeID);
+    return ((h1 ^ (h2 << 1)) >> 1) ^ (h3 << 1);
   }
 };
 
@@ -173,7 +173,7 @@ void SimCalorimeterHitProcessor::process(const SimCalorimeterHitProcessor::Input
   // (hence modify) contributions which is not supported for PodIO VectorMembers. Using
   // reasonable contribution merging, at least the intermediary structure should be
   // quite a bit smaller than the original hit collection.
-  using HitIndex = std::tuple<edm4hep::MCParticle, uint64_t /* cellID */>;
+  using HitIndex = std::tuple<edm4hep::MCParticle, uint64_t /* cellID */, int /* timeID */>;
   std::unordered_map<HitIndex,
                      std::unordered_map<uint64_t /* cellID */, HitContributionAccumulator>>
       hit_map;
@@ -184,21 +184,22 @@ void SimCalorimeterHitProcessor::process(const SimCalorimeterHitProcessor::Input
         (ih.getCellID() & m_hit_id_mask.value() & m_contribution_id_mask.value());
     // the cell ID of this particular contribution (we are using contributions to store
     // the hits making up this "superhit" with more segmentation)
-    const uint64_t newcontrib_cellID = (ih.getCellID() & m_hit_id_mask.value());
+    const uint64_t newcontrib_cellID = (ih.getCellID() & m_contribution_id_mask.value());
     // Optional attenuation
     const double attFactor =
         m_attenuationReferencePosition ? get_attenuation(ih.getPosition().z) : 1.;
     // Use primary particle (traced back through parents) to group contributions
     for (const auto& contrib : ih.getContributions()) {
       edm4hep::MCParticle primary = lookup_primary(contrib);
-      auto& hit_accum             = hit_map[{primary, newhit_cellID}][newcontrib_cellID];
       const double propagationTime =
           m_attenuationReferencePosition
               ? std::abs(m_attenuationReferencePosition.value() - ih.getPosition().z) *
                     m_cfg.inversePropagationSpeed
               : 0.;
-      hit_accum.add(contrib.getEnergy() * attFactor,
-                    contrib.getTime() + propagationTime + m_cfg.fixedTimeDelay, ih.getPosition());
+      const double totalTime  = contrib.getTime() + propagationTime + m_cfg.fixedTimeDelay;
+      const int newhit_timeID = std::floor(totalTime / m_cfg.timeWindow);
+      auto& hit_accum         = hit_map[{primary, newhit_cellID, newhit_timeID}][newcontrib_cellID];
+      hit_accum.add(contrib.getEnergy() * attFactor, totalTime, ih.getPosition());
     }
   }
 
@@ -208,7 +209,7 @@ void SimCalorimeterHitProcessor::process(const SimCalorimeterHitProcessor::Input
 
     auto out_hit = out_hits->create();
 
-    const auto& [particle, cellID] = hit_idx;
+    const auto& [particle, cellID, timeID] = hit_idx;
     HitContributionAccumulator new_hit;
     for (const auto& [contrib_idx, contrib] : contribs) {
       // Aggregate contributions to for the global hit

src/algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h

@@ -5,6 +5,7 @@
 
 #include <string>
 #include <vector>
+#include <edm4eic/unit_system.h>
 
 namespace eicrecon {
 
@@ -27,6 +28,8 @@ struct SimCalorimeterHitProcessorConfig {
   double inversePropagationSpeed{};
   // detector-related time delay (e.g., scintillation)
   double fixedTimeDelay{};
+  // time window for grouping contributions
+  double timeWindow{100 * edm4eic::unit::ns};
 };
 
 } // namespace eicrecon

src/algorithms/digi/SiliconPulseGeneration.cc → src/algorithms/digi/PulseGeneration.cc

@@ -1,22 +1,27 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
-// Copyright (C) 2024-2025 Simon Gardner, Chun Yuen Tsang, Prithwish Tribedy
+// Copyright (C) 2024-2025 Simon Gardner, Chun Yuen Tsang, Prithwish Tribedy,
+//                         Minho Kim, Sylvester Joosten, Wouter Deconinck, Dmitry Kalinkin
 //
 // Convert energy deposition into ADC pulses
 // ADC pulses are assumed to follow the shape of landau function
 
-#include "SiliconPulseGeneration.h"
+#include "PulseGeneration.h"
 
 #include <RtypesCore.h>
 #include <TMath.h>
 #include <algorithms/service.h>
-#include <edm4hep/Vector3d.h>
+#include <edm4hep/CaloHitContribution.h>
+#include <edm4hep/MCParticle.h>
 #include <edm4hep/Vector3f.h>
+#include <podio/RelationRange.h>
+#include <stdlib.h>
+#include <algorithm>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <functional>
-#include <gsl/pointers>
 #include <stdexcept>
+#include <string>
 #include <unordered_map>
 #include <vector>
 
@@ -127,7 +132,37 @@ class PulseShapeFactory {
   }
 };
 
-void SiliconPulseGeneration::init() {
+std::tuple<double, double>
+HitAdapter<edm4hep::SimTrackerHit>::getPulseSources(const edm4hep::SimTrackerHit& hit) {
+  return {hit.getTime(), hit.getEDep()};
+}
+
+#if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
+void HitAdapter<edm4hep::SimTrackerHit>::addRelations(MutablePulseType& pulse,
+                                                      const edm4hep::SimTrackerHit& hit) {
+  pulse.addToTrackerHits(hit);
+  pulse.addToParticles(hit.getParticle());
+}
+#endif
+
+std::tuple<double, double>
+HitAdapter<edm4hep::SimCalorimeterHit>::getPulseSources(const edm4hep::SimCalorimeterHit& hit) {
+  const auto& contribs = hit.getContributions();
+  auto earliest_contrib =
+      std::min_element(contribs.begin(), contribs.end(),
+                       [](const auto& a, const auto& b) { return a.getTime() < b.getTime(); });
+  return {earliest_contrib->getTime(), hit.getEnergy()};
+}
+
+#if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
+void HitAdapter<edm4hep::SimCalorimeterHit>::addRelations(MutablePulseType& pulse,
+                                                          const edm4hep::SimCalorimeterHit& hit) {
+  pulse.addToCalorimeterHits(hit);
+  pulse.addToParticles(hit.getContributions(0).getParticle());
+}
+#endif
+
+template <typename HitT> void PulseGeneration<HitT>::init() {
   m_pulse =
       PulseShapeFactory::createPulseShape(m_cfg.pulse_shape_function, m_cfg.pulse_shape_params);
   m_min_sampling_time = m_cfg.min_sampling_time;
@@ -137,17 +172,15 @@ void SiliconPulseGeneration::init() {
   }
 }
 
-void SiliconPulseGeneration::process(const SiliconPulseGeneration::Input& input,
-                                     const SiliconPulseGeneration::Output& output) const {
+template <typename HitT>
+void PulseGeneration<HitT>::process(
+    const typename PulseGenerationAlgorithm<HitT>::Input& input,
+    const typename PulseGenerationAlgorithm<HitT>::Output& output) const {
   const auto [simhits] = input;
   auto [rawPulses]     = output;
 
   for (const auto& hit : *simhits) {
-
-    auto cellID   = hit.getCellID();
-    double time   = hit.getTime();
-    double charge = hit.getEDep();
-
+    const auto [time, charge] = HitAdapter<HitT>::getPulseSources(hit);
     // Calculate nearest timestep to the hit time rounded down (assume clocks aligned with time 0)
     double signal_time = m_cfg.timestep * std::floor(time / m_cfg.timestep);
 
@@ -178,7 +211,7 @@ void SiliconPulseGeneration::process(const SiliconPulseGeneration::Input& input,
     }
 
     auto time_series = rawPulses->create();
-    time_series.setCellID(cellID);
+    time_series.setCellID(hit.getCellID());
     time_series.setInterval(m_cfg.timestep);
     time_series.setTime(signal_time + skip_bins * m_cfg.timestep);
 
@@ -190,11 +223,13 @@ void SiliconPulseGeneration::process(const SiliconPulseGeneration::Input& input,
     time_series.setIntegral(integral);
     time_series.setPosition(
         edm4hep::Vector3f(hit.getPosition().x, hit.getPosition().y, hit.getPosition().z));
-    time_series.addToTrackerHits(hit);
-    time_series.addToParticles(hit.getParticle());
+    HitAdapter<HitT>::addRelations(time_series, hit);
 #endif
   }
 
-} // SiliconPulseGeneration:process
+} // PulseGeneration:process
+
+template class PulseGeneration<edm4hep::SimTrackerHit>;
+template class PulseGeneration<edm4hep::SimCalorimeterHit>;
 
 } // namespace eicrecon

src/algorithms/digi/PulseGeneration.h

@@ -0,0 +1,76 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024-2025 Simon Gardner, Chun Yuen Tsang, Prithwish Tribedy,
+//                         Minho Kim, Sylvester Joosten, Wouter Deconinck, Dmitry Kalinkin
+//
+// Convert energy deposition into analog pulses
+
+#pragma once
+
+#include <algorithms/algorithm.h>
+#include <edm4eic/EDM4eicVersion.h>
+#include <edm4eic/unit_system.h>
+#include <edm4hep/SimCalorimeterHitCollection.h>
+#include <edm4hep/SimTrackerHitCollection.h>
+#if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
+#include <edm4eic/SimPulseCollection.h>
+#else
+#include <edm4hep/TimeSeriesCollection.h>
+#endif
+#include <memory>
+#include <string_view>
+#include <tuple>
+#include <variant>
+
+#include "algorithms/digi/PulseGenerationConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+#if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
+using PulseType        = edm4eic::SimPulse;
+using MutablePulseType = edm4eic::MutableSimPulse;
+#else
+using PulseType        = edm4hep::TimeSeries;
+using MutablePulseType = edm4hep::MutableTimeSeries;
+#endif
+
+template <typename HitT> struct HitAdapter;
+
+template <> struct HitAdapter<edm4hep::SimTrackerHit> {
+  static std::tuple<double, double> getPulseSources(const edm4hep::SimTrackerHit& hit);
+#if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
+  static void addRelations(MutablePulseType& pulse, const edm4hep::SimTrackerHit& hit);
+#endif
+};
+
+template <> struct HitAdapter<edm4hep::SimCalorimeterHit> {
+  static std::tuple<double, double> getPulseSources(const edm4hep::SimCalorimeterHit& hit);
+#if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
+  static void addRelations(MutablePulseType& pulse, const edm4hep::SimCalorimeterHit& hit);
+#endif
+};
+
+template <typename HitT>
+using PulseGenerationAlgorithm =
+    algorithms::Algorithm<algorithms::Input<typename HitT::collection_type>,
+                          algorithms::Output<PulseType::collection_type>>;
+
+class SignalPulse;
+
+template <typename HitT>
+class PulseGeneration : public PulseGenerationAlgorithm<HitT>,
+                        public WithPodConfig<PulseGenerationConfig> {
+
+public:
+  PulseGeneration(std::string_view name)
+      : PulseGenerationAlgorithm<HitT>{name, {"RawHits"}, {"OutputPulses"}, {}} {}
+  void init() final;
+  void process(const typename PulseGenerationAlgorithm<HitT>::Input&,
+               const typename PulseGenerationAlgorithm<HitT>::Output&) const final;
+
+private:
+  std::shared_ptr<SignalPulse> m_pulse;
+  float m_min_sampling_time = 0 * edm4eic::unit::ns;
+};
+
+} // namespace eicrecon

src/algorithms/digi/SiliconPulseGenerationConfig.h → src/algorithms/digi/PulseGenerationConfig.h

@@ -7,7 +7,7 @@
 
 namespace eicrecon {
 
-struct SiliconPulseGenerationConfig {
+struct PulseGenerationConfig {
   // Parameters of Silicon signal generation
   std::string pulse_shape_function       = "LandauPulse"; // Pulse shape function
   std::vector<double> pulse_shape_params = {1.0, 0.1};    // Parameters of the pulse shape function