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