Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
Public Member Functions | Protected Attributes | Private Attributes | List of all members
JCOMPAREHISTOGRAMS::JTestPoissonLogLikelihoodRatio Class Reference

Implementation of the Poisson log-likelihood ratio test. More...

#include <JTestPoissonLogLikelihoodRatio.hh>

Inheritance diagram for JCOMPAREHISTOGRAMS::JTestPoissonLogLikelihoodRatio:
JCOMPAREHISTOGRAMS::JTest_t std::vector< JTestResult >

Public Member Functions

 JTestPoissonLogLikelihoodRatio ()
 Default constructor.
 
void test (const TObject *o1, const TObject *o2) override
 Applies a Poissonian log-likelihood ratio test to the two given histograms.
 
std::istream & read (std::istream &in) override
 Read test parameters from input.
 
virtual std::ostream & write (std::ostream &out, const char delimiter=' ', const bool onlyFailures=false) const
 Write test result to output.
 
virtual void save (TFile *f, const std::string &path, const bool onlyFailures=false) const
 Writes the test result to root file.
 
const std::string & getTestName () const
 Get test name.
 
const std::string & getResultType () const
 Get result type.
 

Protected Attributes

const std::string testName
 test name
 
const std::string resultType
 test result type
 

Private Attributes

double threshold
 threshold chi-square to decide if test is passed.
 

Detailed Description

Implementation of the Poisson log-likelihood ratio test.

The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis, which is compared to an alternative hypothesis given by the second histogram. The uncertainty on the expectation is taken into account as an additional Poissonian constraint term in the likelihood of each bin ( $\Lambda_i$):

\[ \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}, \]

where:

  • $\gamma_i$ corresponds to a scaling parameter,
  • $m_i$ corresponds to the expectation in bin $i$,
  • $n_i$ corresponds to the observed number of events in bin $i$ of the Asimov dataset and
  • $M_i = \frac{m_i^2}{\sigma_{m_i}^2}$ corresponds to the effective number of exected events in bin $i$,
    where $\sigma_{m_i}$ corresponds to the standard deviation of the expectation $m_i$.

This is inspired on the Beeston-Barlow method, described in the following publication: Computer Physics Communications, Volume 77, Issue 2, 1993 "Fitting using finite Monte Carlo samples", Roger Barlow and Christine Beeston. DOI: 10.1016/0010-4655(93)90005-W

Definition at line 43 of file JTestPoissonLogLikelihoodRatio.hh.

Constructor & Destructor Documentation

◆ JTestPoissonLogLikelihoodRatio()

JCOMPAREHISTOGRAMS::JTestPoissonLogLikelihoodRatio::JTestPoissonLogLikelihoodRatio ( )
inline

Default constructor.

Definition at line 51 of file JTestPoissonLogLikelihoodRatio.hh.

51 :
52 JTest_t("Poisson_NLLR", "NLLR"),
53 threshold(0.0)
54 {}
JCOMPAREHISTOGRAMS::JTestPoissonLogLikelihoodRatio::threshold
double threshold
threshold chi-square to decide if test is passed.
Definition JTestPoissonLogLikelihoodRatio.hh:179
JCOMPAREHISTOGRAMS::JTest_t::JTest_t
JTest_t(const std::string &testName, const std::string &resultType)
Constructor.
Definition JTest_t.hh:51

Member Function Documentation

◆ test()

void JCOMPAREHISTOGRAMS::JTestPoissonLogLikelihoodRatio::test ( const TObject * o1,
const TObject * o2 )
inlineoverridevirtual

Applies a Poissonian log-likelihood ratio test to the two given histograms.


The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis.
The second histogram is treated as the expectation for the alternative hypothesis, to which the null hypothesis is compared.

Parameters
o1First histogram
o2Second histogram

Implements JCOMPAREHISTOGRAMS::JTest_t.

Definition at line 65 of file JTestPoissonLogLikelihoodRatio.hh.

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 }
THROW
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition JException.hh:712
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
MAKE_STRING
#define MAKE_STRING(A)
Make string.
Definition JPrint.hh:63
JCOMPAREHISTOGRAMS::JResultTitle
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
Definition JResultTitle.hh:26
JCOMPAREHISTOGRAMS::JTest_t::resultType
const std::string resultType
test result type
Definition JTest_t.hh:181
JCOMPAREHISTOGRAMS::JTest_t::testName
const std::string testName
test name
Definition JTest_t.hh:180
JGIZMO::JRootObjectID
Auxiliary class to handle file name, ROOT directory and object name.
Definition JRootObjectID.hh:37
JLANG::JDivisionByZero
Exception for division by zero.
Definition JException.hh:288
JLANG::JValueOutOfRange
Exception for accessing a value in a collection that is outside of its range.
Definition JException.hh:180
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JTOOLS::j
int j
Definition JPolint.hh:801
JCOMPAREHISTOGRAMS::JTestResult
Structure containing the result of the test.
Definition JTestResult.hh:30

◆ read()

std::istream & JCOMPAREHISTOGRAMS::JTestPoissonLogLikelihoodRatio::read ( std::istream & in)
inlineoverridevirtual

Read test parameters from input.

Parameters
ininput stream
Returns
input stream

Implements JCOMPAREHISTOGRAMS::JTest_t.

Definition at line 164 of file JTestPoissonLogLikelihoodRatio.hh.

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 }

◆ write()

virtual std::ostream & JCOMPAREHISTOGRAMS::JTest_t::write ( std::ostream & out,
const char delimiter = ' ',
const bool onlyFailures = false ) const
inlinevirtualinherited

Write test result to output.

Parameters
outoutput stream
delimiterfield delimiter
onlyFailuresif true, write only failures.
Returns
output stream

Definition at line 84 of file JTest_t.hh.

87 {
88 using namespace std;
89 using namespace JPP;
90
91 for (vector<JTestResult>::const_iterator r = this->begin() ; r != this->end() ; ++r) {
92
93 if (onlyFailures && r->passed) { continue; }
94
95 print(out, *r, delimiter, true);
96 }
97
98 return out;
99 }
JCOMPAREHISTOGRAMS::print
std::ostream & print(std::ostream &out, const JTestSummary &summary, T __begin, T __end, const bool useColors=true, const JFormat_t &formatting=JFormat_t(18, 3, std::ios::fixed))
Print test summary.
Definition JTestSummary.hh:254

◆ save()

virtual void JCOMPAREHISTOGRAMS::JTest_t::save ( TFile * f,
const std::string & path,
const bool onlyFailures = false ) const
inlinevirtualinherited

Writes the test result to root file.

Parameters
fA ROOT file
pathPath in root file.
onlyFailuresIf true, write only failures.

Definition at line 108 of file JTest_t.hh.

111 {
112 using namespace std;
113 using namespace JPP;
114
115 if (f->GetDirectory(path.c_str())==0) {
116 f->mkdir(path.c_str());
117 }
118
119 f->cd(path.c_str());
120
121 for (vector<JTestResult>::const_iterator r = this->begin() ; r != this->end() ; ++r) {
122
123 if (onlyFailures && r->passed) { continue; }
124
125 r->obj->Write();
126 }
127 }

◆ getTestName()

const std::string & JCOMPAREHISTOGRAMS::JTest_t::getTestName ( ) const
inlineinherited

Get test name.

Returns
test name

Definition at line 135 of file JTest_t.hh.

136 {
137 return testName;
138 }

◆ getResultType()

const std::string & JCOMPAREHISTOGRAMS::JTest_t::getResultType ( ) const
inlineinherited

Get result type.

Returns
result type

Definition at line 146 of file JTest_t.hh.

147 {
148 return resultType;
149 }

Member Data Documentation

◆ threshold

double JCOMPAREHISTOGRAMS::JTestPoissonLogLikelihoodRatio::threshold
private

threshold chi-square to decide if test is passed.

Definition at line 179 of file JTestPoissonLogLikelihoodRatio.hh.

◆ testName

const std::string JCOMPAREHISTOGRAMS::JTest_t::testName
protectedinherited

test name

Definition at line 180 of file JTest_t.hh.

◆ resultType

const std::string JCOMPAREHISTOGRAMS::JTest_t::resultType
protectedinherited

test result type

Definition at line 181 of file JTest_t.hh.

The documentation for this class was generated from the following file:
Generated by doxygen 1.12.0