algorithms/digi/PhotoMultiplierHitDigi.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022, 2023, Chao Peng, Thomas Britton, Christopher Dilks, Luigi Dello Stritto
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  *  Ported from Juggler by Thomas Britton (JLab)
0013  */
0014
0015
0016 #include "PhotoMultiplierHitDigi.h"
0017
0018 #include <Evaluator/DD4hepUnits.h>
0019 #include <algorithms/logger.h>
0020 #include <edm4eic/EDM4eicVersion.h>
0021 #include <edm4hep/Vector3d.h>
0022 #include <fmt/core.h>
0023 #include <math.h>
0024 #include <podio/ObjectID.h>
0025 #include <algorithm>
0026 #include <gsl/pointers>
0027 #include <iterator>
0028
0029 #include "algorithms/digi/PhotoMultiplierHitDigiConfig.h"
0030
0031 namespace eicrecon {
0032
0033 //------------------------
0034 // init
0035 //------------------------
0036 void PhotoMultiplierHitDigi::init()
0037 {
0038     // print the configuration parameters
0039     debug() << m_cfg << endmsg;
0040
0041     /* warn if using potentially thread-unsafe seed
0042      * FIXME: remove this warning when this issue is resolved:
0043      *        https://github.com/eic/EICrecon/issues/539
0044      */
0045     if(m_cfg.seed==0) warning("using seed=0 may cause thread-unsafe behavior of TRandom (EICrecon issue 539)");
0046
0047     // random number generators
0048     m_random.SetSeed(m_cfg.seed);
0049     m_rngNorm = [&](){
0050         return m_random.Gaus(0., 1.0);
0051     };
0052     m_rngUni = [&](){
0053         return m_random.Uniform(0., 1.0);
0054     };
0055     //auto randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
0056     auto sc1 = m_rngUni;//m_rngUni.initialize(randSvc, Rndm::Flat(0., 1.));
0057     auto sc2 = m_rngNorm;//m_rngNorm.initialize(randSvc, Rndm::Gauss(0., 1.));
0058     //if (!sc1.isSuccess() || !sc2.isSuccess()) {
0059     if (!sc1 || !sc2) {
0060         throw std::runtime_error("Cannot initialize random generator!");
0061     }
0062
0063     // initialize quantum efficiency table
0064     qe_init();
0065 }
0066
0067
0068 //------------------------
0069 // process
0070 //------------------------
0071 void PhotoMultiplierHitDigi::process(
0072       const PhotoMultiplierHitDigi::Input& input,
0073       const PhotoMultiplierHitDigi::Output& output) const
0074 {
0075         const auto [sim_hits] = input;
0076         auto [raw_hits, hit_assocs] = output;
0077
0078         trace("{:=^70}"," call PhotoMultiplierHitDigi::process ");
0079         std::unordered_map<CellIDType, std::vector<HitData>> hit_groups;
0080         // collect the photon hit in the same cell
0081         // calculate signal
0082         trace("{:-<70}","Loop over simulated hits ");
0083         for(std::size_t sim_hit_index = 0; sim_hit_index < sim_hits->size(); sim_hit_index++) {
0084             const auto& sim_hit = sim_hits->at(sim_hit_index);
0085             auto edep_eV = sim_hit.getEDep() * 1e9; // [GeV] -> [eV] // FIXME: use common unit converters, when available
0086             auto id      = sim_hit.getCellID();
0087             trace("hit: pixel id={:#018X}  edep = {} eV", id, edep_eV);
0088
0089             // overall safety factor
0090             if (m_rngUni() > m_cfg.safetyFactor) continue;
0091
0092             // quantum efficiency
0093             if (!qe_pass(edep_eV, m_rngUni())) continue;
0094
0095             // pixel gap cuts
0096             if(m_cfg.enablePixelGaps) {
0097               auto pos = sim_hit.getPosition();
0098               if( ! m_PixelGapMask(id, dd4hep::Position(pos.x*dd4hep::mm, pos.y*dd4hep::mm, pos.z*dd4hep::mm)) )
0099                 continue;
0100             }
0101
0102             // cell time, signal amplitude, truth photon
0103             trace(" -> hit accepted");
0104             trace(" -> MC hit id={}", sim_hit.getObjectID().index);
0105             auto   time = sim_hit.getTime();
0106             double amp  = m_cfg.speMean + m_rngNorm() * m_cfg.speError;
0107
0108             // insert hit to `hit_groups`
0109             InsertHit(
0110                 hit_groups,
0111                 id,
0112                 amp,
0113                 time,
0114                 sim_hit_index
0115                 );
0116         }
0117
0118         // print `hit_groups`
0119         if(level() <= algorithms::LogLevel::kTrace) {
0120           trace("{:-<70}","Accepted hit groups ");
0121           for(auto &[id,hitVec] : hit_groups)
0122             for(auto &hit : hitVec) {
0123               trace("hit_group: pixel id={:#018X} -> npe={} signal={} time={}", id, hit.npe, hit.signal, hit.time);
0124               for(auto i : hit.sim_hit_indices)
0125                 trace(" - MC hit: EDep={}, id={}", sim_hits->at(i).getEDep(), sim_hits->at(i).getObjectID().index);
0126             }
0127         }
0128
0129         //build noise raw hits
0130         if (m_cfg.enableNoise) {
0131           trace("{:=^70}"," BEGIN NOISE INJECTION ");
0132           float p = m_cfg.noiseRate*m_cfg.noiseTimeWindow;
0133           auto cellID_action = [this,&hit_groups] (auto id) {
0134
0135             // cell time, signal amplitude
0136             double   amp  = m_cfg.speMean + m_rngNorm()*m_cfg.speError;
0137             TimeType time = m_cfg.noiseTimeWindow*m_rngUni() / dd4hep::ns;
0138             dd4hep::Position pos_hit_global = m_converter->position(id);
0139
0140             // insert in `hit_groups`, or if the pixel already has a hit, update `npe` and `signal`
0141             this->InsertHit(
0142                 hit_groups,
0143                 id,
0144                 amp,
0145                 time,
0146                 0, // not used
0147                 true
0148                 );
0149
0150           };
0151           m_VisitRngCellIDs(cellID_action, p);
0152         }
0153
0154         // build output `RawTrackerHit` and `MCRecoTrackerHitAssociation` collections
0155         trace("{:-<70}","Digitized raw hits ");
0156         for (auto &it : hit_groups) {
0157             for (auto &data : it.second) {
0158
0159                 // build `RawTrackerHit`
0160                 auto raw_hit = raw_hits->create();
0161                 raw_hit.setCellID(it.first);
0162                 raw_hit.setCharge(    static_cast<decltype(edm4eic::RawTrackerHitData::charge)>    (data.signal)                    );
0163                 raw_hit.setTimeStamp( static_cast<decltype(edm4eic::RawTrackerHitData::timeStamp)> (data.time/m_cfg.timeResolution) );
0164                 trace("raw_hit: cellID={:#018X} -> charge={} timeStamp={}",
0165                     raw_hit.getCellID(),
0166                     raw_hit.getCharge(),
0167                     raw_hit.getTimeStamp()
0168                     );
0169
0170                 // build `MCRecoTrackerHitAssociation` (for non-noise hits only)
0171                 if(!data.sim_hit_indices.empty()) {
0172                   for(auto i : data.sim_hit_indices) {
0173                     auto hit_assoc = hit_assocs->create();
0174                     hit_assoc.setWeight(1.0 / data.sim_hit_indices.size()); // not used
0175                     hit_assoc.setRawHit(raw_hit);
0176 #if EDM4EIC_VERSION_MAJOR >= 6
0177                     hit_assoc.setSimHit(sim_hits->at(i));
0178 #else
0179                     hit_assoc.addToSimHits(sim_hits->at(i));
0180 #endif
0181                   }
0182                 }
0183             }
0184         }
0185 }
0186
0187 void PhotoMultiplierHitDigi::qe_init()
0188 {
0189         // get quantum efficiency table
0190         qeff.clear();
0191         auto hc = dd4hep::h_Planck * dd4hep::c_light / (dd4hep::eV * dd4hep::nm); // [eV*nm]
0192         for(const auto &[wl,qe] : m_cfg.quantumEfficiency) {
0193           qeff.emplace_back( hc / wl, qe ); // convert wavelength [nm] -> energy [eV]
0194         }
0195
0196         // sort quantum efficiency data first
0197         std::sort(qeff.begin(), qeff.end(),
0198             [] (const std::pair<double, double> &v1, const std::pair<double, double> &v2) {
0199                 return v1.first < v2.first;
0200             });
0201
0202         // print the table
0203         debug("{:-^60}"," Quantum Efficiency vs. Energy ");
0204         for(auto& [en,qe] : qeff)
0205           debug("  {:>10.4} {:<}",en,qe);
0206         trace("{:=^60}","");
0207
0208         // sanity checks
0209         if (qeff.empty()) {
0210             qeff = {{2.6, 0.3}, {7.0, 0.3}};
0211             warning("Invalid quantum efficiency data provided, using default values {} {:.2f} {} {:.2f} {} {:.2f} {} {:.2f} {}","{{", qeff.front().first, ",", qeff.front().second, "},{",qeff.back().first,",",qeff.back().second,"}}");
0212         }
0213         if (qeff.front().first > 3.0) {
0214             warning("Quantum efficiency data start from {:.2f} {}", qeff.front().first, " eV, maybe you are using wrong units?");
0215         }
0216         if (qeff.back().first < 3.0) {
0217             warning("Quantum efficiency data end at {:.2f} {}", qeff.back().first, " eV, maybe you are using wrong units?");
0218         }
0219 }
0220
0221
0222 template<class RndmIter, typename T, class Compare> RndmIter PhotoMultiplierHitDigi::interval_search(RndmIter beg, RndmIter end, const T &val, Compare comp) const
0223 {
0224         // special cases
0225         auto dist = std::distance(beg, end);
0226         if ((dist < 2) || (comp(*beg, val) > 0) || (comp(*std::prev(end), val) < 0)) {
0227             return end;
0228         }
0229         auto mid = std::next(beg, dist / 2);
0230
0231         while (mid != end) {
0232             if (comp(*mid, val) == 0) {
0233                 return mid;
0234             } else if (comp(*mid, val) > 0) {
0235                 end = mid;
0236             } else {
0237                 beg = std::next(mid);
0238             }
0239             mid = std::next(beg, std::distance(beg, end)/2);
0240         }
0241
0242         if (mid == end || comp(*mid, val) > 0) {
0243             return std::prev(mid);
0244         }
0245         return mid;
0246 }
0247
0248 bool PhotoMultiplierHitDigi::qe_pass(double ev, double rand) const
0249 {
0250         auto it = interval_search(qeff.begin(), qeff.end(), ev,
0251                     [] (const std::pair<double, double> &vals, double val) {
0252                         return vals.first - val;
0253                     });
0254
0255         if (it == qeff.end()) {
0256             // warning("{} eV is out of QE data range, assuming 0\% efficiency",ev);
0257             return false;
0258         }
0259
0260         double prob = it->second;
0261         auto itn = std::next(it);
0262         if (itn != qeff.end() && (itn->first - it->first != 0)) {
0263             prob = (it->second*(itn->first - ev) + itn->second*(ev - it->first)) / (itn->first - it->first);
0264         }
0265
0266         // trace("{} eV, QE: {}\%",ev,prob*100.);
0267         return rand <= prob;
0268 }
0269
0270
0271 // add a hit to local `hit_groups` data structure
0272 // NOLINTBEGIN(bugprone-easily-swappable-parameters)
0273 void PhotoMultiplierHitDigi::InsertHit(
0274     std::unordered_map<CellIDType, std::vector<HitData>> &hit_groups,
0275     CellIDType       id,
0276     double           amp,
0277     TimeType         time,
0278     std::size_t      sim_hit_index,
0279     bool             is_noise_hit
0280     ) const // NOLINTEND(bugprone-easily-swappable-parameters)
0281 {
0282   auto it = hit_groups.find(id);
0283   if (it != hit_groups.end()) {
0284     std::size_t i = 0;
0285     for (auto ghit = it->second.begin(); ghit != it->second.end(); ++ghit, ++i) {
0286       if (std::abs(time - ghit->time) <= (m_cfg.hitTimeWindow)) {
0287         // hit group found, update npe, signal, and list of MC hits
0288         ghit->npe += 1;
0289         ghit->signal += amp;
0290         if(!is_noise_hit) ghit->sim_hit_indices.push_back(sim_hit_index);
0291         trace(" -> add to group @ {:#018X}: signal={}", id, ghit->signal);
0292         break;
0293       }
0294     }
0295     // no hits group found
0296     if (i >= it->second.size()) {
0297       auto sig = amp + m_cfg.pedMean + m_cfg.pedError * m_rngNorm();
0298       decltype(HitData::sim_hit_indices) indices;
0299       if(!is_noise_hit) indices.push_back(sim_hit_index);
0300       hit_groups.insert({ id, {HitData{1, sig, time, indices}} });
0301       trace(" -> no group found,");
0302       trace("    so new group @ {:#018X}: signal={}", id, sig);
0303     }
0304   } else {
0305     auto sig = amp + m_cfg.pedMean + m_cfg.pedError * m_rngNorm();
0306     decltype(HitData::sim_hit_indices) indices;
0307     if(!is_noise_hit) indices.push_back(sim_hit_index);
0308     hit_groups.insert({ id, {HitData{1, sig, time, indices}} });
0309     trace(" -> new group @ {:#018X}: signal={}", id, sig);
0310   }
0311 }
0312
0313 } // namespace eicrecon