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