Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JTestPoissonLogLikelihoodRatioBeestonBarlow.hh
Go to the documentation of this file.
1#ifndef __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIOBEESTONBARLOW__
2#define __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIOBEESTONBARLOW__
3
4#include <istream>
5#include <ostream>
6
7#include "JLang/JException.hh"
8
11
12#include "TH1.h"
13#include "TH2.h"
14
15#include "RooRealVar.h"
16#include "RooRealConstant.h"
17#include "RooParamHistFunc.h"
18#include "RooHistConstraint.h"
19#include "RooProdPdf.h"
20#include "RooDataSet.h"
21#include "RooDataHist.h"
22#include "RooHistFunc.h"
23#include "RooHistPdf.h"
24#include "RooRealSumPdf.h"
25#include "RooFitResult.h"
26
27
28/**
29 * \author bjung
30 */
31namespace JCOMPAREHISTOGRAMS {
32
33 /**
34 * Implementation of the Poisson log-likelihood ratio test.\n
35 * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis,\n
36 * which is compared to an alternative hypothesis given by the second histogram.\n\n
37 *
38 * This class is derived from the abstract class JTest_t(). For a general description of the implementation of this and other tests derived from JTest_t(), see its documentation.\n
39 */
41 public JTest_t
42 {
43 public:
44
45 /**
46 * Default constructor.
47 */
49 JTest_t("Poissonian Log-Likelihood Ratio with Beeston-Barlow Correction", "NLLR")
50 {}
51
52
53 /**
54 * Applies a Poissonian log-likelihood ratio test to the two given histograms.\n
55 * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis.\n
56 * The second histogram is treated as the expectation for the alternative hypothesis, to which the null hypothesis is compared.
57 *
58 * \param o1 First histogram
59 * \param o2 Second histogram
60 */
61 void test(const TObject* o1, const TObject* o2) override
62 {
63 using namespace std;
64 using namespace JPP;
65
66 static const double penalty = -1e9; //!< Penalty factor
67
68 const TH1* h1 = dynamic_cast<const TH1*>(o1);
69 const TH1* h2 = dynamic_cast<const TH1*>(o2);
70
71 if (h1 == NULL || h2 == NULL) {
72 THROW(JValueOutOfRange, "JTestPoissonLogLikelihood::test(): Could not cast given TObjects to TH1.");
73 }
74
75 if(h1->GetNbinsX() != h2->GetNbinsX() ||
76 h1->GetNbinsY() != h2->GetNbinsY() ||
77 h1->GetNbinsZ() != h2->GetNbinsZ()) {
78 THROW(JValueOutOfRange, "JTestPoissonLogLikelihood::test(): Histograms with different binning. The objects: " <<
79 h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
80 }
81
82 RooRealVar x("x", "x", h1->GetBinLowEdge(1), h1->GetBinLowEdge(h1->GetNbinsX()) + h1->GetBinWidth(h1->GetNbinsX()));
83
84 RooDataHist dh1("dh1", h1->GetTitle(), x, h1);
85 RooDataHist dh2("dh2", h2->GetTitle(), x, h2);
86
87#if ROOT_VERSION_CODE < ROOT_VERSION(6,32,0)
88 RooParamHistFunc ph1("ph1", "ph1", dh1);
89#else
90 RooParamHistFunc ph1("ph1", "ph1", dh1,x);
91#endif
92 RooHistConstraint hc1("hc1", "hc1", ph1);
93
94 RooRealSumPdf model_tmp1("model_tmp1", "model_tmp1", ph1, RooRealConstant::value(1.0));
95 RooProdPdf model1 ("model1", "model1", hc1, RooFit::Conditional(model_tmp1,x));
96
97 RooFitResult* result1 = (RooFitResult*) model1.fitTo(dh2,
98 RooFit::SumW2Error(false),
99 RooFit::PrintLevel(-1),
100 RooFit::Verbose(false),
101 RooFit::Save());
102
103 if (result1 == NULL) {
104 THROW(JValueOutOfRange, "JTestPoissonLogLikelihoodRatioBeestonBarlow::test(): Unable to retrieve fit results");
105 }
106
107 double nll2 = 0.0; // Poissonian negative log-likelihood of Asimov data-set onto itself
108
109 for (int i=1 ; i <= h1->GetNbinsX() ; ++i) {
110 for (int j=1 ; j <= h1->GetNbinsY() ; ++j) {
111 for (int k=1 ; k <= h1->GetNbinsZ() ; ++k) {
112
113 const double n = h2->GetBinContent(i,j,k);
114
115 if (n > 0.0) {
116 nll2 += -n*log(n) + log(tgamma(n+1)) + n;
117 } else {
118 nll2 -= penalty;
119 }
120 }
121 }
122 }
123
124 const double nllratio = 2 * (result1->minNll() - nll2);
125 const bool passed = (nllratio < threshold);
126
127 const JResultTitle title(testName, resultType, passed, nllratio);
128
129 TH2* h3 = result1->correlationHist();
130 h3->SetTitle(title.getTitle().c_str());
131
132 const JTestResult r (testName,
133 JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
134 JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
135 resultType, nllratio, threshold, h3, passed);
136
137 this->push_back(r);
138 }
139
140
141 /**
142 * Read test parameters from input.
143 *
144 * \param in input stream
145 * \return input stream
146 */
147 std::istream& read(std::istream& in) override
148 {
149 using namespace JPP;
150
151 in >> threshold;
152
153 if (!(threshold > 0.0)) {
154 THROW(JValueOutOfRange, "JTestPoissonLogLikelihoodRatioBeestonBarlow::read(): Given chi-square threshold " << threshold << " is invalid");
155 }
156
157 return in;
158 }
159
160 private:
161
162 double threshold; //!< threshold chi-square to decide if test is passed.
163 };
164}
165
166#endif
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
#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.
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 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.