70 const TH1* h1 =
dynamic_cast<const TH1*
>(o1);
71 const TH1* h2 =
dynamic_cast<const TH1*
>(o2);
73 if (h1 == NULL || h2 == NULL) {
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);
84 TH1* h3 = (TH1*) h1->Clone(h1->GetName() == h2->GetName() ?
90 double LnLratio = 0.0;
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) {
96 double contribution = 0.0;
98 const double n1 = h1->GetBinContent(i,j,k);
99 const double n2 = h2->GetBinContent(i,j,k);
101 const double e1 = h1->GetBinError(i,j,k);
102 const double e2 = h2->GetBinError(i,j,k);
107 THROW(
JDivisionByZero,
"JTestPoissonLogLikelihoodRatio::test(): Invalid effective bin content for histogram " << h2->GetName() <<
" at bin (" << i <<
',' << j <<
',' << k <<
')');
110 const double N2 = n2*n2/e2/e2;
115 THROW(
JDivisionByZero,
"JTestPoissonLogLikelihoodRatio::test(): Invalid effective bin content for histogram " << h1->GetName() <<
" at bin (" << i <<
',' << j <<
',' << k <<
')');
118 const double f1 = (n1 + N2) / (n2 + N2);
120 const double N1 = n1*n1/e1/e1;
122 contribution = 2 * (N2 - N1 +
123 lgamma(N2 + 1) - lgamma(N1 + 1) +
124 n1*log(n1/n2/f1) - N2*log(f1*N2) +
129 const double f1 = N2 / (n2 + N2);
131 contribution = 2 * (N2 + lgamma(N2 + 1) - N2*log(f1*N2));
135 h3->SetBinContent(i,j,k, contribution);
137 LnLratio += contribution;
142 const bool passed = (LnLratio <
threshold);
146 h3->SetTitle(title.
getTitle().c_str());