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