Belle II Software light-2406-ragdoll
PIDDetectorWeights.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <framework/gearbox/Unit.h>
10#include <analysis/dbobjects/PIDDetectorWeights.h>
11
12using namespace Belle2;
13
14
16{
17
18 auto dummyBinEdges = std::set<double>();
19
20 // Flatten the RDataFrame content into the tabular structure for each std charged hypothesis.
21 for (const auto& hypo : Const::chargedStableSet) {
22
23 WeightsTable weightsTable;
24
25 auto cut_pdgId = [&](double pdg) { return (pdg == hypo.getPDGCode()); };
26 auto filteredRDF = m_weightsRDataFrame.Filter(cut_pdgId, {"pdgId"});
27
28 if (*filteredRDF.Count()) {
29
30 // Get the lower and upper bin edges from the RDF as C++ vectors
31 auto pMinEdges = filteredRDF.Take<double>("p_min").GetValue();
32 auto pMaxEdges = filteredRDF.Take<double>("p_max").GetValue();
33 // Get the bin indexes.
34 auto pBinIdxs = filteredRDF.Take<double>("p_bin_idx").GetValue();
35 // Convert to a set and take the union. This way duplicates are removed.
36 std::set<double> pMinEdges_set(std::begin(pMinEdges), std::end(pMinEdges));
37 std::set<double> pMaxEdges_set(std::begin(pMaxEdges), std::end(pMaxEdges));
38 std::set<double> pBinEdges;
39 std::set_union(std::begin(pMinEdges_set), std::end(pMinEdges_set),
40 std::begin(pMaxEdges_set), std::end(pMaxEdges_set),
41 std::inserter(pBinEdges, std::begin(pBinEdges)));
42 // Store the set of bin edges into the internal tabular structure for this particle hypothesis.
43 weightsTable.m_pBinEdges = pBinEdges;
44 // Store the nr. of p bins. Used for bin idx linearisation.
45 weightsTable.m_nPBins = pBinEdges.size() - 1;
46 // Check thqt the min, max bin edges have no gaps.
47 // Revert to a std::vector for easy index-based iteration.
48 std::vector<double> pMinEdgesVec(std::begin(pMinEdges_set), std::end(pMinEdges_set));
49 std::vector<double> pMaxEdgesVec(std::begin(pMaxEdges_set), std::end(pMaxEdges_set));
50 for (unsigned int iEdgeIdx(0); iEdgeIdx < pMaxEdgesVec.size() - 1; ++iEdgeIdx) {
51 if (pMaxEdgesVec[iEdgeIdx] != pMinEdgesVec[iEdgeIdx + 1]) {
52 B2FATAL("The p bins in the weights table are not contiguous. Please check the input CSV file.");
53 }
54 }
55
56 auto thetaMinEdges = filteredRDF.Take<double>("theta_min").GetValue();
57 auto thetaMaxEdges = filteredRDF.Take<double>("theta_max").GetValue();
58 auto thetaBinIdxs = filteredRDF.Take<double>("theta_bin_idx").GetValue();
59 std::set<double> thetaMinEdges_set(std::begin(thetaMinEdges), std::end(thetaMinEdges));
60 std::set<double> thetaMaxEdges_set(std::begin(thetaMaxEdges), std::end(thetaMaxEdges));
61 std::set<double> thetaBinEdges;
62 std::set_union(std::begin(thetaMinEdges_set), std::end(thetaMinEdges_set),
63 std::begin(thetaMaxEdges_set), std::end(thetaMaxEdges_set),
64 std::inserter(thetaBinEdges, std::begin(thetaBinEdges)));
65 weightsTable.m_thetaBinEdges = thetaBinEdges;
66 weightsTable.m_nThetaBins = thetaBinEdges.size() - 1;
67 std::vector<double> thetaMinEdgesVec(std::begin(thetaMinEdges_set), std::end(thetaMinEdges_set));
68 std::vector<double> thetaMaxEdgesVec(std::begin(thetaMaxEdges_set), std::end(thetaMaxEdges_set));
69 for (unsigned int iEdge(0); iEdge < thetaMaxEdgesVec.size() - 1; ++iEdge) {
70 if (thetaMaxEdgesVec[iEdge] != thetaMinEdgesVec[iEdge + 1]) {
71 B2FATAL("The theta bins in the weights table are not contiguous. Please check the input CSV file.");
72 }
73 }
74
75 // Get the p, theta bin index columns and zip them together.
76 std::vector<std::tuple<double, double>> pAndThetaIdxs;
77 std::transform(std::begin(pBinIdxs), std::end(pBinIdxs),
78 std::begin(thetaBinIdxs),
79 std::back_inserter(pAndThetaIdxs), [](const auto & pIdx, const auto & thetaIdx) { return std::make_tuple(pIdx, thetaIdx); });
80
81 // Fill a std::unordered_map with the linearised (p, theta) bin index as key and the vector index
82 // (i.e. the filtered RDataFrame row idx) as value for fast lookup.
83 unsigned int iRow(0);
84 for (const auto& tup : pAndThetaIdxs) {
85 double linBinIdx = (std::get<0>(tup) - 1.0) + (std::get<1>(tup) - 1.0) * weightsTable.m_nPBins;
86 weightsTable.m_linBinIdxsToRowIdxs.insert(std::make_pair(linBinIdx, iRow));
87 ++iRow;
88 }
89
90 // Store the vector (column) of weights per detector.
91 for (const Const::EDetector& det : Const::PIDDetectorSet::set()) {
92 auto detName = Const::parseDetectors(det);
93 auto colName = "ablat_s_" + detName;
94 auto weights = filteredRDF.Take<double>(colName.c_str()).GetValue();
95 weightsTable.m_weightsPerDet.insert(std::make_pair(detName, weights));
96 }
97
98 } else {
99
100 B2WARNING("Couldn't find detector weights in input ROOT::RDataFrame for std charged particle hypothesis: " << hypo.getPDGCode());
101 weightsTable.m_isEmpty = true;
102
103 }
104
105 m_weightsTablePerHypo.insert(std::make_pair(hypo.getPDGCode(), weightsTable));
106
107 };
108
109}
110
111
112double PIDDetectorWeights::getWeight(Const::ChargedStable hypo, Const::EDetector det, double p, double theta) const
113{
114
115 const auto weightsTable = &m_weightsTablePerHypo.at(hypo.getPDGCode());
116
117 if (weightsTable->m_isEmpty) {
118 // Weights not available for this pdgId.
119 return std::numeric_limits<float>::quiet_NaN();
120 }
121
122 // Get the bin index of the input p value.
123 auto pBinUpperEdge_it = std::lower_bound(std::begin(weightsTable->m_pBinEdges),
124 std::end(weightsTable->m_pBinEdges),
125 p);
126 auto pBinIdx = std::distance(std::begin(weightsTable->m_pBinEdges), pBinUpperEdge_it);
127 // Get the bin index of the input theta value.
128 auto thetaBinUpperEdge_it = std::lower_bound(std::begin(weightsTable->m_thetaBinEdges),
129 std::end(weightsTable->m_thetaBinEdges),
130 theta);
131 auto thetaBinIdx = std::distance(std::begin(weightsTable->m_thetaBinEdges), thetaBinUpperEdge_it);
132
133 // Get the linearised bin index to look up.
134 double linBinIdx = (pBinIdx - 1.0) + (thetaBinIdx - 1.0) * weightsTable->m_nPBins;
135
136 if (!weightsTable->m_linBinIdxsToRowIdxs.count(linBinIdx)) {
137 // Out-of-bin range p or theta
138 B2DEBUG(11, "\n"
139 << "p = " << p << " [GeV/c], theta = " << theta << " [rad].\n"
140 << "Bin indexes: (" << pBinIdx << ", " << thetaBinIdx << ").\n"
141 << "Either input value is outside of bin range for PID detector weights:\n"
142 << "p : [" << *weightsTable->m_pBinEdges.begin() << ", " << *weightsTable->m_pBinEdges.rbegin() << "],\n"
143 << "theta: [" << *weightsTable->m_thetaBinEdges.begin() << ", " << *weightsTable->m_thetaBinEdges.rbegin() << "]");
144 return std::numeric_limits<float>::quiet_NaN();
145 }
146
147 auto rowIdx = weightsTable->m_linBinIdxsToRowIdxs.at(linBinIdx);
148
149 return weightsTable->m_weightsPerDet.at(Const::parseDetectors(det)).at(rowIdx);
150
151};
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:589
int getPDGCode() const
PDG code.
Definition: Const.h:473
static DetectorSet set()
Accessor for the set of valid detector IDs.
Definition: Const.h:333
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:618
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:42
static std::string parseDetectors(EDetector det)
Converts Const::EDetector object to string.
Definition: UnitConst.cc:162
Nested class acting as a container the per-detector weights.
std::unordered_map< double, unsigned int > m_linBinIdxsToRowIdxs
Map the linearised index of the (p, theta) bin indexes to the row index in the (filtered) table.
std::map< std::string, std::vector< double > > m_weightsPerDet
Map each detector to its vector of weights.
std::set< double > m_pBinEdges
Set of p bins edges.
std::set< double > m_thetaBinEdges
Set of theta bins edges.
void fillWeightsTablePerHypoFromRDF()
Fill the internal weights container class per particle hypo, based on the content of the RDataFrame.
ROOT::RDataFrame m_weightsRDataFrame
The RDataFrame containing the detector weights per particle hypo, per phase space bin.
double getWeight(Const::ChargedStable hypo, Const::EDetector det, double p, double theta) const
Lookup the weight from the internal map structures.
std::map< int, WeightsTable > m_weightsTablePerHypo
Map containing a WeightsTable object per particle hypo.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24