algorithms/digi/PhotoMultiplierHitDigi.h

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 #pragma once
0016
0017 #include <DD4hep/Detector.h>
0018 #include <DD4hep/Objects.h>
0019 #include <DDRec/CellIDPositionConverter.h>
0020 #include <Math/GenVector/Cartesian3D.h>
0021 #include <Math/GenVector/DisplacementVector3D.h>
0022 #include <algorithms/algorithm.h>
0023 #include <algorithms/geo.h>
0024 #include <edm4eic/MCRecoTrackerHitAssociationCollection.h>
0025 #include <edm4eic/RawTrackerHitCollection.h>
0026 #include <edm4hep/EventHeaderCollection.h>
0027 #include <edm4hep/SimTrackerHitCollection.h>
0028 #include <fmt/format.h> // IWYU pragma: keep
0029 #include <stdint.h>
0030 #include <cstddef>
0031 #include <functional>
0032 #include <gsl/pointers>
0033 #include <iterator>
0034 #include <random>
0035 #include <stdexcept>
0036 #include <string>
0037 #include <string_view>
0038 #include <unordered_map>
0039 #include <utility>
0040 #include <vector>
0041
0042 #include "PhotoMultiplierHitDigiConfig.h"
0043 #include "algorithms/interfaces/UniqueIDGenSvc.h"
0044 #include "algorithms/interfaces/WithPodConfig.h"
0045
0046 namespace eicrecon {
0047
0048 using PhotoMultiplierHitDigiAlgorithm = algorithms::Algorithm<
0049     algorithms::Input<edm4hep::EventHeaderCollection, edm4hep::SimTrackerHitCollection>,
0050     algorithms::Output<edm4eic::RawTrackerHitCollection,
0051                        edm4eic::MCRecoTrackerHitAssociationCollection>>;
0052
0053 class PhotoMultiplierHitDigi : public PhotoMultiplierHitDigiAlgorithm,
0054                                public WithPodConfig<PhotoMultiplierHitDigiConfig> {
0055
0056 public:
0057   PhotoMultiplierHitDigi(std::string_view name)
0058       : PhotoMultiplierHitDigiAlgorithm{name,
0059                                         {"eventHeaderCollection", "inputHitCollection"},
0060                                         {"outputRawHitCollection", "outputRawHitAssociations"},
0061                                         "Digitize within ADC range, add pedestal, convert time "
0062                                         "with smearing resolution."} {}
0063
0064   void init() final;
0065   void process(const Input&, const Output&) const final;
0066
0067   // EDM datatype member types
0068   using CellIDType = decltype(edm4hep::SimTrackerHitData::cellID);
0069   using TimeType   = decltype(edm4hep::SimTrackerHitData::time);
0070
0071   // local structure to hold data for a hit
0072   struct HitData {
0073     uint32_t npe;
0074     double signal;
0075     TimeType time;
0076     std::vector<std::size_t> sim_hit_indices;
0077   };
0078
0079   // set `m_VisitAllRngPixels`, a visitor to run an action (type
0080   // `function<void(cellID)>`) on a selection of random CellIDs; must be
0081   // defined externally, since this would be detector-specific
0082   void SetVisitRngCellIDs(std::function<void(std::function<void(CellIDType)>, float)> visitor) {
0083     m_VisitRngCellIDs = visitor;
0084   }
0085
0086   // set `m_PixelGapMask`, which takes `cellID` and MC hit position, returning
0087   // true if the hit position is on a pixel, or false if on a pixel gap; must be
0088   // defined externally, since this would be detector-specific
0089   void SetPixelGapMask(std::function<bool(CellIDType, dd4hep::Position)> mask) {
0090     m_PixelGapMask = mask;
0091   }
0092
0093 protected:
0094   // visitor of all possible CellIDs (set with SetVisitRngCellIDs)
0095   std::function<void(std::function<void(CellIDType)>, float)> m_VisitRngCellIDs =
0096       [](std::function<void(CellIDType)> /* visitor_action */,
0097          float /* p */) { /* default no-op */ };
0098
0099   // pixel gap mask
0100   std::function<bool(CellIDType, dd4hep::Position)> m_PixelGapMask =
0101       [](CellIDType /* cellID */, dd4hep::Position /* pos_hit_global */) {
0102         throw std::runtime_error("pixel gap cuts enabled, but none defined");
0103         return false;
0104       };
0105
0106 private:
0107   // add a hit to local `hit_groups` data structure
0108   void InsertHit(std::unordered_map<CellIDType, std::vector<HitData>>& hit_groups, CellIDType id,
0109                  double amp, TimeType time, std::size_t sim_hit_index,
0110                  std::default_random_engine& generator, std::normal_distribution<double>& gaussian,
0111                  bool is_noise_hit = false) const;
0112
0113   const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
0114   const dd4hep::rec::CellIDPositionConverter* m_converter{
0115       algorithms::GeoSvc::instance().cellIDPositionConverter()};
0116
0117   const algorithms::UniqueIDGenSvc& m_uid = algorithms::UniqueIDGenSvc::instance();
0118
0119   // std::default_random_engine generator; // TODO: need something more appropriate here
0120   // std::normal_distribution<double> m_normDist; // defaults to mean=0, sigma=1
0121
0122   std::vector<std::pair<double, double>> qeff;
0123   void qe_init();
0124   template <class RndmIter, typename T, class Compare>
0125   RndmIter interval_search(RndmIter beg, RndmIter end, const T& val, Compare comp) const;
0126   bool qe_pass(double ev, double rand) const;
0127 };
0128 } // namespace eicrecon