algorithms/digi/CFDROCDigitization.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Chun Yuen Tsang
0003 //
0004 // Convert ADC pulses into ADC and TDC values using constant fraction
0005
0006 #include <podio/RelationRange.h>
0007 #include <algorithm>
0008 #include <cmath>
0009 #include <cstdlib>
0010 #include <gsl/pointers>
0011 #include <stack>
0012 #include <utility>
0013
0014 #include "CFDROCDigitization.h"
0015 #include "algorithms/digi/CFDROCDigitizationConfig.h"
0016
0017 namespace eicrecon {
0018
0019 void CFDROCDigitization::process(const CFDROCDigitization::Input& input,
0020                                  const CFDROCDigitization::Output& output) const {
0021   const auto [simhits] = input;
0022   auto [rawhits]       = output;
0023
0024   // The real CFD compares delayed pulse with an inverted scaled pulse
0025   // This code is doing none of that, it's simply finding pulse height at a fraction of peak
0026   // more sophisticaed algorithm TBD
0027   //
0028   for (const auto& pulse : *simhits) {
0029     auto adcs = pulse.getAdcCounts();
0030     if (adcs.size() == 0)
0031       continue;
0032     int n_CFDROC_cycle = static_cast<int>(std::floor(pulse.getTime() / m_cfg.tMax));
0033
0034     // first we find all the peaks and store their location
0035     // Then we find the time corresponding to fraction of the peak height
0036     // use stack to store peak height and time
0037     std::stack<std::pair<int, int>> peakTimeAndHeight;
0038
0039     for (size_t time_bin = 1; time_bin < adcs.size() - 1; ++time_bin) {
0040       auto V      = adcs[time_bin];
0041       auto prev_V = adcs[time_bin - 1];
0042       auto next_V = adcs[time_bin + 1];
0043       if ((std::abs(prev_V) < std::abs(V)) &&
0044           (std::abs(V) >= std::abs(next_V))) { // To get peak of the Analog signal
0045         peakTimeAndHeight.push({time_bin, V});
0046       }
0047     }
0048
0049     // scan the peak in reverse time to find TDC values at fraction of peaks height
0050     // start from the last adc bin
0051     int time_bin = static_cast<int>(adcs.size() - 2);
0052     // find time corresponding to each peak one by one
0053     while (!peakTimeAndHeight.empty()) {
0054       auto peak = peakTimeAndHeight.top();
0055       if (peak.first >= time_bin) {
0056         // peaks that are situated later than current time are all discarded
0057         peakTimeAndHeight.pop();
0058         continue;
0059       }
0060       time_bin            = peak.first;
0061       int target_height_V = static_cast<int>(peak.second * m_cfg.fraction);
0062       int prev_V          = adcs[time_bin];
0063       --time_bin;
0064       // scan the peak in reverse time to find TDC values at fraction of peaks height
0065       for (; time_bin >= 0; --time_bin) {
0066         double V = adcs[time_bin];
0067         // check voltage threshold
0068         if (std::abs(V) <= std::abs(target_height_V) &&
0069             std::abs(target_height_V) <= std::abs(prev_V)) {
0070           int tdc = time_bin + n_CFDROC_cycle * m_cfg.tdc_range;
0071           // limit the range of adc values
0072           int adc = std::min(m_cfg.adc_range, std::abs(peak.second));
0073           rawhits->create(pulse.getCellID(), adc, tdc);
0074           // the peak is found. Break the time scan to the next peak
0075           break;
0076         }
0077         prev_V = V;
0078       }
0079       // discard current peak from stack to look for the next peak
0080       peakTimeAndHeight.pop();
0081     }
0082   }
0083 } // CFDROCDigitization:process
0084 } // namespace eicrecon