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