algorithms/digi/SiliconPulseGeneration.cc

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