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