src/todo/PhotoMultiplierDigi.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng
0003
0004 /*  General PhotoMultiplier Digitization
0005  *
0006  *  Apply the given quantum efficiency for photon detection
0007  *  Converts the number of detected photons to signal amplitude
0008  *
0009  *  Author: Chao Peng (ANL)
0010  *  Date: 10/02/2020
0011  */
0012
0013 #include <iterator>
0014 #include <algorithm>
0015 #include <unordered_map>
0016 #include <cmath>
0017
0018 #include "GaudiAlg/Transformer.h"
0019 #include "GaudiAlg/GaudiTool.h"
0020 #include "GaudiKernel/RndmGenerators.h"
0021 #include "GaudiKernel/PhysicalConstants.h"
0022
0023 #include "JugBase/DataHandle.h"
0024
0025 // Event Model related classes
0026 #include "edm4eic/RawTrackerHitCollection.h"
0027 #include "edm4hep/MCParticleCollection.h"
0028 #include "edm4hep/SimTrackerHitCollection.h"
0029
0030
0031 using namespace Gaudi::Units;
0032
0033 namespace Jug::Digi {
0034
0035 /** PhotoMultiplierDigi.
0036  *
0037  * \ingroup digi
0038  */
0039 class PhotoMultiplierDigi : public GaudiAlgorithm
0040 {
0041 public:
0042     DataHandle<edm4hep::SimTrackerHitCollection>
0043         m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0044     DataHandle<edm4eic::RawTrackerHitCollection>
0045         m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
0046     Gaudi::Property<std::vector<std::pair<double, double>>>
0047         u_quantumEfficiency{this, "quantumEfficiency", {{2.6*eV, 0.3}, {7.0*eV, 0.3}}};
0048     Gaudi::Property<double> m_hitTimeWindow{this, "hitTimeWindow", 20.0*ns};
0049     Gaudi::Property<double> m_timeStep{this, "timeStep", 0.0625*ns};
0050     Gaudi::Property<double> m_speMean{this, "speMean", 80.0};
0051     Gaudi::Property<double> m_speError{this, "speError", 16.0};
0052     Gaudi::Property<double> m_pedMean{this, "pedMean", 200.0};
0053     Gaudi::Property<double> m_pedError{this, "pedError", 3.0};
0054     Rndm::Numbers m_rngUni, m_rngNorm;
0055
0056     // constructor
0057     PhotoMultiplierDigi(const std::string& name, ISvcLocator* svcLoc)
0058         : GaudiAlgorithm(name, svcLoc)
0059     {
0060         declareProperty("inputHitCollection", m_inputHitCollection,"");
0061         declareProperty("outputHitCollection", m_outputHitCollection, "");
0062     }
0063
0064     StatusCode initialize() override
0065     {
0066         if (GaudiAlgorithm::initialize().isFailure()) {
0067             return StatusCode::FAILURE;
0068         }
0069
0070         auto randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
0071         auto sc1 = m_rngUni.initialize(randSvc, Rndm::Flat(0., 1.));
0072         auto sc2 = m_rngNorm.initialize(randSvc, Rndm::Gauss(0., 1.));
0073         if (!sc1.isSuccess() || !sc2.isSuccess()) {
0074             error() << "Cannot initialize random generator!" << endmsg;
0075             return StatusCode::FAILURE;
0076         }
0077
0078         qe_init();
0079
0080         return StatusCode::SUCCESS;
0081     }
0082
0083     StatusCode execute() override
0084     {
0085         // input collection
0086         const auto &sim = *m_inputHitCollection.get();
0087         // Create output collections
0088         auto &raw = *m_outputHitCollection.createAndPut();
0089
0090         struct HitData { int npe; double signal; double time; };
0091         std::unordered_map<decltype(edm4eic::RawTrackerHitData::cellID), std::vector<HitData>> hit_groups;
0092         // collect the photon hit in the same cell
0093         // calculate signal
0094         for(const auto& ahit : sim) {
0095             // quantum efficiency
0096             if (!qe_pass(ahit.getEDep(), m_rngUni())) {
0097                 continue;
0098             }
0099             // cell id, time, signal amplitude
0100             uint64_t id = ahit.getCellID();
0101             double time = ahit.getMCParticle().getTime();
0102             double amp = m_speMean + m_rngNorm()*m_speError;
0103
0104             // group hits
0105             auto it = hit_groups.find(id);
0106             if (it != hit_groups.end()) {
0107                 size_t i = 0;
0108                 for (auto git = it->second.begin(); git != it->second.end(); ++git, ++i) {
0109                     if (std::abs(time - git->time) <= (m_hitTimeWindow/ns)) {
0110                         git->npe += 1;
0111                         git->signal += amp;
0112                         break;
0113                     }
0114                 }
0115                 // no hits group found
0116                 if (i >= it->second.size()) {
0117                     it->second.emplace_back(HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time});
0118                 }
0119             } else {
0120                 hit_groups[id] = {HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time}};
0121             }
0122         }
0123
0124         // build hit
0125         for (auto &it : hit_groups) {
0126             for (auto &data : it.second) {
0127                 edm4eic::RawTrackerHit hit{
0128                   it.first,
0129                   static_cast<decltype(edm4eic::RawTrackerHitData::charge)>(data.signal),
0130                   static_cast<decltype(edm4eic::RawTrackerHitData::timeStamp)>(data.time/(m_timeStep/ns))};
0131                 raw.push_back(hit);
0132             }
0133         }
0134
0135         return StatusCode::SUCCESS;
0136     }
0137
0138 private:
0139     void qe_init()
0140     {
0141         auto &qeff = u_quantumEfficiency.value();
0142
0143         // sort quantum efficiency data first
0144         std::sort(qeff.begin(), qeff.end(),
0145             [] (const std::pair<double, double> &v1, const std::pair<double, double> &v2) {
0146                 return v1.first < v2.first;
0147             });
0148
0149         // sanity checks
0150         if (qeff.empty()) {
0151             qeff = {{2.6*eV, 0.3}, {7.0*eV, 0.3}};
0152             warning() << "Invalid quantum efficiency data provided, using default values: " << qeff << endmsg;
0153         }
0154         if (qeff.front().first > 3.0*eV) {
0155             warning() << "Quantum efficiency data start from " << qeff.front().first/eV
0156                       << " eV, maybe you are using wrong units?" << endmsg;
0157         }
0158         if (qeff.back().first < 6.0*eV) {
0159             warning() << "Quantum efficiency data end at " << qeff.back().first/eV
0160                       << " eV, maybe you are using wrong units?" << endmsg;
0161         }
0162     }
0163
0164     // helper function for linear interpolation
0165     // Comp return is defined as: equal, 0;  greater, > 0; less, < 0
0166     template<class RndmIter, typename T, class Compare>
0167     RndmIter interval_search(RndmIter beg, RndmIter end, const T &val, Compare comp) const
0168     {
0169         // special cases
0170         auto dist = std::distance(beg, end);
0171         if ((dist < 2) || (comp(*beg, val) > 0) || (comp(*std::prev(end), val) < 0)) {
0172             return end;
0173         }
0174         auto mid = std::next(beg, dist / 2);
0175
0176         while (mid != end) {
0177             if (comp(*mid, val) == 0) {
0178                 return mid;
0179             } else if (comp(*mid, val) > 0) {
0180                 end = mid;
0181             } else {
0182                 beg = std::next(mid);
0183             }
0184             mid = std::next(beg, std::distance(beg, end)/2);
0185         }
0186
0187         if (mid == end || comp(*mid, val) > 0) {
0188             return std::prev(mid);
0189         }
0190         return mid;
0191     }
0192
0193     bool qe_pass(double ev, double rand) const
0194     {
0195         const auto &qeff = u_quantumEfficiency.value();
0196         auto it = interval_search(qeff.begin(), qeff.end(), ev,
0197                     [] (const std::pair<double, double> &vals, double val) {
0198                         return vals.first - val;
0199                     });
0200
0201         if (it == qeff.end()) {
0202             // info() << ev/eV << " eV is out of QE data range, assuming 0% efficiency" << endmsg;
0203             return false;
0204         }
0205
0206         double prob = it->second;
0207         auto itn = std::next(it);
0208         if (itn != qeff.end() && (itn->first - it->first != 0)) {
0209             prob = (it->second*(itn->first - ev) + itn->second*(ev - it->first)) / (itn->first - it->first);
0210         }
0211
0212         // info() << ev/eV << " eV, QE: "  << prob*100. << "%" << endmsg;
0213         return rand <= prob;
0214     }
0215 };
0216
0217 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0218 DECLARE_COMPONENT(PhotoMultiplierDigi)
0219
0220 } // namespace Jug::Digi