algorithms/digi/PulseGeneration.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024-2025 Simon Gardner, Chun Yuen Tsang, Prithwish Tribedy,
0003 //                         Minho Kim, Sylvester Joosten, Wouter Deconinck, Dmitry Kalinkin
0004 //
0005 // Convert energy deposition into ADC pulses
0006 // ADC pulses are assumed to follow the shape of landau function
0007
0008 #include "PulseGeneration.h"
0009
0010 #include <RtypesCore.h>
0011 #include <TMath.h>
0012 #include <algorithms/service.h>
0013 #include <edm4hep/CaloHitContribution.h>
0014 #include <edm4hep/MCParticle.h>
0015 #include <edm4hep/Vector3f.h>
0016 #include <podio/RelationRange.h>
0017 #include <stdlib.h>
0018 #include <algorithm>
0019 #include <cmath>
0020 #include <cstddef>
0021 #include <cstdint>
0022 #include <functional>
0023 #include <stdexcept>
0024 #include <string>
0025 #include <unordered_map>
0026 #include <vector>
0027
0028 #include "services/evaluator/EvaluatorSvc.h"
0029
0030 namespace eicrecon {
0031
0032 class SignalPulse {
0033
0034 public:
0035   virtual ~SignalPulse() = default; // Virtual destructor
0036
0037   virtual double operator()(double time, double charge) = 0;
0038
0039   virtual double getMaximumTime() const = 0;
0040 };
0041
0042 // ----------------------------------------------------------------------------
0043 // Landau Pulse Shape Functor
0044 // ----------------------------------------------------------------------------
0045 class LandauPulse : public SignalPulse {
0046 public:
0047   LandauPulse(std::vector<double> params) {
0048
0049     if ((params.size() != 2) && (params.size() != 3)) {
0050       throw std::runtime_error(
0051           "LandauPulse requires 2 or 3 parameters, gain, sigma_analog, [hit_sigma_offset], got " +
0052           std::to_string(params.size()));
0053     }
0054
0055     m_gain         = params[0];
0056     m_sigma_analog = params[1];
0057     if (params.size() == 3) {
0058       m_hit_sigma_offset = params[2];
0059     }
0060   };
0061
0062   double operator()(double time, double charge) override {
0063     return charge * m_gain *
0064            TMath::Landau(time, m_hit_sigma_offset * m_sigma_analog, m_sigma_analog, kTRUE);
0065   }
0066
0067   double getMaximumTime() const override { return m_hit_sigma_offset * m_sigma_analog; }
0068
0069 private:
0070   double m_gain             = 1.0;
0071   double m_sigma_analog     = 1.0;
0072   double m_hit_sigma_offset = 3.5;
0073 };
0074
0075 // EvaluatorSvc Pulse
0076 class EvaluatorPulse : public SignalPulse {
0077 public:
0078   EvaluatorPulse(const std::string& expression, const std::vector<double>& params) {
0079
0080     std::vector<std::string> keys = {"time", "charge"};
0081     for (std::size_t i = 0; i < params.size(); i++) {
0082       std::string p = "param" + std::to_string(i);
0083       //Check the expression contains the parameter
0084       if (expression.find(p) == std::string::npos) {
0085         throw std::runtime_error("Parameter " + p + " not found in expression");
0086       }
0087       keys.push_back(p);
0088       param_map[p] = params[i];
0089     }
0090
0091     // Check the expression is contains time and charge
0092     if (expression.find("time") == std::string::npos) {
0093       throw std::runtime_error("Parameter [time] not found in expression");
0094     }
0095     if (expression.find("charge") == std::string::npos) {
0096       throw std::runtime_error("Parameter [charge] not found in expression");
0097     }
0098
0099     auto& serviceSvc = algorithms::ServiceSvc::instance();
0100     m_evaluator      = serviceSvc.service<EvaluatorSvc>("EvaluatorSvc")->_compile(expression, keys);
0101   };
0102
0103   double operator()(double time, double charge) override {
0104     param_map["time"]   = time;
0105     param_map["charge"] = charge;
0106     return m_evaluator(param_map);
0107   }
0108
0109   double getMaximumTime() const override { return 0; }
0110
0111 private:
0112   std::unordered_map<std::string, double> param_map;
0113   std::function<double(const std::unordered_map<std::string, double>&)> m_evaluator;
0114 };
0115
0116 class PulseShapeFactory {
0117 public:
0118   static std::unique_ptr<SignalPulse> createPulseShape(const std::string& type,
0119                                                        const std::vector<double>& params) {
0120     if (type == "LandauPulse") {
0121       return std::make_unique<LandauPulse>(params);
0122     }
0123     //
0124     // Add more pulse shape variants here as needed
0125
0126     // If type not found, try and make a function using the ElavulatorSvc
0127     try {
0128       return std::make_unique<EvaluatorPulse>(type, params);
0129     } catch (...) {
0130       throw std::invalid_argument("Unable to make pulse shape type: " + type);
0131     }
0132   }
0133 };
0134
0135 std::tuple<double, double>
0136 HitAdapter<edm4hep::SimTrackerHit>::getPulseSources(const edm4hep::SimTrackerHit& hit) {
0137   return {hit.getTime(), hit.getEDep()};
0138 }
0139
0140 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
0141 void HitAdapter<edm4hep::SimTrackerHit>::addRelations(MutablePulseType& pulse,
0142                                                       const edm4hep::SimTrackerHit& hit) {
0143   pulse.addToTrackerHits(hit);
0144   pulse.addToParticles(hit.getParticle());
0145 }
0146 #endif
0147
0148 std::tuple<double, double>
0149 HitAdapter<edm4hep::SimCalorimeterHit>::getPulseSources(const edm4hep::SimCalorimeterHit& hit) {
0150   const auto& contribs = hit.getContributions();
0151   auto earliest_contrib =
0152       std::min_element(contribs.begin(), contribs.end(),
0153                        [](const auto& a, const auto& b) { return a.getTime() < b.getTime(); });
0154   return {earliest_contrib->getTime(), hit.getEnergy()};
0155 }
0156
0157 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
0158 void HitAdapter<edm4hep::SimCalorimeterHit>::addRelations(MutablePulseType& pulse,
0159                                                           const edm4hep::SimCalorimeterHit& hit) {
0160   pulse.addToCalorimeterHits(hit);
0161   pulse.addToParticles(hit.getContributions(0).getParticle());
0162 }
0163 #endif
0164
0165 template <typename HitT> void PulseGeneration<HitT>::init() {
0166   m_pulse =
0167       PulseShapeFactory::createPulseShape(m_cfg.pulse_shape_function, m_cfg.pulse_shape_params);
0168   m_min_sampling_time = m_cfg.min_sampling_time;
0169
0170   if (m_pulse->getMaximumTime() > m_min_sampling_time) {
0171     m_min_sampling_time = m_pulse->getMaximumTime();
0172   }
0173 }
0174
0175 template <typename HitT>
0176 void PulseGeneration<HitT>::process(
0177     const typename PulseGenerationAlgorithm<HitT>::Input& input,
0178     const typename PulseGenerationAlgorithm<HitT>::Output& output) const {
0179   const auto [simhits] = input;
0180   auto [rawPulses]     = output;
0181
0182   for (const auto& hit : *simhits) {
0183     const auto [time, charge] = HitAdapter<HitT>::getPulseSources(hit);
0184     // Calculate nearest timestep to the hit time rounded down (assume clocks aligned with time 0)
0185     double signal_time = m_cfg.timestep * std::floor(time / m_cfg.timestep);
0186
0187     bool passed_threshold   = false;
0188     std::uint32_t skip_bins = 0;
0189     float integral          = 0;
0190     std::vector<float> pulse;
0191
0192     for (std::uint32_t i = 0; i < m_cfg.max_time_bins; i++) {
0193       double t    = signal_time + i * m_cfg.timestep - time;
0194       auto signal = (*m_pulse)(t, charge);
0195       if (std::abs(signal) < m_cfg.ignore_thres) {
0196         if (!passed_threshold) {
0197           skip_bins = i;
0198           continue;
0199         }
0200         if (t > m_min_sampling_time) {
0201           break;
0202         }
0203       }
0204       passed_threshold = true;
0205       pulse.push_back(signal);
0206       integral += signal;
0207     }
0208
0209     if (!passed_threshold) {
0210       continue;
0211     }
0212
0213     auto time_series = rawPulses->create();
0214     time_series.setCellID(hit.getCellID());
0215     time_series.setInterval(m_cfg.timestep);
0216     time_series.setTime(signal_time + skip_bins * m_cfg.timestep);
0217
0218     for (const auto& value : pulse) {
0219       time_series.addToAmplitude(value);
0220     }
0221
0222 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
0223     time_series.setIntegral(integral);
0224     time_series.setPosition(
0225         edm4hep::Vector3f(hit.getPosition().x, hit.getPosition().y, hit.getPosition().z));
0226     HitAdapter<HitT>::addRelations(time_series, hit);
0227 #endif
0228   }
0229
0230 } // PulseGeneration:process
0231
0232 template class PulseGeneration<edm4hep::SimTrackerHit>;
0233 template class PulseGeneration<edm4hep::SimCalorimeterHit>;
0234
0235 } // namespace eicrecon