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