Applies a Poissonian log-likelihood ratio test to the two given histograms.
62 {
65
66 static const double penalty = -1e9;
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) {
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;
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
128
129 TH2* h3 = result1->correlationHist();
130 h3->SetTitle(title.getTitle().c_str());
131
136
137 this->push_back(r);
138 }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
#define MAKE_STRING(A)
Make string.
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
double threshold
threshold chi-square to decide if test is passed.
const std::string resultType
test result type
const std::string testName
test name
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.