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 <algorithms/logger.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 namespace eicrecon {
0026
0027 // Tools namespace, filled with miscellaneous helper functions
0028 namespace Tools {
0029
0030   // -------------------------------------------------------------------------------------
0031   // Radiator IDs
0032
0033   inline std::unordered_map<int, std::string> GetRadiatorIDs() {
0034     return std::unordered_map<int, std::string>{{0, "Aerogel"}, {1, "Gas"}};
0035   }
0036
0037   inline std::string GetRadiatorName(int id) {
0038     std::string name;
0039     try {
0040       name = GetRadiatorIDs().at(id);
0041     } catch (const std::out_of_range& e) {
0042       throw std::runtime_error(fmt::format(
0043           "RUNTIME ERROR: unknown radiator ID={} in algorithms/pid/Tools::GetRadiatorName", id));
0044     }
0045     return name;
0046   }
0047
0048   inline int GetRadiatorID(std::string name) {
0049     for (auto& [id, name_tmp] : GetRadiatorIDs())
0050       if (name == name_tmp)
0051         return id;
0052     throw std::runtime_error(fmt::format(
0053         "RUNTIME ERROR: unknown radiator '{}' in algorithms/pid/Tools::GetRadiatorID", name));
0054     return -1;
0055   }
0056
0057   // -------------------------------------------------------------------------------------
0058   // Table rebinning and lookup
0059
0060   // Rebin input table `input` to have `nbins+1` equidistant bins; returns the rebinned table
0061   inline std::vector<std::pair<double, double>>
0062   ApplyFineBinning(const std::vector<std::pair<double, double>>& input, unsigned nbins) {
0063     std::vector<std::pair<double, double>> ret;
0064
0065     // Well, could have probably just reordered the initial vector;
0066     std::map<double, double> buffer;
0067
0068     for (auto entry : input)
0069       buffer[entry.first] = entry.second;
0070
0071     // Sanity checks; return empty map in case do not pass them;
0072     if (buffer.size() < 2 || nbins < 2)
0073       return ret;
0074
0075     double from = (*buffer.begin()).first;
0076     double to   = (*buffer.rbegin()).first;
0077     // Will be "nbins+1" equidistant entries;
0078     double step = (to - from) / nbins;
0079
0080     for (auto entry : buffer) {
0081       double e1  = entry.first;
0082       double qe1 = entry.second;
0083
0084       if (!ret.size())
0085         ret.push_back(std::make_pair(e1, qe1));
0086       else {
0087         const auto& prev = ret[ret.size() - 1];
0088
0089         double e0  = prev.first;
0090         double qe0 = prev.second;
0091         double a   = (qe1 - qe0) / (e1 - e0);
0092         double b   = qe0 - a * e0;
0093         // FIXME: check floating point accuracy when moving to a next point; do we actually
0094         // care whether the overall number of bins will be "nbins+1" or more?;
0095         for (double e = e0 + step; e < e1; e += step)
0096           ret.push_back(std::make_pair(e, a * e + b));
0097       } //if
0098     }   //for entry
0099
0100     return ret;
0101   }
0102
0103   // Find the bin in table `table` that contains entry `argument` in the first column and
0104   // sets `entry` to the corresponding element of the second column; returns true if successful
0105   inline bool GetFinelyBinnedTableEntry(const std::vector<std::pair<double, double>>& table,
0106                                         double argument, double* entry) {
0107     // Get the tabulated table reference; perform sanity checks;
0108     //const std::vector<std::pair<double, double>> &qe = u_quantumEfficiency.value();
0109     unsigned dim = table.size();
0110     if (dim < 2)
0111       return false;
0112
0113     // Find a proper bin; no tricks, they are all equidistant;
0114     auto const& from = table[0];
0115     auto const& to   = table[dim - 1];
0116     double emin      = from.first;
0117     double emax      = to.first;
0118     double step      = (emax - emin) / (dim - 1);
0119     int ibin         = (int)floor((argument - emin) / step);
0120
0121     //printf("%f vs %f, %f -> %d\n", ev, from.first, to. first, ibin);
0122
0123     // Out of range check;
0124     if (ibin < 0 || ibin >= int(dim))
0125       return false;
0126
0127     *entry = table[ibin].second;
0128     return true;
0129   }
0130
0131   // -------------------------------------------------------------------------------------
0132   // convert PODIO vector datatype to ROOT TVector3
0133   template <class PodioVector3> TVector3 PodioVector3_to_TVector3(const PodioVector3 v) {
0134     return TVector3(v.x, v.y, v.z);
0135   }
0136   // convert ROOT::Math::Vector to ROOT TVector3
0137   template <class MathVector3> TVector3 MathVector3_to_TVector3(MathVector3 v) {
0138     return TVector3(v.x(), v.y(), v.z());
0139   }
0140
0141   // -------------------------------------------------------------------------------------
0142
0143   // printing: vectors
0144   inline std::string PrintTVector3(std::string name, TVector3 vec, int nameBuffer = 30) {
0145     return fmt::format("{:>{}} = ( {:>10.2f} {:>10.2f} {:>10.2f} )", name, nameBuffer, vec.x(),
0146                        vec.y(), vec.z());
0147   }
0148
0149   // printing: hypothesis tables
0150   inline std::string HypothesisTableHead(int indent = 4) {
0151     return fmt::format("{:{}}{:>6}  {:>10}  {:>10}", "", indent, "PDG", "Weight", "NPE");
0152   }
0153   inline std::string HypothesisTableLine(edm4eic::CherenkovParticleIDHypothesis hyp,
0154                                          int indent = 4) {
0155     return fmt::format("{:{}}{:>6}  {:>10.8}  {:>10.8}", "", indent, hyp.PDG, hyp.weight, hyp.npe);
0156   }
0157   inline std::string HypothesisTableLine(edm4hep::ParticleID hyp, int indent = 4) {
0158     float npe =
0159         hyp.parameters_size() > 0 ? hyp.getParameters(0) : -1; // assume NPE is the first parameter
0160     return fmt::format("{:{}}{:>6}  {:>10.8}  {:>10.8}", "", indent, hyp.getPDG(),
0161                        hyp.getLikelihood(), npe);
0162   }
0163
0164 }; // namespace Tools
0165
0166 } // namespace eicrecon