Belle II Software  release-05-01-25
PIDLikelihood.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Marko Staric, Thomas Kuhr *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <mdst/dataobjects/PIDLikelihood.h>
12 #include <framework/logging/Logger.h>
13 
14 #include <TROOT.h>
15 #include <TColor.h>
16 #include <TParticlePDG.h>
17 
18 #include <cmath>
19 #include <iostream>
20 #include <iomanip>
21 
22 using namespace std;
23 using namespace Belle2;
24 
25 PIDLikelihood::PIDLikelihood()
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 
35 void PIDLikelihood::setLogLikelihood(Const::EDetector det,
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. (" << Const::CDC << "=CDC, " << Const::TOP << "=TOP, " << Const::ARICH << "=ARICH, " <<
47  Const::ECL << "=ECL)");
48 
49  return;
50  }
51  m_detectors += det;
52  m_logl[index][part.getIndex()] = logl;
53 }
54 
55 
56 float PIDLikelihood::getLogL(const Const::ChargedStable& part,
57  Const::PIDDetectorSet set) const
58 {
59  float result = 0;
60  for (unsigned int index = 0; index < Const::PIDDetectorSet::set().size(); ++index) {
61  if (set.contains(Const::PIDDetectorSet::set()[index]))
62  result += m_logl[index][part.getIndex()];
63  }
64  return result;
65 }
66 
67 
68 double PIDLikelihood::getProbability(const Const::ChargedStable& p1,
69  const Const::ChargedStable& p2,
70  double ratio,
71  Const::PIDDetectorSet set) const
72 {
73  if (ratio < 0) {
74  B2ERROR("PIDLikelihood::probability argument 'ratio' is given with negative value");
75  return 0;
76  }
77  if (ratio == 0) return 0;
78 
79  double dlogl = getLogL(p2, set) - getLogL(p1, set);
80  double res;
81  if (dlogl < 0) {
82  double elogl = exp(dlogl);
83  res = ratio / (ratio + elogl);
84  } else {
85  double elogl = exp(-dlogl) * ratio; // to prevent overflow for very large dlogl
86  res = elogl / (1.0 + elogl);
87  }
88  //TODO: only necessary if one wants to use mcprod1405 MC sample. Remove when there's a good replacement.
89  if (std::isfinite(res))
90  return res;
91  return 0;
92 }
93 
94 double PIDLikelihood::getProbability(const Const::ChargedStable& part,
95  const double* fractions,
96  Const::PIDDetectorSet detSet) const
97 {
98  double prob[Const::ChargedStable::c_SetSize];
99  probability(prob, fractions, detSet);
100 
101  int k = part.getIndex();
102  if (k < 0) return 0;
103 
104  return prob[k];
105 
106 }
107 
108 Const::ChargedStable PIDLikelihood::getMostLikely(const double* fractions,
109  Const::PIDDetectorSet detSet) const
110 {
111  const unsigned int n = Const::ChargedStable::c_SetSize;
112  double prob[n];
113  probability(prob, fractions, detSet);
114 
115  int k = 0;
116  double maxProb = prob[k];
117  for (unsigned i = 0; i < n; ++i) {
118  if (prob[i] > maxProb) {maxProb = prob[i]; k = i;}
119  }
120  return Const::chargedStableSet.at(k);
121 
122 }
123 
124 
125 void PIDLikelihood::probability(double probabilities[],
126  const double* fractions,
127  Const::PIDDetectorSet detSet) const
128 {
129  const unsigned int n = Const::ChargedStable::c_SetSize;
130  double frac[n];
131  if (!fractions) {
132  for (unsigned int i = 0; i < n; ++i) frac[i] = 1.0; // normalization not needed
133  fractions = frac;
134  }
135 
136  double logL[n];
137  double logLmax = 0;
138  bool hasMax = false;
139  for (unsigned i = 0; i < n; ++i) {
140  logL[i] = 0;
141  if (fractions[i] > 0) {
142  logL[i] = getLogL(Const::chargedStableSet.at(i), detSet);
143  if (!hasMax || logL[i] > logLmax) {
144  logLmax = logL[i];
145  hasMax = true;
146  }
147  }
148  }
149 
150  double norm = 0;
151  for (unsigned i = 0; i < n; ++i) {
152  probabilities[i] = 0;
153  if (fractions[i] > 0) probabilities[i] = exp(logL[i] - logLmax) * fractions[i];
154  norm += probabilities[i];
155  }
156  if (norm == 0) return;
157 
158  for (unsigned i = 0; i < n; ++i) {
159  probabilities[i] /= norm;
160  }
161 
162 }
163 
164 
165 void PIDLikelihood::printArray() const
166 {
167 
168  string detectorName[Const::PIDDetectors::c_size] =
169  {"SVD", "CDC", "TOP", "ARICH", "ECL", "KLM"};
170  string hline("-------");
171  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
172  hline += "-----------";
173  string Hline("=======");
174  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
175  Hline += "===========";
176 
177  cout << Hline << endl;
178 
179  cout << "PDGcode";
180  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
181  cout << setw(10) << Const::chargedStableSet.at(i).getPDGCode() << " ";
182  cout << endl;
183 
184  cout << Hline << endl;
185 
186  float sum_logl[Const::ChargedStable::c_SetSize];
187  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) sum_logl[i] = 0;
188 
189  for (unsigned k = 0; k < Const::PIDDetectors::c_size; k++) {
190  cout << setw(7) << detectorName[k];
191  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
192  cout << setw(10) << setprecision(4) << m_logl[k][i] << " ";
193  sum_logl[i] += m_logl[k][i];
194  }
195  cout << endl;
196  }
197 
198  cout << hline << endl;
199 
200  cout << setw(7) << "sum";
201  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
202  cout << setw(10) << setprecision(4) << sum_logl[i] << " ";
203  cout << endl;
204 
205  if (isAvailable(Const::SVD) or isAvailable(Const::CDC) or isAvailable(Const::TOP) or
206  isAvailable(Const::ARICH) or isAvailable(Const::ECL) or isAvailable(Const::KLM)) {
207  unsigned k = 0;
208  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++)
209  if (sum_logl[i] > sum_logl[k]) k = i;
210  unsigned pos = 11 * (k + 1) + 3;
211  Hline.replace(pos, 1, "^");
212  }
213  cout << Hline << endl;
214 
215 }
216 
217 std::string PIDLikelihood::getInfoHTML() const
218 {
219  string detectorName[Const::PIDDetectors::c_size] = {"SVD", "CDC", "TOP", "ARICH", "ECL", "KLM"};
220 
221  std::stringstream stream;
222  stream << std::setprecision(4);
223 
224  //colors from event display
225  std::string colour[Const::ChargedStable::c_SetSize];
226  colour[0] = gROOT->GetColor(kAzure)->AsHexString();
227  colour[1] = gROOT->GetColor(kCyan + 1)->AsHexString();
228  colour[2] = gROOT->GetColor(kGray + 1)->AsHexString();
229  colour[3] = gROOT->GetColor(kRed + 1)->AsHexString();
230  colour[4] = gROOT->GetColor(kOrange - 2)->AsHexString();
231  colour[5] = gROOT->GetColor(kMagenta)->AsHexString();
232 
233  stream << "<b>Probabilities</b><br>";
234  stream << "<table border=0 width=100%>";
235  stream << "<tr><td>All</td><td>";
236  stream << "<table cellspacing=1 border=0 width=100%><tr>";
237  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
238  double p = getProbability(Const::chargedStableSet.at(i));
239  stream << "<td bgcolor=\"" << colour[i] << "\" width=" << p * 100 << "%></td>";
240  }
241  stream << "</tr></table></td></tr>";
242 
243  for (unsigned k = 0; k < Const::PIDDetectors::c_size; k++) {
244  auto det = Const::PIDDetectorSet::set()[k];
245  stream << "<tr><td>" << detectorName[k] << "</td><td>";
246  if (!isAvailable(det)) {
247  stream << "</td></tr>";
248  continue;
249  }
250  stream << "<table cellspacing=1 border=0 width=100%><tr>";
251  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
252  double p = getProbability(Const::chargedStableSet.at(i), 0, det);
253  stream << "<td bgcolor=\"" << colour[i] << "\" width=" << p * 100 << "%></td>";
254  }
255  stream << "</tr></table></td></tr>";
256  }
257  stream << "</table>\n";
258 
259  stream << "<b>Log-Likelihoods</b><br>";
260  stream << "<table>";
261  stream << "<tr><th>PID / Detector</th>";
262  for (unsigned k = 0; k < Const::PIDDetectors::c_size; k++)
263  stream << "<th>" << detectorName[k] << "</th>";
264  stream << "</tr>";
265  for (unsigned i = 0; i < Const::ChargedStable::c_SetSize; i++) {
266  stream << "<tr>";
267  stream << "<td bgcolor=\"" << colour[i] << "\">" << Const::chargedStableSet.at(i).getParticlePDG()->GetName() << "</td>";
268  for (unsigned k = 0; k < Const::PIDDetectors::c_size; k++) {
269  stream << "<td>" << m_logl[k][i] << "</td>";
270  }
271  stream << "</tr>";
272  }
273  stream << "</table>";
274  stream << "<br>";
275 
276  return stream.str();
277 }
278 
279 
Belle2::Const::EDetector
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:44
Belle2::Const::RestrictedDetectorSet
A class for sets of detector IDs whose content is limited to restricted set of valid detector IDs.
Definition: Const.h:172
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::Const::DetectorSet::contains
bool contains(const DetectorSet &set) const
Check whether this set contains another set.
Definition: Const.h:116