algorithms/pid/Tools.h

0001 // Copyright 2023, Alexander Kiselev, Christopher Dilks
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 //
0004 // Common functions for PID algorithms
0005 // Several methods ported from Juggler's JugPID `IRTAlgorithmServices`
0006 //
0007
0008 #pragma once
0009
0010 // general
0011 #include <map>
0012 #include <math.h>
0013 #include <spdlog/spdlog.h>
0014
0015 // ROOT
0016 #include <TVector3.h>
0017
0018 // DD4hep
0019 #include <Evaluator/DD4hepUnits.h>
0020
0021 // data model
0022 #include <edm4eic/CherenkovParticleIDCollection.h>
0023 #include <edm4hep/ParticleIDCollection.h>
0024
0025
0026 namespace eicrecon {
0027
0028   // Tools class, filled with miscellaneous helper functions
0029   class Tools {
0030     public:
0031
0032       // -------------------------------------------------------------------------------------
0033
0034       // h*c constant, for wavelength <=> energy conversion [GeV*nm]
0035       static constexpr double HC = dd4hep::h_Planck * dd4hep::c_light / (dd4hep::GeV * dd4hep::nm);
0036
0037
0038       // -------------------------------------------------------------------------------------
0039       // Radiator IDs
0040
0041       static std::unordered_map<int,std::string> GetRadiatorIDs() {
0042         return std::unordered_map<int,std::string>{
0043           { 0, "Aerogel" },
0044           { 1, "Gas" }
0045         };
0046       }
0047
0048       static std::string GetRadiatorName(int id) {
0049         std::string name;
0050         try { name = GetRadiatorIDs().at(id); }
0051         catch(const std::out_of_range& e) {
0052           throw std::runtime_error(fmt::format("RUNTIME ERROR: unknown radiator ID={} in algorithms/pid/Tools::GetRadiatorName",id));
0053         }
0054         return name;
0055       }
0056
0057       static int GetRadiatorID(std::string name) {
0058         for(auto& [id,name_tmp] : GetRadiatorIDs())
0059           if(name==name_tmp) return id;
0060         throw std::runtime_error(fmt::format("RUNTIME ERROR: unknown radiator '{}' in algorithms/pid/Tools::GetRadiatorID",name));
0061         return -1;
0062       }
0063
0064
0065       // -------------------------------------------------------------------------------------
0066       // Table rebinning and lookup
0067
0068       // Rebin input table `input` to have `nbins+1` equidistant bins; returns the rebinned table
0069       static std::vector<std::pair<double, double>> ApplyFineBinning(
0070           const std::vector<std::pair<double,double>> &input,
0071           unsigned nbins
0072           )
0073       {
0074         std::vector<std::pair<double, double>> ret;
0075
0076         // Well, could have probably just reordered the initial vector;
0077         std::map<double, double> buffer;
0078
0079         for(auto entry: input)
0080           buffer[entry.first] = entry.second;
0081
0082         // Sanity checks; return empty map in case do not pass them;
0083         if (buffer.size() < 2 || nbins < 2) return ret;
0084
0085         double from = (*buffer.begin()).first;
0086         double to   = (*buffer.rbegin()).first;
0087         // Will be "nbins+1" equidistant entries;
0088         double step = (to - from) / nbins;
0089
0090         for(auto entry: buffer) {
0091           double e1 = entry.first;
0092           double qe1 = entry.second;
0093
0094           if (!ret.size())
0095             ret.push_back(std::make_pair(e1, qe1));
0096           else {
0097             const auto &prev = ret[ret.size()-1];
0098
0099             double e0 = prev.first;
0100             double qe0 = prev.second;
0101             double a = (qe1 - qe0) / (e1 - e0);
0102             double b = qe0 - a*e0;
0103             // FIXME: check floating point accuracy when moving to a next point; do we actually
0104             // care whether the overall number of bins will be "nbins+1" or more?;
0105             for(double e = e0+step; e<e1; e+=step)
0106               ret.push_back(std::make_pair(e, a*e + b));
0107           } //if
0108         } //for entry
0109
0110         return ret;
0111       }
0112
0113       // Find the bin in table `table` that contains entry `argument` in the first column and
0114       // sets `entry` to the corresponding element of the second column; returns true if successful
0115       static bool GetFinelyBinnedTableEntry(
0116           const std::vector<std::pair<double, double>> &table,
0117           double argument,
0118           double *entry
0119           )
0120       {
0121         // Get the tabulated table reference; perform sanity checks;
0122         //const std::vector<std::pair<double, double>> &qe = u_quantumEfficiency.value();
0123         unsigned dim = table.size(); if (dim < 2) return false;
0124
0125         // Find a proper bin; no tricks, they are all equidistant;
0126         auto const &from = table[0];
0127         auto const &to = table[dim-1];
0128         double emin = from.first;
0129         double emax = to.first;
0130         double step = (emax - emin) / (dim - 1);
0131         int ibin = (int)floor((argument - emin) / step);
0132
0133         //printf("%f vs %f, %f -> %d\n", ev, from.first, to. first, ibin);
0134
0135         // Out of range check;
0136         if (ibin < 0 || ibin >= int(dim)) return false;
0137
0138         *entry = table[ibin].second;
0139         return true;
0140       }
0141
0142       // -------------------------------------------------------------------------------------
0143       // convert PODIO vector datatype to ROOT TVector3
0144       template<class PodioVector3>
0145         static TVector3 PodioVector3_to_TVector3(const PodioVector3 v) {
0146           return TVector3(v.x, v.y, v.z);
0147         }
0148       // convert ROOT::Math::Vector to ROOT TVector3
0149       template<class MathVector3>
0150         static TVector3 MathVector3_to_TVector3(MathVector3 v) {
0151           return TVector3(v.x(), v.y(), v.z());
0152         }
0153
0154       // -------------------------------------------------------------------------------------
0155
0156       // printing: vectors
0157       static void PrintTVector3(
0158           std::shared_ptr<spdlog::logger> m_log,
0159           std::string name, TVector3 vec,
0160           int nameBuffer=30, spdlog::level::level_enum lvl=spdlog::level::trace
0161           )
0162       {
0163         m_log->log(lvl, "{:>{}} = ( {:>10.2f} {:>10.2f} {:>10.2f} )", name, nameBuffer, vec.x(), vec.y(), vec.z());
0164       }
0165
0166       // printing: hypothesis tables
0167       static void PrintHypothesisTableHead(
0168           std::shared_ptr<spdlog::logger> m_log,
0169           int indent=4, spdlog::level::level_enum lvl=spdlog::level::trace
0170           )
0171       {
0172         m_log->log(lvl, "{:{}}{:>6}  {:>10}  {:>10}", "", indent, "PDG", "Weight", "NPE");
0173       }
0174       static void PrintHypothesisTableLine(
0175           std::shared_ptr<spdlog::logger> m_log,
0176           edm4eic::CherenkovParticleIDHypothesis hyp,
0177           int indent=4, spdlog::level::level_enum lvl=spdlog::level::trace
0178           )
0179       {
0180         m_log->log(lvl, "{:{}}{:>6}  {:>10.8}  {:>10.8}", "", indent, hyp.PDG, hyp.weight, hyp.npe);
0181       }
0182       static void PrintHypothesisTableLine(
0183           std::shared_ptr<spdlog::logger> m_log,
0184           edm4hep::ParticleID hyp,
0185           int indent=4, spdlog::level::level_enum lvl=spdlog::level::trace
0186           )
0187       {
0188         float npe = hyp.parameters_size() > 0 ? hyp.getParameters(0) : -1; // assume NPE is the first parameter
0189         m_log->log(lvl, "{:{}}{:>6}  {:>10.8}  {:>10.8}", "", indent, hyp.getPDG(), hyp.getLikelihood(), npe);
0190       }
0191
0192       // printing: Cherenkov angle estimate
0193       static void PrintCherenkovEstimate(
0194           std::shared_ptr<spdlog::logger> m_log,
0195           edm4eic::CherenkovParticleID pid,
0196           bool printHypotheses = true,
0197           int indent=2, spdlog::level::level_enum lvl=spdlog::level::trace
0198           )
0199       {
0200         if(m_log->level() <= lvl) {
0201           double thetaAve = 0;
0202           if(pid.getNpe() > 0)
0203             for(const auto& [theta,phi] : pid.getThetaPhiPhotons())
0204               thetaAve += theta / pid.getNpe();
0205           m_log->log(lvl, "{:{}}Cherenkov Angle Estimate:", "", indent);
0206           m_log->log(lvl, "{:{}}  {:>16}:  {:>10}",         "", indent, "NPE",      pid.getNpe());
0207           m_log->log(lvl, "{:{}}  {:>16}:  {:>10.8} mrad",  "", indent, "<theta>",  thetaAve*1e3); // [rad] -> [mrad]
0208           m_log->log(lvl, "{:{}}  {:>16}:  {:>10.8}",       "", indent, "<rindex>", pid.getRefractiveIndex());
0209           m_log->log(lvl, "{:{}}  {:>16}:  {:>10.8} eV",    "", indent, "<energy>", pid.getPhotonEnergy()*1e9); // [GeV] -> [eV]
0210           if(printHypotheses) {
0211             m_log->log(lvl, "{:{}}Mass Hypotheses:",          "", indent);
0212             Tools::PrintHypothesisTableHead(m_log, indent+2);
0213             for(const auto& hyp : pid.getHypotheses())
0214               Tools::PrintHypothesisTableLine(m_log, hyp, indent+2);
0215           }
0216         }
0217       }
0218
0219
0220   }; // class Tools
0221 } // namespace eicrecon