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 <limits.h>
0009 #include <podio/RelationRange.h>
0010 #include <cmath>
0011 #include <gsl/pointers>
0012 #include <unordered_map>
0013
0014 #include "SiliconPulseDiscretization.h"
0015 // use TGraph for interpolation
0016 #include "TGraph.h"
0017
0018 namespace eicrecon {
0019
0020 void SiliconPulseDiscretization::init() {}
0021
0022 double SiliconPulseDiscretization::_interpolateOrZero(const TGraph& graph, double t, double tMin,
0023                                                       double tMax) const {
0024   // return 0 if t is outside of tMin - tMax range
0025   // otherwise, return graph interpolation value
0026   if (t < tMin || t > tMax)
0027     return 0;
0028   else {
0029     double height = graph.Eval(t, nullptr, "S"); // spline interpolation
0030     if (!std::isfinite(height))
0031       error("Pulse interpolation returns nan. This happen mostly because there are multiple pulse "
0032             "height values at the same time. Did you call PulseCombiner?");
0033     return height;
0034   }
0035 }
0036
0037 void SiliconPulseDiscretization::process(const SiliconPulseDiscretization::Input& input,
0038                                          const SiliconPulseDiscretization::Output& output) const {
0039   const auto [inPulses] = input;
0040   auto [outPulses]      = output;
0041
0042   std::unordered_map<dd4hep::rec::CellID, TGraph> Graph4Cells;
0043
0044   for (const auto& pulse : *inPulses) {
0045
0046     auto cellID   = pulse.getCellID();
0047     auto time     = pulse.getTime();
0048     auto interval = pulse.getInterval();
0049
0050     // one TGraph per pulse
0051     // Interpolate the pulse with TGraph
0052     auto& graph = Graph4Cells[cellID];
0053     for (unsigned int i = 0; i < pulse.getAmplitude().size(); i++) {
0054       auto currTime = time + i * interval;
0055       graph.SetPoint(graph.GetN(), currTime + m_cfg.global_offset, pulse.getAmplitude()[i]);
0056     }
0057   }
0058
0059   // sort all pulses data points to avoid TGraph::Eval giving nan due to non-monotonic data
0060   for (auto& [cellID, graph] : Graph4Cells)
0061     graph.Sort();
0062
0063   // sum all digitized pulses
0064   for (const auto& [cellID, graph] : Graph4Cells) {
0065     double tMin, tMax;
0066     double temp; // do not use
0067     graph.ComputeRange(tMin, temp, tMax, temp);
0068
0069     // time beings at an EICROC cycle
0070     double currTime = std::floor(tMin / m_cfg.EICROC_period) * m_cfg.EICROC_period;
0071     // ensure that the outPulse -> create is called in the first cycle
0072     int iEICRocCycle = INT_MIN;
0073
0074     edm4hep::MutableRawTimeSeries outPulse;
0075     for (; currTime <= tMax; currTime += m_cfg.local_period) {
0076       // find current EICROC cycle NO to see if we arrive at the next cycle
0077       int curriEICRocCycle = std::floor(currTime / m_cfg.EICROC_period);
0078       if (curriEICRocCycle != iEICRocCycle) {
0079         // new pulse for each EICROC cycle
0080         iEICRocCycle = curriEICRocCycle;
0081         outPulse     = outPulses->create();
0082         outPulse.setCellID(cellID);
0083         outPulse.setInterval(m_cfg.local_period);
0084         outPulse.setTime(iEICRocCycle * m_cfg.EICROC_period);
0085       }
0086       outPulse.addToAdcCounts(this->_interpolateOrZero(graph, currTime, tMin, tMax));
0087     }
0088   }
0089 } // SiliconPulseDiscretization:process
0090 } // namespace eicrecon