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 "Gaudi/Algorithm.h"
0019 #include "GaudiKernel/RndmGenerators.h"
0020 #include "GaudiKernel/PhysicalConstants.h"
0021
0022 #include <k4FWCore/DataHandle.h>
0023
0024 // Event Model related classes
0025 #include "edm4eic/RawTrackerHitCollection.h"
0026 #include "edm4hep/EDM4hepVersion.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 Gaudi::Algorithm
0040 {
0041 public:
0042     mutable DataHandle<edm4hep::SimTrackerHitCollection>
0043         m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0044     mutable 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         : Gaudi::Algorithm(name, svcLoc)
0059     {
0060         declareProperty("inputHitCollection", m_inputHitCollection,"");
0061         declareProperty("outputHitCollection", m_outputHitCollection, "");
0062     }
0063
0064     StatusCode initialize() override
0065     {
0066         if (Gaudi::Algorithm::initialize().isFailure()) {
0067             return StatusCode::FAILURE;
0068         }
0069
0070         auto randSvc = Gaudi::svcLocator()->service<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(const EventContext&) const 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 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
0102             double time = ahit.getParticle().getTime();
0103 #else
0104             double time = ahit.getMCParticle().getTime();
0105 #endif
0106             double amp = m_speMean + m_rngNorm()*m_speError;
0107
0108             // group hits
0109             auto it = hit_groups.find(id);
0110             if (it != hit_groups.end()) {
0111                 size_t i = 0;
0112                 for (auto git = it->second.begin(); git != it->second.end(); ++git, ++i) {
0113                     if (std::abs(time - git->time) <= (m_hitTimeWindow/ns)) {
0114                         git->npe += 1;
0115                         git->signal += amp;
0116                         break;
0117                     }
0118                 }
0119                 // no hits group found
0120                 if (i >= it->second.size()) {
0121                     it->second.emplace_back(HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time});
0122                 }
0123             } else {
0124                 hit_groups[id] = {HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time}};
0125             }
0126         }
0127
0128         // build hit
0129         for (auto &it : hit_groups) {
0130             for (auto &data : it.second) {
0131                 raw.create(
0132                   it.first,
0133                   static_cast<decltype(edm4eic::RawTrackerHitData::charge)>(data.signal),
0134                   static_cast<decltype(edm4eic::RawTrackerHitData::timeStamp)>(data.time/(m_timeStep/ns))
0135                 );
0136             }
0137         }
0138
0139         return StatusCode::SUCCESS;
0140     }
0141
0142 private:
0143     void qe_init()
0144     {
0145         auto &qeff = u_quantumEfficiency.value();
0146
0147         // sort quantum efficiency data first
0148         std::sort(qeff.begin(), qeff.end(),
0149             [] (const std::pair<double, double> &v1, const std::pair<double, double> &v2) {
0150                 return v1.first < v2.first;
0151             });
0152
0153         // sanity checks
0154         if (qeff.empty()) {
0155             qeff = {{2.6*eV, 0.3}, {7.0*eV, 0.3}};
0156             warning() << "Invalid quantum efficiency data provided, using default values: " << qeff << endmsg;
0157         }
0158         if (qeff.front().first > 3.0*eV) {
0159             warning() << "Quantum efficiency data start from " << qeff.front().first/eV
0160                       << " eV, maybe you are using wrong units?" << endmsg;
0161         }
0162         if (qeff.back().first < 6.0*eV) {
0163             warning() << "Quantum efficiency data end at " << qeff.back().first/eV
0164                       << " eV, maybe you are using wrong units?" << endmsg;
0165         }
0166     }
0167
0168     // helper function for linear interpolation
0169     // Comp return is defined as: equal, 0;  greater, > 0; less, < 0
0170     template<class RndmIter, typename T, class Compare>
0171     RndmIter interval_search(RndmIter beg, RndmIter end, const T &val, Compare comp) const
0172     {
0173         // special cases
0174         auto dist = std::distance(beg, end);
0175         if ((dist < 2) || (comp(*beg, val) > 0) || (comp(*std::prev(end), val) < 0)) {
0176             return end;
0177         }
0178         auto mid = std::next(beg, dist / 2);
0179
0180         while (mid != end) {
0181             if (comp(*mid, val) == 0) {
0182                 return mid;
0183             } else if (comp(*mid, val) > 0) {
0184                 end = mid;
0185             } else {
0186                 beg = std::next(mid);
0187             }
0188             mid = std::next(beg, std::distance(beg, end)/2);
0189         }
0190
0191         if (mid == end || comp(*mid, val) > 0) {
0192             return std::prev(mid);
0193         }
0194         return mid;
0195     }
0196
0197     bool qe_pass(double ev, double rand) const
0198     {
0199         const auto &qeff = u_quantumEfficiency.value();
0200         auto it = interval_search(qeff.begin(), qeff.end(), ev,
0201                     [] (const std::pair<double, double> &vals, double val) {
0202                         return vals.first - val;
0203                     });
0204
0205         if (it == qeff.end()) {
0206             // info() << ev/eV << " eV is out of QE data range, assuming 0% efficiency" << endmsg;
0207             return false;
0208         }
0209
0210         double prob = it->second;
0211         auto itn = std::next(it);
0212         if (itn != qeff.end() && (itn->first - it->first != 0)) {
0213             prob = (it->second*(itn->first - ev) + itn->second*(ev - it->first)) / (itn->first - it->first);
0214         }
0215
0216         // info() << ev/eV << " eV, QE: "  << prob*100. << "%" << endmsg;
0217         return rand <= prob;
0218     }
0219 };
0220
0221 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0222 DECLARE_COMPONENT(PhotoMultiplierDigi)
0223
0224 } // namespace Jug::Digi