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.empty()) {
0031       continue;
0032     }
0033     int n_CFDROC_cycle = static_cast<int>(std::floor(pulse.getTime() / m_cfg.tMax));
0034
0035     // first we find all the peaks and store their location
0036     // Then we find the time corresponding to fraction of the peak height
0037     // use stack to store peak height and time
0038     std::stack<std::pair<int, int>> peakTimeAndHeight;
0039
0040     for (size_t time_bin = 1; time_bin < adcs.size() - 1; ++time_bin) {
0041       auto V      = adcs[time_bin];
0042       auto prev_V = adcs[time_bin - 1];
0043       auto next_V = adcs[time_bin + 1];
0044       if ((std::abs(prev_V) < std::abs(V)) &&
0045           (std::abs(V) >= std::abs(next_V))) { // To get peak of the Analog signal
0046         peakTimeAndHeight.emplace(time_bin, V);
0047       }
0048     }
0049
0050     // scan the peak in reverse time to find TDC values at fraction of peaks height
0051     // start from the last adc bin
0052     int time_bin = static_cast<int>(adcs.size() - 2);
0053     // find time corresponding to each peak one by one
0054     while (!peakTimeAndHeight.empty()) {
0055       auto peak = peakTimeAndHeight.top();
0056       if (peak.first >= time_bin) {
0057         // peaks that are situated later than current time are all discarded
0058         peakTimeAndHeight.pop();
0059         continue;
0060       }
0061       time_bin            = peak.first;
0062       int target_height_V = static_cast<int>(peak.second * m_cfg.fraction);
0063       int prev_V          = adcs[time_bin];
0064       --time_bin;
0065       // scan the peak in reverse time to find TDC values at fraction of peaks height
0066       for (; time_bin >= 0; --time_bin) {
0067         double V = adcs[time_bin];
0068         // check voltage threshold
0069         if (std::abs(V) <= std::abs(target_height_V) &&
0070             std::abs(target_height_V) <= std::abs(prev_V)) {
0071           int tdc = time_bin + n_CFDROC_cycle * m_cfg.tdc_range;
0072           // limit the range of adc values
0073           int adc = std::min(m_cfg.adc_range, std::abs(peak.second));
0074           rawhits->create(pulse.getCellID(), adc, tdc);
0075           // the peak is found. Break the time scan to the next peak
0076           break;
0077         }
0078         prev_V = V;
0079       }
0080       // discard current peak from stack to look for the next peak
0081       peakTimeAndHeight.pop();
0082     }
0083   }
0084 } // CFDROCDigitization:process
0085 } // namespace eicrecon