Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JTestPoissonLogLikelihoodRatio.hh
Go to the documentation of this file.
1#ifndef __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIO__
2#define __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIO__
3
4#include <istream>
5#include <ostream>
6
7#include "JLang/JException.hh"
8
11
12#include "TH1.h"
13
14
15/**
16 * \author bjung
17 */
18namespace JCOMPAREHISTOGRAMS {
19
20 /**
21 * Implementation of the Poisson log-likelihood ratio test.
22 * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis,
23 * which is compared to an alternative hypothesis given by the second histogram.
24 * The uncertainty on the expectation is taken into account as an additional Poissonian constraint term in the likelihood of each bin (\f$\Lambda_i\f$):
25 *
26 * \f[
27 * \Lambda_i = \frac{(\gamma_i m_i)^{n_i}}{n_i!} e^{-\gamma_i m_i} \cdot \frac{(\gamma_i M_i)^{M_i}}{M_i!} e^{-\gamma_i M_i},
28 * \f]
29 *
30 * where:
31 *
32 * - \f$\gamma_i\f$ corresponds to a scaling parameter,
33 * - \f$m_i\f$ corresponds to the expectation in bin \f$i\f$,
34 * - \f$n_i\f$ corresponds to the observed number of events in bin \f$i\f$ of the Asimov dataset and
35 * - \f$M_i = \frac{m_i^2}{\sigma_{m_i}^2}\f$ corresponds to the effective number of exected events in bin \f$i\f$,\n where \f$\sigma_{m_i}\f$ corresponds to the standard deviation of the expectation \f$m_i\f$.
36 *
37 * This is inspired on the Beeston-Barlow method, described in the following publication:
38 * Computer Physics Communications, Volume 77, Issue 2, 1993
39 * "Fitting using finite Monte Carlo samples",
40 * Roger Barlow and Christine Beeston.
41 * DOI: 10.1016/0010-4655(93)90005-W
42 */
44 public JTest_t
45 {
46 public:
47
48 /**
49 * Default constructor.
50 */
52 JTest_t("Poisson_NLLR", "NLLR"),
53 threshold(0.0)
54 {}
55
56
57 /**
58 * Applies a Poissonian log-likelihood ratio test to the two given histograms.\n
59 * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis.\n
60 * The second histogram is treated as the expectation for the alternative hypothesis, to which the null hypothesis is compared.
61 *
62 * \param o1 First histogram
63 * \param o2 Second histogram
64 */
65 void test(const TObject* o1, const TObject* o2) override
66 {
67 using namespace std;
68 using namespace JPP;
69
70 const TH1* h1 = dynamic_cast<const TH1*>(o1);
71 const TH1* h2 = dynamic_cast<const TH1*>(o2);
72
73 if (h1 == NULL || h2 == NULL) {
74 THROW(JValueOutOfRange, "JTestPoissonLogLikelihood::test(): Could not cast given TObjects to TH1.");
75 }
76
77 if(h1->GetNbinsX() != h2->GetNbinsX() ||
78 h1->GetNbinsY() != h2->GetNbinsY() ||
79 h1->GetNbinsZ() != h2->GetNbinsZ()) {
80 THROW(JValueOutOfRange, "JTestPoissonLogLikelihood::test(): Histograms with different binning. The objects: " <<
81 h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
82 }
83
84 TH1* h3 = (TH1*) h1->Clone(h1->GetName() == h2->GetName() ?
85 MAKE_CSTRING(h1->GetName() << "_" << testName) :
86 MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_" << testName));
87
88 h3->Reset();
89
90 double LnLratio = 0.0;
91
92 for (int i=1 ; i <= h1->GetNbinsX() ; ++i) {
93 for (int j=1 ; j <= h1->GetNbinsY() ; ++j) {
94 for (int k=1 ; k <= h1->GetNbinsZ() ; ++k) {
95
96 double contribution = 0.0;
97
98 const double n1 = h1->GetBinContent(i,j,k);
99 const double n2 = h2->GetBinContent(i,j,k);
100
101 const double e1 = h1->GetBinError(i,j,k);
102 const double e2 = h2->GetBinError(i,j,k);
103
104 if (n2 > 0.0) {
105
106 if (!(e2 > 0.0)) {
107 THROW(JDivisionByZero, "JTestPoissonLogLikelihoodRatio::test(): Invalid effective bin content for histogram " << h2->GetName() << " at bin (" << i << ',' << j << ',' << k << ')');
108 }
109
110 const double N2 = n2*n2/e2/e2; // Effective bin content of alternative model
111
112 if (n1 > 0.0) {
113
114 if (!(e1 > 0.0)) {
115 THROW(JDivisionByZero, "JTestPoissonLogLikelihoodRatio::test(): Invalid effective bin content for histogram " << h1->GetName() << " at bin (" << i << ',' << j << ',' << k << ')');
116 }
117
118 const double f1 = (n1 + N2) / (n2 + N2); // Optimal scale factor for alternative model bin
119
120 const double N1 = n1*n1/e1/e1; // Effective bin content of Asimov dataset
121
122 contribution = 2 * (N2 - N1 +
123 lgamma(N2 + 1) - lgamma(N1 + 1) +
124 n1*log(n1/n2/f1) - N2*log(f1*N2) +
125 N1*log(N1));
126
127 } else {
128
129 const double f1 = N2 / (n2 + N2); // Optimal scale factor for alternative model bin
130
131 contribution = 2 * (N2 + lgamma(N2 + 1) - N2*log(f1*N2));
132 }
133 }
134
135 h3->SetBinContent(i,j,k, contribution);
136
137 LnLratio += contribution;
138 }
139 }
140 }
141
142 const bool passed = (LnLratio < threshold);
143
144 const JResultTitle title(testName, resultType, passed, LnLratio);
145
146 h3->SetTitle(title.getTitle().c_str());
147 h3->GetYaxis()->SetTitle(resultType.c_str());
148
149 const JTestResult r (testName,
150 JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
151 JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
152 resultType, LnLratio, threshold, h3, passed);
153
154 this->push_back(r);
155 }
156
157
158 /**
159 * Read test parameters from input.
160 *
161 * \param in input stream
162 * \return input stream
163 */
164 std::istream& read(std::istream& in) override
165 {
166 using namespace JPP;
167
168 in >> threshold;
169
170 if (!(threshold > 0.0)) {
171 THROW(JValueOutOfRange, "JTestPoissonLogLikelihoodRatio::read(): Given chi-square threshold " << threshold << " is invalid");
172 }
173
174 return in;
175 }
176
177 private:
178
179 double threshold; //!< threshold chi-square to decide if test is passed.
180 };
181}
182
183#endif
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
#define MAKE_STRING(A)
Make string.
Definition JPrint.hh:63
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
std::string getTitle() const
Returns a standard string to be used as title of a graphical root object.
Implementation of the Poisson log-likelihood ratio test.
double threshold
threshold chi-square to decide if test is passed.
std::istream & read(std::istream &in) override
Read test parameters from input.
void test(const TObject *o1, const TObject *o2) override
Applies a Poissonian log-likelihood ratio test to the two given histograms.
Interface to read input and write output for TObject tests.
Definition JTest_t.hh:42
const std::string resultType
test result type
Definition JTest_t.hh:181
const std::string testName
test name
Definition JTest_t.hh:180
Auxiliary class to handle file name, ROOT directory and object name.
Exception for division by zero.
Exception for accessing a value in a collection that is outside of its range.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Structure containing the result of the test.