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