9#include <mdst/dataobjects/PIDLikelihood.h>
10#include <framework/logging/Logger.h>
14#include <TParticlePDG.h>
41 B2ERROR(
"PIDLikelihood::setLogLikelihood: detector is not a PID device");
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)");
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;
69 for (
const auto& detector : set) {
95 B2ERROR(
"PIDLikelihood::probability argument 'ratio' is given with negative value");
98 if (ratio == 0)
return 0;
103 double elogl = exp(dlogl);
104 result = ratio / (ratio + elogl);
106 double elogl = exp(-dlogl) * ratio;
107 result = elogl / (1.0 + elogl);
115 const double* fractions,
123 const double* fractions,
130 if (k < 0)
return -INFINITY;
132 if (fractions and fractions[k] <= 0)
return -INFINITY;
134 double logL_part =
getLogL(part, detSet);
135 std::vector<double> logLs;
136 std::vector<double> priorRatios;
138 if (i == k)
continue;
139 double ratio = fractions ? fractions[i] / fractions[k] : 1;
140 if (ratio <= 0)
continue;
141 priorRatios.push_back(ratio);
144 if (logLs.empty())
return +INFINITY;
146 double logL_max = logLs[0];
147 for (
auto logl : logLs) logL_max = std::max(logl, logL_max);
150 for (
size_t i = 0; i < logLs.size(); i++) sum += priorRatios[i] * exp(logLs[i] - logL_max);
151 if (sum == 0)
return +INFINITY;
153 return logL_part - logL_max - log(sum);
165 double maxProb = prob[k];
166 for (
unsigned i = 0; i < n; ++i) {
167 if (prob[i] > maxProb) {maxProb = prob[i]; k = i;}
175 const double* fractions,
181 for (
unsigned int i = 0; i < n; ++i) frac[i] = 1.0;
188 for (
unsigned i = 0; i < n; ++i) {
190 if (fractions[i] > 0) {
192 if (!hasMax || logL[i] > logLmax) {
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];
205 if (norm == 0)
return;
207 for (
unsigned i = 0; i < n; ++i) {
208 probabilities[i] /= norm;
218 {
"SVD",
"CDC",
"TOP",
"ARICH",
"ECL",
"KLM"};
219 string hline(
"-------");
221 hline +=
"-----------";
222 string Hline(
"=======");
224 Hline +=
"===========";
226 cout << Hline << endl;
233 cout << Hline << endl;
239 cout << setw(7) << detectorName[k];
241 cout << setw(10) << setprecision(4) <<
m_logl[k][i] <<
" ";
242 sum_logl[i] +=
m_logl[k][i];
247 cout << hline << endl;
249 cout << setw(7) <<
"sum";
251 cout << setw(10) << setprecision(4) << sum_logl[i] <<
" ";
258 if (sum_logl[i] > sum_logl[k]) k = i;
259 unsigned pos = 11 * (k + 1) + 3;
260 Hline.replace(pos, 1,
"^");
262 cout << Hline << endl;
270 std::stringstream stream;
271 stream << std::setprecision(4);
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();
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>";
288 stream <<
"<td bgcolor=\"" << colour[i] <<
"\" width=" << p * 100 <<
"%></td>";
290 stream <<
"</tr></table></td></tr>";
295 stream <<
"<tr><td>" << detectorName[it.getIndex()] <<
"</td><td>";
297 stream <<
"</td></tr>";
300 stream <<
"<table cellspacing=1 border=0 width=100%><tr>";
303 stream <<
"<td bgcolor=\"" << colour[i] <<
"\" width=" << p * 100 <<
"%></td>";
305 stream <<
"</tr></table></td></tr>";
307 stream <<
"</table>\n";
309 stream <<
"<b>Log-Likelihoods</b><br>";
311 stream <<
"<tr><th>PID / Detector</th>";
313 stream <<
"<th>" << detectorName[k] <<
"</th>";
317 stream <<
"<td bgcolor=\"" << colour[i] <<
"\">" <<
Const::chargedStableSet.at(i).getParticlePDG()->GetName() <<
"</td>";
319 stream <<
"<td>" <<
m_logl[k][i] <<
"</td>";
323 stream <<
"</table>";
337 B2WARNING(
"PIDLikelihood::getPreOfficialLikelihood: preOfficialIdentifier " << preOfficialIdentifier <<
" does not exist. ");
Provides a type-safe way to pass members of the chargedStableSet set.
static const unsigned int c_SetSize
Number of elements (for use in array bounds etc.)
int getIndex(EDetector det) const
Getter for the index of a given detector in this set.
Iterator end() const
Ending iterator.
static DetectorSet set()
Accessor function for the set of valid detectors.
static const size_t c_size
Number of PID detectors, temporary workaround.
int getIndex() const
This particle's index in the associated set.
RestrictedDetectorSet< PIDDetectors > PIDDetectorSet
Typedef for set of PID detectors.
static const ParticleSet chargedStableSet
set of charged stable particles
EDetector
Enum for identifying the detector components (detector and subdetector).
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.