Belle II Software  release-05-02-19
Cluster.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Torben Ferber *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <mdst/dataobjects/Cluster.h>
12 #include <framework/logging/Logger.h>
13 
14 #include <cmath>
15 #include <iostream>
16 #include <iomanip>
17 
18 using namespace std;
19 using namespace Belle2;
20 
21 Cluster::Cluster()
22 {
23  for (unsigned short i = 0; i < Const::ClusterDetectors::c_size; i++) {
24  for (unsigned int k = 0; k < Const::Cluster::c_SetSize; k++) {
25  m_logl[i][k] = 0.0;
26  }
27  }
28 }
29 
30 
31 void Cluster::setLogLikelihood(Const::EDetector det,
32  const Const::Cluster& part,
33  float logl)
34 {
35  int index = Const::ClusterDetectors::set().getIndex(det);
36  if (index < 0) {
37  B2ERROR("ClusterLikelihood::setLogLikelihood: detector is not a cluster device");
38  return;
39  }
40  if (logl != logl or logl == INFINITY) {
41  B2ERROR("ClusterLikelihood::setLogLikelihood: log-likelihood for detector " << det << " is " << logl <<
42  " (i.e. +inf or NaN)! Ignoring this value. (" << Const::ECL << "=ECL, " << Const::KLM << "=KLM");
43 
44  return;
45  }
46  m_detectors += det;
47  m_logl[index][part.getIndex()] = logl;
48 }
49 
50 
51 float Cluster::getLogL(const Const::Cluster& cluster,
52  Const::ClusterDetectorSet set) const
53 {
54  float result = 0;
55  for (unsigned int index = 0; index < Const::ClusterDetectorSet::set().size(); ++index) {
56  if (set.contains(Const::ClusterDetectorSet::set()[index]))
57  result += m_logl[index][cluster.getIndex()];
58  }
59  return result;
60 }
61 
62 
63 Double32_t Cluster::getProbability(const Const::Cluster& c1,
64  const Const::Cluster& c2,
65  double ratio,
66  Const::ClusterDetectorSet set) const
67 {
68  if (ratio < 0) {
69  B2ERROR("Cluster::probability argument 'ratio' is given with negative value");
70  return 0;
71  }
72  if (ratio == 0) return 0;
73 
74  double dlogl = getLogL(c2, set) - getLogL(c1, set);
75  double res;
76  if (dlogl < 0) {
77  double elogl = exp(dlogl);
78  res = ratio / (ratio + elogl);
79  } else {
80  double elogl = exp(-dlogl) * ratio; // to prevent overflow for very large dlogl
81  res = elogl / (1.0 + elogl);
82  }
83  //TODO: only necessary if one wants to use mcprod1405 MC sample. Remove when there's a good replacement.
84  if (std::isfinite(res))
85  return res;
86  return 0.0;
87 }
88 
89 double Cluster::getProbability(const Const::Cluster& cluster,
90  const double* fractions,
91  Const::ClusterDetectorSet detSet) const
92 {
93  double prob[Const::Cluster::c_SetSize];
94  probability(prob, fractions, detSet);
95 
96  int k = cluster.getIndex();
97  if (k < 0) return 0;
98 
99  return prob[k];
100 
101 }
102 
103 Const::Cluster Cluster::getMostLikely(const double* fractions,
104  Const::ClusterDetectorSet detSet) const
105 {
106  const unsigned int n = Const::Cluster::c_SetSize;
107  double prob[n];
108  probability(prob, fractions, detSet);
109 
110  int k = 0;
111  double maxProb = prob[k];
112  for (unsigned i = 0; i < n; ++i) {
113  if (prob[i] > maxProb) {maxProb = prob[i]; k = i;}
114  }
115  return Const::clusterSet.at(k);
116 
117 }
118 
119 
120 void Cluster::probability(double probabilities[],
121  const double* fractions,
122  Const::ClusterDetectorSet detSet) const
123 {
124  const unsigned int n = Const::Cluster::c_SetSize;
125  double frac[n];
126  if (!fractions) {
127  for (unsigned int i = 0; i < n; ++i) frac[i] = 1.0; // normalization not needed
128  fractions = frac;
129  }
130 
131  double logL[n];
132  double logLmax = 0;
133  bool hasMax = false;
134  for (unsigned i = 0; i < n; ++i) {
135  logL[i] = 0;
136  if (fractions[i] > 0) {
137  logL[i] = getLogL(Const::clusterSet.at(i), detSet);
138  if (!hasMax || logL[i] > logLmax) {
139  logLmax = logL[i];
140  hasMax = true;
141  }
142  }
143  }
144 
145  double norm = 0;
146  for (unsigned i = 0; i < n; ++i) {
147  probabilities[i] = 0;
148  if (fractions[i] > 0) probabilities[i] = exp(logL[i] - logLmax) * fractions[i];
149  norm += probabilities[i];
150  }
151  if (norm == 0) return;
152 
153  for (unsigned i = 0; i < n; ++i) {
154  probabilities[i] /= norm;
155  }
156 
157 }
158 
159 
160 void Cluster::printArray() const
161 {
162 
163  string detectorName[Const::ClusterDetectors::c_size] = {"ECL", "KLM"};
164  string hline("-------");
165  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++)
166  hline += "-----------";
167  string Hline("=======");
168  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++)
169  Hline += "===========";
170 
171  cout << Hline << endl;
172 
173  cout << "PDGcode";
174  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++)
175  cout << setw(10) << Const::clusterSet.at(i).getPDGCode() << " ";
176  cout << endl;
177 
178  cout << Hline << endl;
179 
180  float sum_logl[Const::Cluster::c_SetSize];
181  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++) sum_logl[i] = 0;
182 
183  for (unsigned k = 0; k < Const::ClusterDetectors::c_size; k++) {
184  cout << setw(7) << detectorName[k];
185  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++) {
186  cout << setw(10) << setprecision(4) << m_logl[k][i] << " ";
187  sum_logl[i] += m_logl[k][i];
188  }
189  cout << endl;
190  }
191 
192  cout << hline << endl;
193 
194  cout << setw(7) << "sum";
195  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++)
196  cout << setw(10) << setprecision(4) << sum_logl[i] << " ";
197  cout << endl;
198 
199  if (isAvailable(Const::ECL) or isAvailable(Const::KLM)) {
200  unsigned k = 0;
201  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++)
202  if (sum_logl[i] > sum_logl[k]) k = i;
203  unsigned pos = 11 * (k + 1) + 3;
204  Hline.replace(pos, 1, "^");
205  }
206  cout << Hline << endl;
207 
208 }
209 
210 std::string Cluster::getInfoHTML() const
211 {
212  string detectorName[Const::ClusterDetectors::c_size] = {"ECL", "KLM"};
213 
214  std::stringstream stream;
215  stream << std::setprecision(4);
216  stream << "<b>Likelihoods</b><br>";
217  stream << "<table>";
218  stream << "<tr><th>Cluster / Detector</th>";
219  for (unsigned k = 0; k < Const::ClusterDetectors::c_size; k++)
220  stream << "<th>" << detectorName[k] << "</th>";
221  stream << "</tr>";
222  for (unsigned i = 0; i < Const::Cluster::c_SetSize; i++) {
223  stream << "<tr>";
224  stream << "<td>" << Const::clusterSet.at(i).getPDGCode() << "</td>";
225  for (unsigned k = 0; k < Const::ClusterDetectors::c_size; k++) {
226  stream << "<td>" << m_logl[k][i] << "</td>";
227  }
228  stream << "</tr>";
229  }
230  stream << "</table>";
231  stream << "<br>";
232 
233  return stream.str();
234 }
235 
236 
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::Const::Cluster
Provides a type-safe way to pass members of the clusterSet set.
Definition: Const.h:502
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Const::DetectorSet::contains
bool contains(const DetectorSet &set) const
Check whether this set contains another set.
Definition: Const.h:116