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