algorithms/digi/CALOROCDigitization.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Minho Kim, Wouter Deconinck, Dmitry Kalinkin, Derek Anderson, Simon Gardner, Sylvester Joosten, Maria Zurek
0003 //
0004
0005 #include <edm4eic/EDM4eicVersion.h>
0006 #include <tuple>
0007
0008 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 7)
0009 #include <edm4eic/CALOROC1ASample.h>
0010 #include <edm4eic/CALOROC1BSample.h>
0011 #include <podio/RelationRange.h>
0012 #include <algorithm>
0013 #include <cmath>
0014 #include <cstddef>
0015 #include <vector>
0016
0017 #include "CALOROCDigitization.h"
0018
0019 namespace {
0020
0021 struct RawEntry {
0022   double adc{0};
0023   double toa{0};
0024   double tot{0};
0025 };
0026 } // namespace
0027
0028 namespace eicrecon {
0029
0030 void CALOROCDigitization::init() {}
0031
0032 void CALOROCDigitization::process(const CALOROCDigitization::Input& input,
0033                                   const CALOROCDigitization::Output& output) const {
0034   const auto [in_pulses] = input;
0035   auto [out_digi_hits]   = output;
0036
0037   for (const auto& pulse : *in_pulses) {
0038     double pulse_t     = pulse.getTime();
0039     double pulse_dt    = pulse.getInterval();
0040     std::size_t n_amps = pulse.getAmplitude().size();
0041     std::size_t time_stamp =
0042         static_cast<std::size_t>(std::ceil((pulse_t - m_cfg.adc_phase) / m_cfg.time_window));
0043     std::size_t idx_amp_first = static_cast<std::size_t>(
0044         (m_cfg.adc_phase + time_stamp * m_cfg.time_window - pulse_t) / pulse_dt);
0045     std::size_t sample_tick = static_cast<std::size_t>(m_cfg.time_window / pulse_dt);
0046
0047     std::vector<RawEntry> raw_samples(m_cfg.n_samples);
0048
0049     // ADCs are filled in advance because the measurement indices
0050     // are already determined.
0051     // CALOROC measures pulse amplitude for ADC.
0052     for (std::size_t i = 0; i < m_cfg.n_samples; i++) {
0053       std::size_t idx_amp = idx_amp_first + i * sample_tick;
0054       if (idx_amp < n_amps)
0055         raw_samples[i].adc = pulse.getAmplitude()[idx_amp];
0056       else
0057         break;
0058     }
0059
0060     std::size_t idx_sample  = 0;
0061     std::size_t idx_toa     = 0;
0062     double t_upcross        = 0;
0063     bool is_above_threshold = false;
0064
0065     // Measure the TOAs and TOTs while scanning the amplitudes.
0066     // Start from i = 1 since amplitude[i-1] is used to calculate the crossing time.
0067     for (std::size_t i = 1; i < n_amps; i++) {
0068       double t = pulse_t + i * pulse_dt;
0069       if (i > idx_amp_first)
0070         idx_sample = (i + sample_tick - idx_amp_first - 1) / sample_tick;
0071       if (idx_sample == m_cfg.n_samples)
0072         break;
0073
0074       // Measure up-crossing time for TOA
0075       if (!is_above_threshold && pulse.getAmplitude()[i] > m_cfg.toa_thres) {
0076         idx_toa   = idx_sample;
0077         t_upcross = get_crossing_time(m_cfg.toa_thres, pulse_dt, t, pulse.getAmplitude()[i],
0078                                       pulse.getAmplitude()[i - 1]);
0079         raw_samples[idx_toa].toa =
0080             m_cfg.adc_phase + (time_stamp + idx_toa) * m_cfg.time_window - t_upcross;
0081         is_above_threshold = true;
0082       }
0083
0084       // Measure down-crossing time for TOT
0085       if (is_above_threshold && pulse.getAmplitude()[i] < m_cfg.tot_thres) {
0086         raw_samples[idx_toa].tot =
0087             get_crossing_time(m_cfg.tot_thres, pulse_dt, t, pulse.getAmplitude()[i],
0088                               pulse.getAmplitude()[i - 1]) -
0089             t_upcross;
0090         is_above_threshold = false;
0091       }
0092     }
0093
0094     // Fill CALOROCSamples and RawCALOROCHit
0095     auto out_digi_hit = out_digi_hits->create();
0096     out_digi_hit.setCellID(pulse.getCellID());
0097     out_digi_hit.setSamplePhase(std::llround(m_cfg.adc_phase / m_cfg.dyRangeTOA * m_cfg.capTOA));
0098     out_digi_hit.setTimeStamp(time_stamp);
0099
0100     for (const auto& raw_sample : raw_samples) {
0101       auto adc =
0102           std::clamp(std::llround(raw_sample.adc / m_cfg.dyRangeSingleGainADC * m_cfg.capADC), 0LL,
0103                      static_cast<long long>(m_cfg.capADC) - 1);
0104       auto toa = std::clamp(std::llround(raw_sample.toa / m_cfg.dyRangeTOA * m_cfg.capTOA), 0LL,
0105                             static_cast<long long>(m_cfg.capTOA) - 1);
0106       auto tot = std::clamp(std::llround(raw_sample.tot / m_cfg.dyRangeTOT * m_cfg.capTOT), 0LL,
0107                             static_cast<long long>(m_cfg.capTOT) - 1);
0108
0109       out_digi_hit.addToASamples([&]() {
0110         edm4eic::CALOROC1ASample aSample;
0111         aSample.ADC               = adc;
0112         aSample.timeOfArrival     = toa;
0113         aSample.timeOverThreshold = tot;
0114         return aSample;
0115       }());
0116
0117       auto high_adc =
0118           std::clamp(std::llround(raw_sample.adc / m_cfg.dyRangeHighGainADC * m_cfg.capADC), 0LL,
0119                      static_cast<long long>(m_cfg.capADC) - 1);
0120       auto low_adc =
0121           std::clamp(std::llround(raw_sample.adc / m_cfg.dyRangeLowGainADC * m_cfg.capADC), 0LL,
0122                      static_cast<long long>(m_cfg.capADC) - 1);
0123
0124       out_digi_hit.addToBSamples([&]() {
0125         edm4eic::CALOROC1BSample bSample;
0126         bSample.highGainADC   = high_adc;
0127         bSample.lowGainADC    = low_adc;
0128         bSample.timeOfArrival = toa;
0129         return bSample;
0130       }());
0131     }
0132   }
0133 } // CALOROCDigitization:process
0134
0135 double CALOROCDigitization::get_crossing_time(double thres, double dt, double t, double amp1,
0136                                               double amp2) const {
0137   double numerator   = (amp1 - thres) * dt;
0138   double denominator = amp2 - amp1;
0139   double added       = t;
0140   return (numerator / denominator) + added;
0141 }
0142 } // namespace eicrecon
0143 #endif