66 static const double penalty = -1e9;
68 const TH1* h1 =
dynamic_cast<const TH1*
>(o1);
69 const TH1* h2 =
dynamic_cast<const TH1*
>(o2);
71 if (h1 == NULL || h2 == NULL) {
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);
82 RooRealVar x(
"x",
"x", h1->GetBinLowEdge(1), h1->GetBinLowEdge(h1->GetNbinsX()) + h1->GetBinWidth(h1->GetNbinsX()));
84 RooDataHist dh1(
"dh1", h1->GetTitle(), x, h1);
85 RooDataHist dh2(
"dh2", h2->GetTitle(), x, h2);
87#if ROOT_VERSION_CODE < ROOT_VERSION(6,32,0)
88 RooParamHistFunc ph1(
"ph1",
"ph1", dh1);
90 RooParamHistFunc ph1(
"ph1",
"ph1", dh1,x);
92 RooHistConstraint hc1(
"hc1",
"hc1", ph1);
94 RooRealSumPdf model_tmp1(
"model_tmp1",
"model_tmp1", ph1, RooRealConstant::value(1.0));
95 RooProdPdf model1 (
"model1",
"model1", hc1, RooFit::Conditional(model_tmp1,x));
97 RooFitResult* result1 = (RooFitResult*) model1.fitTo(dh2,
98 RooFit::SumW2Error(
false),
99 RooFit::PrintLevel(-1),
100 RooFit::Verbose(
false),
103 if (result1 == NULL) {
104 THROW(
JValueOutOfRange,
"JTestPoissonLogLikelihoodRatioBeestonBarlow::test(): Unable to retrieve fit results");
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) {
113 const double n = h2->GetBinContent(i,j,k);
116 nll2 += -n*log(n) + log(tgamma(n+1)) + n;
124 const double nllratio = 2 * (result1->minNll() - nll2);
125 const bool passed = (nllratio <
threshold);
129 TH2* h3 = result1->correlationHist();
130 h3->SetTitle(title.
getTitle().c_str());