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 <cmath>
0013 #include <functional>
0014 #include <gsl/pointers>
0015 #include <stdexcept>
0016 #include <unordered_map>
0017 #include <vector>
0018
0019 #include "services/evaluator/EvaluatorSvc.h"
0020
0021 namespace eicrecon {
0022
0023
0024 class SignalPulse {
0025
0026 public:
0027     virtual ~SignalPulse() = default; // Virtual destructor
0028
0029     virtual double operator()(double time, double charge) = 0;
0030
0031     virtual double getMaximumTime() const = 0;
0032
0033 protected:
0034     double m_charge;
0035
0036 };
0037
0038 // ----------------------------------------------------------------------------
0039 // Landau Pulse Shape Functor
0040 // ----------------------------------------------------------------------------
0041 class LandauPulse: public SignalPulse {
0042 public:
0043
0044 LandauPulse(std::vector<double> params) {
0045
0046     if ((params.size() != 2) && (params.size() != 3)) {
0047         throw std::runtime_error("LandauPulse requires 2 or 3 parameters, gain, sigma_analog, [hit_sigma_offset], got " + 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
0058 double operator()(double time, double charge) {
0059     return charge * m_gain * TMath::Landau(time, m_hit_sigma_offset*m_sigma_analog, m_sigma_analog, kTRUE);
0060 }
0061
0062 double getMaximumTime() const {
0063     return m_hit_sigma_offset*m_sigma_analog;
0064 }
0065
0066 private:
0067     double m_gain = 1.0;
0068     double m_sigma_analog = 1.0;
0069     double m_hit_sigma_offset = 3.5;
0070 };
0071
0072 // EvaluatorSvc Pulse
0073 class EvaluatorPulse: public SignalPulse {
0074 public:
0075
0076     EvaluatorPulse(const std::string& expression, const std::vector<double>& params) {
0077
0078         std::vector<std::string> keys = {"time", "charge"};
0079         for(int i=0; i<params.size(); i++) {
0080             std::string p = "param" + std::to_string(i);
0081             //Check the expression contains the parameter
0082             if(expression.find(p) == std::string::npos) {
0083                 throw std::runtime_error("Parameter " + p + " not found in expression");
0084             }
0085             keys.push_back(p);
0086             param_map[p] = params[i];
0087         }
0088
0089         // Check the expression is contains time and charge
0090         if(expression.find("time") == std::string::npos) {
0091             throw std::runtime_error("Parameter [time] not found in expression");
0092         }
0093         if(expression.find("charge") == std::string::npos) {
0094             throw std::runtime_error("Parameter [charge] not found in expression");
0095         }
0096
0097         auto& serviceSvc = algorithms::ServiceSvc::instance();
0098         m_evaluator = serviceSvc.service<EvaluatorSvc>("EvaluatorSvc")->_compile(expression, keys);
0099     };
0100
0101     double operator()(double time, double charge) {
0102         param_map["time"] = time;
0103         param_map["charge"] = charge;
0104         return m_evaluator(param_map);
0105     }
0106
0107     double getMaximumTime() const {
0108         return 0;
0109     }
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
0117
0118
0119 class PulseShapeFactory {
0120 public:
0121     static std::unique_ptr<SignalPulse> createPulseShape(const std::string& type, const std::vector<double>& params) {
0122         if (type == "LandauPulse") {
0123             return std::make_unique<LandauPulse>(params);
0124         }
0125         //
0126         // Add more pulse shape variants here as needed
0127
0128         // If type not found, try and make a function using the ElavulatorSvc
0129         try {
0130             return std::make_unique<EvaluatorPulse>(type,params);
0131         } catch (...) {
0132             throw std::invalid_argument("Unable to make pulse shape type: " + type);
0133         }
0134
0135     }
0136 };
0137
0138
0139 void SiliconPulseGeneration::init() {
0140     m_pulse             = PulseShapeFactory::createPulseShape(m_cfg.pulse_shape_function, m_cfg.pulse_shape_params);
0141     m_min_sampling_time = m_cfg.min_sampling_time;
0142
0143     if(m_pulse->getMaximumTime()>m_min_sampling_time) {
0144       m_min_sampling_time = m_pulse->getMaximumTime();
0145     }
0146 }
0147
0148 void SiliconPulseGeneration::process(const SiliconPulseGeneration::Input& input,
0149                                      const SiliconPulseGeneration::Output& output) const {
0150   const auto [simhits] = input;
0151   auto [rawPulses]     = output;
0152
0153   for (const auto& hit : *simhits) {
0154
0155     auto   cellID = hit.getCellID();
0156     double time   = hit.getTime();
0157     double charge = hit.getEDep();
0158
0159     // Calculate nearest timestep to the hit time rounded down (assume clocks aligned with time 0)
0160     double signal_time = m_cfg.timestep*std::floor(time / m_cfg.timestep);
0161
0162     auto time_series = rawPulses->create();
0163     time_series.setCellID(cellID);
0164     time_series.setInterval(m_cfg.timestep);
0165
0166     bool passed_threshold = false;
0167     int  skip_bins = 0;
0168
0169     for(int i = 0; i < m_cfg.max_time_bins; i ++) {
0170       double t = signal_time + i*m_cfg.timestep - time;
0171       auto signal = (*m_pulse)(t,charge);
0172       if (std::abs(signal) < m_cfg.ignore_thres) {
0173         if(passed_threshold==false) {
0174           skip_bins = i;
0175           continue;
0176         }
0177         if (t > m_min_sampling_time) {
0178           break;
0179         }
0180       }
0181       passed_threshold = true;
0182       time_series.addToAmplitude(signal);
0183     }
0184
0185     time_series.setTime(signal_time+skip_bins*m_cfg.timestep);
0186
0187   }
0188
0189 } // SiliconPulseGeneration:process
0190
0191 } // namespace eicrecon