algorithms/digi/SiliconPulseDiscretization.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Chun Yuen Tsang
0003 //
0004 // Creates a discrete pulse from a continuous pulse
0005 //
0006
0007 #include <DDRec/CellIDPositionConverter.h>
0008 #include <podio/RelationRange.h>
0009 #include <cmath>
0010 #include <gsl/pointers>
0011 #include <unordered_map>
0012 #include <unordered_set>
0013 #include <utility>
0014
0015 #include "SiliconPulseDiscretization.h"
0016 // use TGraph for interpolation
0017 #include "TGraph.h"
0018
0019 namespace eicrecon {
0020
0021 void SiliconPulseDiscretization::init() {}
0022
0023 double SiliconPulseDiscretization::_interpolateOrZero(const TGraph& graph, double t, double tMin,
0024                                                       double tMax) const {
0025   // return 0 if t is outside of tMin - tMax range
0026   // otherwise, return graph interpolation value
0027   if (t < tMin || t > tMax) {
0028     return 0;
0029   }
0030   double height = graph.Eval(t, nullptr, "S"); // spline interpolation
0031   if (!std::isfinite(height)) {
0032     error("Pulse interpolation returns nan. This happen mostly because there are multiple "
0033           "pulse height values at the same time. Did you call PulseCombiner?");
0034   }
0035   return height;
0036 }
0037
0038 void SiliconPulseDiscretization::process(const SiliconPulseDiscretization::Input& input,
0039                                          const SiliconPulseDiscretization::Output& output) const {
0040   const auto [inPulses] = input;
0041   auto [outPulses]      = output;
0042
0043   // sometimes the first hit arrives early, but the last hit arrive very late
0044   // there is a lot of nothing in between
0045   // If we loop through those time interval with nothing in it, the creation of outPulses will take forever
0046   // Speeds things up by denoting which EICROCcycle contains pulse information
0047   // And only focus on those cycles
0048   std::unordered_map<dd4hep::rec::CellID, std::pair<TGraph, std::unordered_set<int>>>
0049       graphAndCycle4Cells;
0050
0051   for (const auto& pulse : *inPulses) {
0052
0053     auto cellID   = pulse.getCellID();
0054     auto time     = pulse.getTime();
0055     auto interval = pulse.getInterval();
0056
0057     // one TGraph per pulse
0058     // Interpolate the pulse with TGraph
0059     auto& graph        = graphAndCycle4Cells[cellID].first;
0060     auto& activeCycles = graphAndCycle4Cells[cellID].second;
0061     for (unsigned int i = 0; i < pulse.getAmplitude().size(); i++) {
0062       auto currTime = time + i * interval;
0063       // current EICROC cycle
0064       int EICROCCycle = std::floor(currTime / m_cfg.EICROC_period);
0065       activeCycles.insert(EICROCCycle);
0066       graph.SetPoint(graph.GetN(), currTime + m_cfg.global_offset, pulse.getAmplitude()[i]);
0067     }
0068   }
0069
0070   // sort all pulses data points to avoid TGraph::Eval giving nan due to non-monotonic data
0071   for (auto& [_, graphAndCycle] : graphAndCycle4Cells) {
0072     graphAndCycle.first.Sort();
0073     graphAndCycle.first.SetBit(TGraph::kIsSortedX);
0074   }
0075
0076   // sum all digitized pulses
0077   for (const auto& [cellID, graphAndCycle] : graphAndCycle4Cells) {
0078     const auto& graph       = graphAndCycle.first;
0079     const auto& activeCycle = graphAndCycle.second;
0080     double tMin             = NAN;
0081     double tMax             = NAN;
0082     double temp             = NAN; // do not use
0083     graph.ComputeRange(tMin, temp, tMax, temp);
0084
0085     for (int curriEICRocCycle : activeCycle) {
0086       // time beings at an EICROC cycle
0087       double startTime = curriEICRocCycle * m_cfg.EICROC_period;
0088
0089       auto outPulse = outPulses->create();
0090       outPulse.setCellID(cellID);
0091       outPulse.setInterval(m_cfg.local_period);
0092       outPulse.setTime(startTime);
0093
0094       // stop at the next cycle
0095       for (double currTime = startTime; currTime < startTime + m_cfg.EICROC_period;
0096            currTime += m_cfg.local_period)
0097         outPulse.addToAdcCounts(this->_interpolateOrZero(graph, currTime, tMin, tMax));
0098     }
0099   }
0100 } // SiliconPulseDiscretization:process
0101 } // namespace eicrecon