Digitization algorithm that reproduces the HGCROC measurement

CMakeLists.txt

@@ -58,7 +58,7 @@ endif()
 set(Acts_VERSION_MIN 34.0.0)
 set(algorithms_VERSION_MIN 1.0.0)
 set(DD4hep_VERSION_MIN 1.21)
-set(EDM4EIC_VERSION_MIN 8.0)
+set(EDM4EIC_VERSION_MIN 8.3.0)
 set(EDM4HEP_VERSION_MIN 0.7.1)
 set(Eigen3_VERSION_MIN 3.3)
 set(FastJet_VERSION_MIN 3)

src/algorithms/digi/PulseDigi.cc

@@ -0,0 +1,146 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Minho Kim
+//
+
+#include <algorithms/service.h>
+#include <fmt/core.h>
+#include <podio/RelationRange.h>
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <gsl/pointers>
+#include <limits>
+#include <map>
+#include <random>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include "PulseDigi.h"
+
+class HGCROCRawSample {
+public:
+  struct RawEntry {
+    double adc{0};
+    double toa{0};
+    double tot{0};
+    bool totProgress{false};
+    bool totComplete{false};
+  };
+
+  HGCROCRawSample(std::size_t n_samp) { rawEntries.resize(n_samp); }
+
+  void setAmplitude(std::size_t idx, double amp) { rawEntries[idx].adc = amp; }
+  void setTOA(std::size_t idx, double toa) { rawEntries[idx].toa = toa; }
+  void setTOT(std::size_t idx, double cross_t) {
+    rawEntries[idx].tot = cross_t - rawEntries[idx].toa;
+  }
+  void setTotProgress(std::size_t idx) { rawEntries[idx].totProgress = true; }
+  void setTotComplete(std::size_t idx) { rawEntries[idx].totComplete = true; }
+
+  const std::vector<RawEntry>& getEntries() const { return rawEntries; }
+
+private:
+  std::vector<RawEntry> rawEntries;
+};
+
+namespace eicrecon {
+
+void PulseDigi::init() {}
+
+void PulseDigi::process(const PulseDigi::Input& input, const PulseDigi::Output& output) const {
+  const auto [in_pulses] = input;
+  auto [out_digi_hits]   = output;
+
+  for (const auto& pulse : *in_pulses) {
+    double pulse_t     = pulse.getTime();
+    double pulse_dt    = pulse.getInterval();
+    std::size_t n_amps = pulse.getAmplitude().size();
+
+    // Estimate the number of samples.
+    const std::size_t timeIdx_begin =
+        static_cast<std::size_t>(std::floor(pulse_t / m_cfg.time_window));
+    const std::size_t timeIdx_end = static_cast<std::size_t>(
+        std::floor((pulse_t + (n_amps - 1) * pulse_dt) / m_cfg.time_window));
+
+    HGCROCRawSample raw_sample(timeIdx_end - timeIdx_begin + 1);
+
+    // For ADC, amplitude is measured with a fixed phase.
+    // This was reproduced by sample_tick and adc_counter as follows.
+    // Amplitude is measured whenever adc_counter reaches sample_tick.
+    int sample_tick = std::llround(m_cfg.time_window / pulse_dt);
+    int adc_counter =
+        std::llround((timeIdx_begin * m_cfg.time_window + m_cfg.adc_phase - pulse_t) / pulse_dt);
+
+    bool tot_progress  = false;
+    bool tot_complete  = false;
+    std::size_t toaIdx = 0;
+
+    for (std::size_t i = 0; i < n_amps; i++) {
+      double t = pulse_t + i * pulse_dt;
+      const std::size_t sampleIdx =
+          static_cast<std::size_t>(std::floor(t / m_cfg.time_window)) - timeIdx_begin;
+
+      adc_counter++;
+
+      // Measure amplitudes for ADC
+      if (adc_counter == sample_tick) {
+        raw_sample.setAmplitude(sampleIdx, pulse.getAmplitude()[i]);
+        adc_counter = 0;
+        if (tot_progress)
+          raw_sample.setTotProgress(sampleIdx);
+      }
+
+      // Measure up-crossing time for TOA
+      if (!tot_progress && pulse.getAmplitude()[i] > m_cfg.toa_thres) {
+        toaIdx = sampleIdx;
+        raw_sample.setTOA(sampleIdx,
+                          get_crossing_time(m_cfg.toa_thres, pulse_dt, t, pulse.getAmplitude()[i],
+                                            pulse.getAmplitude()[i - 1]));
+        tot_progress = true;
+        tot_complete = false;
+        raw_sample.setTotProgress(sampleIdx);
+      }
+
+      // Measure down-crossing time for TOT
+      if (tot_progress && !tot_complete && pulse.getAmplitude()[i] < m_cfg.tot_thres) {
+        raw_sample.setTOT(toaIdx,
+                          get_crossing_time(m_cfg.tot_thres, pulse_dt, t, pulse.getAmplitude()[i],
+                                            pulse.getAmplitude()[i - 1]));
+        tot_progress = false;
+        tot_complete = true;
+        raw_sample.setTotComplete(sampleIdx);
+      }
+    }
+
+    // Fill HGCROCSamples and RawHGCROCHit
+    auto out_digi_hit = out_digi_hits->create();
+    out_digi_hit.setCellID(pulse.getCellID());
+
+    const auto& entries = raw_sample.getEntries();
+
+    for (const auto& entry : entries) {
+      edm4eic::HGCROCSample sample;
+      auto adc   = std::max(std::llround(entry.adc / m_cfg.dyRangeADC * m_cfg.capADC), 0LL);
+      sample.ADC = adc > m_cfg.capADC ? m_cfg.capADC : adc;
+      auto toa   = std::max(std::llround(entry.toa / m_cfg.dyRangeTOA * m_cfg.capTOA), 0LL);
+      sample.timeOfArrival = toa > m_cfg.capTOA ? m_cfg.capTOA : toa;
+      auto tot = std::max(std::llround(entry.tot / m_cfg.dyRangeTOT * m_cfg.capTOT), 0LL);
+      sample.timeOverThreshold = tot > m_cfg.capTOT ? m_cfg.capTOT : tot;
+      sample.TOTInProgress     = entry.totProgress;
+      sample.TOTComplete       = entry.totComplete;
+      out_digi_hit.addToSamples(sample);
+    }
+  }
+} // PulseDigi:process
+
+double PulseDigi::get_crossing_time(double thres, double dt, double t, double amp1,
+                                    double amp2) const {
+  double numerator   = (thres - amp2) * dt;
+  double denominator = amp2 - amp1;
+  double added       = t;
+  return (numerator / denominator) + added;
+}
+} // namespace eicrecon

src/algorithms/digi/PulseDigi.h

@@ -0,0 +1,40 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Minho Kim
+
+#pragma once
+
+#include <algorithms/algorithm.h>
+#include <edm4eic/SimPulseCollection.h>
+#include <edm4eic/RawHGCROCHitCollection.h>
+#include <edm4eic/HGCROCSample.h>
+#include <stdint.h>
+#include <functional>
+#include <string>
+#include <string_view>
+
+#include "algorithms/digi/PulseDigiConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using PulseDigiAlgorithm =
+    algorithms::Algorithm<algorithms::Input<edm4eic::SimPulseCollection>,
+                          algorithms::Output<edm4eic::RawHGCROCHitCollection>>;
+
+class PulseDigi : public PulseDigiAlgorithm, public WithPodConfig<PulseDigiConfig> {
+
+public:
+  PulseDigi(std::string_view name)
+      : PulseDigiAlgorithm{name,
+                           {"InputPulses"},
+                           {"OutputDigiHits"},
+                           {"ADC, TOA, and TOT are measured referring to the"
+                            "working principle of the HGCROC."}} {}
+  virtual void init() final;
+  void process(const Input&, const Output&) const;
+
+private:
+  double get_crossing_time(double thres, double dt, double t, double amp1, double amp2) const;
+};
+
+} // namespace eicrecon

src/algorithms/digi/PulseDigiConfig.h

@@ -0,0 +1,27 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Minho Kim
+
+#pragma once
+
+#include <edm4eic/unit_system.h>
+
+namespace eicrecon {
+
+struct PulseDigiConfig {
+
+  // Variables for HGCROC measurement
+  double time_window{25 * edm4eic::unit::ns};
+  double adc_phase{10 * edm4eic::unit::ns};
+  double toa_thres{1};
+  double tot_thres{1};
+
+  // Variables for digitization
+  unsigned int capADC{1024};
+  double dyRangeADC{1};
+  unsigned int capTOA{1024};
+  double dyRangeTOA{1 * edm4eic::unit::ns};
+  unsigned int capTOT{1024};
+  double dyRangeTOT{1 * edm4eic::unit::ns};
+};
+
+} // namespace eicrecon

src/factories/digi/PulseDigi_factory.h

@@ -0,0 +1,48 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Minho Kim
+
+#pragma once
+
+#include "algorithms/digi/PulseDigi.h"
+#include "services/algorithms_init/AlgorithmsInit_service.h"
+#include "extensions/jana/JOmniFactory.h"
+
+namespace eicrecon {
+
+class PulseDigi_factory : public JOmniFactory<PulseDigi_factory, PulseDigiConfig> {
+
+public:
+  using AlgoT = eicrecon::PulseDigi;
+
+private:
+  std::unique_ptr<AlgoT> m_algo;
+
+  PodioInput<edm4eic::SimPulse> m_pulse_input{this};
+  PodioOutput<edm4eic::RawHGCROCHit> m_digi_output{this};
+
+  ParameterRef<double> m_time_window{this, "timeWindow", config().time_window};
+  ParameterRef<double> m_adc_phase{this, "adcPhase", config().adc_phase};
+  ParameterRef<double> m_toa_thres{this, "toaThres", config().toa_thres};
+  ParameterRef<double> m_tot_thres{this, "totThres", config().tot_thres};
+  ParameterRef<unsigned int> m_capADC{this, "capADC", config().capADC};
+  ParameterRef<double> m_dyRangeADC{this, "dyRangeADC", config().dyRangeADC};
+  ParameterRef<unsigned int> m_capTOA{this, "capTOA", config().capTOA};
+  ParameterRef<double> m_dyRangeTOA{this, "dyRangeTOA", config().dyRangeTOA};
+  ParameterRef<unsigned int> m_capTOT{this, "capTOT", config().capTOT};
+  ParameterRef<double> m_dyRangeTOT{this, "dyRangeTOT", config().dyRangeTOT};
+
+  Service<AlgorithmsInit_service> m_algorithmsInit{this};
+
+public:
+  void Configure() {
+    m_algo = std::make_unique<AlgoT>(GetPrefix());
+    m_algo->level(static_cast<algorithms::LogLevel>(logger()->level()));
+    m_algo->applyConfig(config());
+    m_algo->init();
+  }
+
+  void Process(int32_t /* run_number */, uint64_t /* event_number */) {
+    m_algo->process({m_pulse_input()}, {m_digi_output().get()});
+  }
+};
+} // namespace eicrecon