11 #include <mdst/dataobjects/Cluster.h>
12 #include <framework/logging/Logger.h>
23 for (
unsigned short i = 0; i < Const::ClusterDetectors::c_size; i++) {
24 for (
unsigned int k = 0; k < Const::Cluster::c_SetSize; k++) {
35 int index = Const::ClusterDetectors::set().getIndex(det);
37 B2ERROR(
"ClusterLikelihood::setLogLikelihood: detector is not a cluster device");
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");
47 m_logl[index][part.getIndex()] = logl;
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()];
69 B2ERROR(
"Cluster::probability argument 'ratio' is given with negative value");
72 if (ratio == 0)
return 0;
74 double dlogl = getLogL(c2, set) - getLogL(c1, set);
77 double elogl = exp(dlogl);
78 res = ratio / (ratio + elogl);
80 double elogl = exp(-dlogl) * ratio;
81 res = elogl / (1.0 + elogl);
84 if (std::isfinite(res))
90 const double* fractions,
93 double prob[Const::Cluster::c_SetSize];
94 probability(prob, fractions, detSet);
96 int k = cluster.getIndex();
106 const unsigned int n = Const::Cluster::c_SetSize;
108 probability(prob, fractions, detSet);
111 double maxProb = prob[k];
112 for (
unsigned i = 0; i < n; ++i) {
113 if (prob[i] > maxProb) {maxProb = prob[i]; k = i;}
115 return Const::clusterSet.at(k);
120 void Cluster::probability(
double probabilities[],
121 const double* fractions,
124 const unsigned int n = Const::Cluster::c_SetSize;
127 for (
unsigned int i = 0; i < n; ++i) frac[i] = 1.0;
134 for (
unsigned i = 0; i < n; ++i) {
136 if (fractions[i] > 0) {
137 logL[i] = getLogL(Const::clusterSet.at(i), detSet);
138 if (!hasMax || logL[i] > logLmax) {
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];
151 if (norm == 0)
return;
153 for (
unsigned i = 0; i < n; ++i) {
154 probabilities[i] /= norm;
160 void Cluster::printArray()
const
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 +=
"===========";
171 cout << Hline << endl;
174 for (
unsigned i = 0; i < Const::Cluster::c_SetSize; i++)
175 cout << setw(10) << Const::clusterSet.at(i).getPDGCode() <<
" ";
178 cout << Hline << endl;
180 float sum_logl[Const::Cluster::c_SetSize];
181 for (
unsigned i = 0; i < Const::Cluster::c_SetSize; i++) sum_logl[i] = 0;
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];
192 cout << hline << endl;
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] <<
" ";
199 if (isAvailable(Const::ECL) or isAvailable(Const::KLM)) {
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,
"^");
206 cout << Hline << endl;
210 std::string Cluster::getInfoHTML()
const
212 string detectorName[Const::ClusterDetectors::c_size] = {
"ECL",
"KLM"};
214 std::stringstream stream;
215 stream << std::setprecision(4);
216 stream <<
"<b>Likelihoods</b><br>";
218 stream <<
"<tr><th>Cluster / Detector</th>";
219 for (
unsigned k = 0; k < Const::ClusterDetectors::c_size; k++)
220 stream <<
"<th>" << detectorName[k] <<
"</th>";
222 for (
unsigned i = 0; i < Const::Cluster::c_SetSize; i++) {
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>";
230 stream <<
"</table>";