rich/include/Tools.h

0001 // Copyright 2023, Christopher Dilks
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003
0004 #pragma once
0005
0006 // general
0007 #include <map>
0008 #include <spdlog/spdlog.h>
0009
0010 // ROOT
0011 #include <TVector3.h>
0012
0013 // IRT
0014 #include <IRT/ParametricSurface.h>
0015
0016 // DD4hep
0017 #include <Evaluator/DD4hepUnits.h>
0018
0019 // data model
0020 #include <edm4hep/EDM4hepVersion.h>
0021 #include <edm4hep/SimTrackerHitCollection.h>
0022 #include <edm4hep/utils/kinematics.h>
0023
0024 namespace benchmarks {
0025
0026   // Tools class, filled with miscellanous helper functions
0027   class Tools {
0028     public:
0029
0030       // -------------------------------------------------------------------------------------
0031       // physics constants
0032
0033       // h*c constant, for wavelength <=> energy conversion [GeV*nm]
0034       // obtained from `dd4hep::h_Planck * dd4hep::c_light / (dd4hep::GeV * dd4hep::nm)`
0035       static constexpr double HC = 1.239841875e-06;
0036
0037       // opticalphoton PDG
0038       static constexpr int PHOTON_PDG = -22;
0039
0040       // -------------------------------------------------------------------------------------
0041       // PDG mass lookup
0042
0043       // local PDG mass database
0044       // FIXME: cannot use `TDatabasePDG` since it is not thread safe; until we
0045       // have a proper PDG database service, we hard-code the masses we need;
0046       // use Tools::GetPDGMass for access
0047       static std::unordered_map<int,double> GetPDGMasses() {
0048         return std::unordered_map<int,double>{
0049           { 11,   0.000510999 },
0050           { 211,  0.13957     },
0051           { 321,  0.493677    },
0052           { 2212, 0.938272    }
0053         };
0054       }
0055
0056       static double GetPDGMass(int pdg) {
0057         double mass;
0058         try { mass = GetPDGMasses().at(std::abs(pdg)); }
0059         catch(const std::out_of_range& e) {
0060           throw std::runtime_error(fmt::format("RUNTIME ERROR: unknown PDG={} in algorithms/pid/Tools::GetPDGMass",pdg));
0061         }
0062         return mass;
0063       }
0064
0065       static int GetNumPDGs() { return GetPDGMasses().size(); };
0066
0067
0068       // -------------------------------------------------------------------------------------
0069       // get photon wavelength
0070       static double GetPhotonWavelength(const edm4hep::SimTrackerHit& hit, std::shared_ptr<spdlog::logger> log) {
0071 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
0072         auto phot = hit.getParticle();
0073 #else
0074         auto phot = hit.getMCParticle();
0075 #endif
0076         if(!phot.isAvailable()) {
0077           log->error("no MCParticle in hit");
0078           return -1.0;
0079         }
0080         if(phot.getPDG() != PHOTON_PDG) {
0081           log->warn("hit MCParticle is not an opticalphoton; PDG is {}; ignoring", phot.getPDG());
0082           return -1.0;
0083         }
0084         auto momentum   = edm4hep::utils::p(phot);
0085         auto wavelength = momentum > 0 ? HC / momentum : 0.0;
0086         return wavelength;
0087       }
0088
0089       // -------------------------------------------------------------------------------------
0090       // binning
0091       static constexpr int    n_bins         = 100;
0092       static constexpr int    momentum_bins  = 500;
0093       static constexpr double momentum_max   = 70;
0094       static constexpr int    eta_bins       = 50;
0095       static constexpr double eta_min        = 0;
0096       static constexpr double eta_max        = 5;
0097       static constexpr int    npe_max        = 50;
0098       static constexpr int    nphot_max      = 3000;
0099       static constexpr int    theta_bins     = 1500;
0100       static constexpr double theta_max      = 300;
0101       static constexpr double thetaResid_max = 100;
0102       static constexpr int    phi_bins       = 100;
0103       static constexpr int    wl_bins        = 120;
0104       static constexpr double wl_min         = 0;
0105       static constexpr double wl_max         = 1200;
0106
0107   }; // class Tools
0108 } // namespace benchmarks