Belle II Software light-2406-ragdoll
PIDLikelihood.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 <mdst/dataobjects/PIDLikelihood.h>
10#include <framework/logging/Logger.h>
11
12#include <TROOT.h>
13#include <TColor.h>
14#include <TParticlePDG.h>
15
16#include <cmath>
17#include <iostream>
18#include <iomanip>
19#include <algorithm>
20#include <vector>
21
22using namespace std;
23using namespace Belle2;
24
26{
27 for (unsigned short i = 0; i < Const::PIDDetectors::c_size; i++) {
28 for (unsigned int k = 0; k < Const::ChargedStable::c_SetSize; k++) {
29 m_logl[i][k] = 0.0;
30 }
31 }
32}
33
34
36 const Const::ChargedStable& part,
37 float logl)
38{
39 int index = Const::PIDDetectors::set().getIndex(det);
40 if (index < 0) {
41 B2ERROR("PIDLikelihood::setLogLikelihood: detector is not a PID device");
42 return;
43 }
44 if (logl != logl or logl == INFINITY) {
45 B2ERROR("PIDLikelihood::setLogLikelihood: log-likelihood for detector " << det << " is " << logl <<
46 " (i.e. +inf or NaN)! Ignoring this value. ("
47 << Const::SVD << "=SVD, " << Const::CDC << "=CDC, " << Const::TOP << "=TOP, "
48 << Const::ARICH << "=ARICH, " << Const::ECL << "=ECL, " << Const::KLM << "=KLM)");
49
50 return;
51 }
52 m_detectors += det;
53 m_logl[index][part.getIndex()] = logl;
54}
55
56
58{
59 for (auto& detectorLLs : m_logl) {
60 auto maxLL = detectorLLs[0];
61 for (auto LL : detectorLLs) maxLL = std::max(maxLL, LL);
62 for (auto& LL : detectorLLs) LL -= maxLL;
63 }
64}
65
66
68{
69 for (const auto& detector : set) {
70 if (m_detectors.contains(detector)) return true;
71 }
72 return false;
73}
74
75
77 Const::PIDDetectorSet set) const
78{
79 float result = 0;
81 it != Const::PIDDetectorSet::set().end(); ++it) {
82 if (set.contains(it))
83 result += m_logl[it.getIndex()][part.getIndex()];
84 }
85 return result;
86}
87
88
90 const Const::ChargedStable& p2,
91 double ratio,
92 Const::PIDDetectorSet set) const
93{
94 if (ratio < 0) {
95 B2ERROR("PIDLikelihood::probability argument 'ratio' is given with negative value");
96 return 0;
97 }
98 if (ratio == 0) return 0;
99
100 double dlogl = getLogL(p2, set) - getLogL(p1, set);
101 double result = 0;
102 if (dlogl < 0) {
103 double elogl = exp(dlogl);
104 result = ratio / (ratio + elogl);
105 } else {
106 double elogl = exp(-dlogl) * ratio; // to prevent overflow for very large dlogl
107 result = elogl / (1.0 + elogl);
108 }
109
110 return result;
111}
112
113
115 const double* fractions,
116 Const::PIDDetectorSet detSet) const
117{
118 return 1 / (1 + exp(-getLogarithmicProbability(part, fractions, detSet)));
119}
120
121
123 const double* fractions,
124 Const::PIDDetectorSet detSet) const
125{
126 // Defined as log(p / (1 - p)) where p is global PID probability.
127 // In order to save the range of return values it must be implemented as below.
128
129 int k = part.getIndex();
130 if (k < 0) return -INFINITY; // p = 0
131
132 if (fractions and fractions[k] <= 0) return -INFINITY; // p = 0
133
134 double logL_part = getLogL(part, detSet);
135 std::vector<double> logLs; // all other log likelihoods with priors > 0
136 std::vector<double> priorRatios;
137 for (int i = 0; i < static_cast<int>(Const::ChargedStable::c_SetSize); i++) {
138 if (i == k) continue;
139 double ratio = fractions ? fractions[i] / fractions[k] : 1;
140 if (ratio <= 0) continue;
141 priorRatios.push_back(ratio);
142 logLs.push_back(getLogL(Const::chargedStableSet.at(i), detSet));
143 }
144 if (logLs.empty()) return +INFINITY; // p = 1
145
146 double logL_max = logLs[0]; // maximal log likelihood of the others
147 for (auto logl : logLs) logL_max = std::max(logl, logL_max);
148
149 double sum = 0;
150 for (size_t i = 0; i < logLs.size(); i++) sum += priorRatios[i] * exp(logLs[i] - logL_max);
151 if (sum == 0) return +INFINITY; // to suppress cppckeck warning although this can never happen (see prev. return)
152
153 return logL_part - logL_max - log(sum);
154}
155
156
158 Const::PIDDetectorSet detSet) const
159{
160 const unsigned int n = Const::ChargedStable::c_SetSize;
161 double prob[n];
162 probability(prob, fractions, detSet);
163
164 int k = 0;
165 double maxProb = prob[k];
166 for (unsigned i = 0; i < n; ++i) {
167 if (prob[i] > maxProb) {maxProb = prob[i]; k = i;}
168 }
169 return Const::chargedStableSet.at(k);
170
171}
172
173
174void PIDLikelihood::probability(double probabilities[],
175 const double* fractions,
176 Const::PIDDetectorSet detSet) const
177{
178 const unsigned int n = Const::ChargedStable::c_SetSize;
179 double frac[n];
180 if (!fractions) {
181 for (unsigned int i = 0; i < n; ++i) frac[i] = 1.0; // normalization not needed
182 fractions = frac;
183 }
184
185 double logL[n];
186 double logLmax = 0;
187 bool hasMax = false;
188 for (unsigned i = 0; i < n; ++i) {
189 logL[i] = 0;
190 if (fractions[i] > 0) {
191 logL[i] = getLogL(Const::chargedStableSet.at(i), detSet);
192 if (!hasMax || logL[i] > logLmax) {
193 logLmax = logL[i];
194 hasMax = true;
195 }
196 }
197 }
198
199 double norm = 0;
200 for (unsigned i = 0; i < n; ++i) {
201 probabilities[i] = 0;
202 if (fractions[i] > 0) probabilities[i] = exp(logL[i] - logLmax) * fractions[i];
203 norm += probabilities[i];
204 }
205 if (norm == 0) return;
206
207 for (unsigned i = 0; i < n; ++i) {
208 probabilities[i] /= norm;
209 }
210
211}
212
213
215{
216
217 const string detectorName[Const::PIDDetectors::c_size] =
218 {"SVD", "CDC", "TOP", "ARICH", "ECL", "KLM"};
219 string hline("-------");
220 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
221 hline += "-----------";
222 string Hline("=======");
223 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
224 Hline += "===========";
225
226 cout << Hline << endl;
227
228 cout << "PDGcode";
229 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
230 cout << setw(10) << Const::chargedStableSet.at(i).getPDGCode() << " ";
231 cout << endl;
232
233 cout << Hline << endl;
234
235 float sum_logl[Const::ChargedStable::c_SetSize];
236 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) sum_logl[i] = 0;
237
238 for (unsigned k = 0; k < Const::PIDDetectors::c_size; k++) {
239 cout << setw(7) << detectorName[k];
240 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
241 cout << setw(10) << setprecision(4) << m_logl[k][i] << " ";
242 sum_logl[i] += m_logl[k][i];
243 }
244 cout << endl;
245 }
246
247 cout << hline << endl;
248
249 cout << setw(7) << "sum";
250 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
251 cout << setw(10) << setprecision(4) << sum_logl[i] << " ";
252 cout << endl;
253
254 if (isAvailable(Const::SVD) or isAvailable(Const::CDC) or isAvailable(Const::TOP) or
255 isAvailable(Const::ARICH) or isAvailable(Const::ECL) or isAvailable(Const::KLM)) {
256 unsigned k = 0;
257 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
258 if (sum_logl[i] > sum_logl[k]) k = i;
259 unsigned pos = 11 * (k + 1) + 3;
260 Hline.replace(pos, 1, "^");
261 }
262 cout << Hline << endl;
263
264}
265
266std::string PIDLikelihood::getInfoHTML() const
267{
268 const string detectorName[Const::PIDDetectors::c_size] = {"SVD", "CDC", "TOP", "ARICH", "ECL", "KLM"};
269
270 std::stringstream stream;
271 stream << std::setprecision(4);
272
273 //colors from event display
274 std::string colour[Const::ChargedStable::c_SetSize];
275 colour[0] = gROOT->GetColor(kAzure)->AsHexString();
276 colour[1] = gROOT->GetColor(kCyan + 1)->AsHexString();
277 colour[2] = gROOT->GetColor(kGray + 1)->AsHexString();
278 colour[3] = gROOT->GetColor(kRed + 1)->AsHexString();
279 colour[4] = gROOT->GetColor(kOrange - 2)->AsHexString();
280 colour[5] = gROOT->GetColor(kMagenta)->AsHexString();
281
282 stream << "<b>Probabilities</b><br>";
283 stream << "<table border=0 width=100%>";
284 stream << "<tr><td>All</td><td>";
285 stream << "<table cellspacing=1 border=0 width=100%><tr>";
286 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
287 double p = getProbability(Const::chargedStableSet.at(i));
288 stream << "<td bgcolor=\"" << colour[i] << "\" width=" << p * 100 << "%></td>";
289 }
290 stream << "</tr></table></td></tr>";
291
293 it != Const::PIDDetectorSet::set().end(); ++it) {
294 auto det = *it;
295 stream << "<tr><td>" << detectorName[it.getIndex()] << "</td><td>";
296 if (!isAvailable(det)) {
297 stream << "</td></tr>";
298 continue;
299 }
300 stream << "<table cellspacing=1 border=0 width=100%><tr>";
301 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
302 double p = getProbability(Const::chargedStableSet.at(i), 0, det);
303 stream << "<td bgcolor=\"" << colour[i] << "\" width=" << p * 100 << "%></td>";
304 }
305 stream << "</tr></table></td></tr>";
306 }
307 stream << "</table>\n";
308
309 stream << "<b>Log-Likelihoods</b><br>";
310 stream << "<table>";
311 stream << "<tr><th>PID / Detector</th>";
312 for (unsigned k = 0; k < Const::PIDDetectors::c_size; k++)
313 stream << "<th>" << detectorName[k] << "</th>";
314 stream << "</tr>";
315 for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
316 stream << "<tr>";
317 stream << "<td bgcolor=\"" << colour[i] << "\">" << Const::chargedStableSet.at(i).getParticlePDG()->GetName() << "</td>";
318 for (unsigned k = 0; k < Const::PIDDetectors::c_size; k++) {
319 stream << "<td>" << m_logl[k][i] << "</td>";
320 }
321 stream << "</tr>";
322 }
323 stream << "</table>";
324 stream << "<br>";
325
326 return stream.str();
327}
328
329void PIDLikelihood::addPreOfficialLikelihood(const std::string& preOfficialIdentifier, const double preOfficialLikelihood)
330{
331 m_preOfficialLikelihoods[preOfficialIdentifier] = preOfficialLikelihood;
332}
333
334double PIDLikelihood::getPreOfficialLikelihood(const std::string& preOfficialIdentifier) const
335{
336 if (m_preOfficialLikelihoods.count(preOfficialIdentifier) == 0) {
337 B2WARNING("PIDLikelihood::getPreOfficialLikelihood: preOfficialIdentifier " << preOfficialIdentifier << " does not exist. ");
338 return -1.0;
339 }
340 return m_preOfficialLikelihoods.at(preOfficialIdentifier);
341}
342
343
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:589
static const unsigned int c_SetSize
Number of elements (for use in array bounds etc.)
Definition: Const.h:615
int getIndex(EDetector det) const
Getter for the index of a given detector in this set.
Definition: UnitConst.cc:278
Iterator end() const
Ending iterator.
Definition: UnitConst.cc:220
bool contains(const DetectorSet &set) const
Check whether this set contains another set.
Definition: Const.h:235
static DetectorSet set()
Accessor function for the set of valid detectors.
Definition: Const.h:374
static const size_t c_size
Number of PID detectors, temporary workaround.
Definition: Const.h:376
const ParticleType & at(unsigned int index) const
Return particle at given index, or end() if out of range.
Definition: Const.h:549
int getPDGCode() const
PDG code.
Definition: Const.h:473
int getIndex() const
This particle's index in the associated set.
Definition: Const.h:461
const TParticlePDG * getParticlePDG() const
Accessor for ROOT TParticlePDG object.
Definition: UnitConst.cc:351
A class for sets of detector IDs whose content is limited to restricted set of valid detector IDs.
Definition: Const.h:296
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
PIDLikelihood()
Default constructor: log likelihoods and flags set to 0.
Const::ChargedStable getMostLikely(const double *fractions=nullptr, Const::PIDDetectorSet set=Const::PIDDetectorSet::set()) const
Return most likely particle among chargedStableSet; if prior fractions not given equal prior probabil...
void subtractMaximum()
Subtract the maximum of log likelihoods of each detector component in order to reduce the range of va...
std::map< std::string, double > m_preOfficialLikelihoods
Internal storage of pre-official likelihood.
double getLogarithmicProbability(const Const::ChargedStable &part, const double *fractions=nullptr, Const::PIDDetectorSet set=Const::PIDDetectorSet::set()) const
Return logarithmic equivalent of likelihood probability defined as log(p/(1-p)), where p is the combi...
bool isAvailable(Const::PIDDetectorSet set=Const::PIDDetectorSet::set()) const
Check whether PID information is available for at least one of the detectors in a given set.
std::string getInfoHTML() const override
Return HTML Info of PID Likelihoods.
float m_logl[Const::PIDDetectors::c_size][Const::ChargedStable::c_SetSize]
log likelihoods
void printArray() const
Prints the content of a private array of log likelihoods.
float getLogL(const Const::ChargedStable &part, Const::PIDDetectorSet set=Const::PIDDetectorSet::set()) const
Return log likelihood for a given detector set and particle.
void setLogLikelihood(Const::EDetector det, const Const::ChargedStable &part, float logl)
Set log likelihood for a given detector and particle.
void addPreOfficialLikelihood(const std::string &preOfficialIdentifier, const double preOfficialLikelihood)
Add the pre-official likelihood.
void probability(double probabilities[], const double *fractions, Const::PIDDetectorSet detSet) const
Calculate likelihood probabilities.
double getProbability(const Const::ChargedStable &p1, const Const::ChargedStable &p2, Const::PIDDetectorSet set=Const::PIDDetectorSet::set()) const
Return combined likelihood probability for a particle being p1 and not p2, assuming equal prior proba...
Const::DetectorSet m_detectors
set of detectors with PID information
double getPreOfficialLikelihood(const std::string &preOfficialIdentifier) const
Get the pre-official likelihood.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.