1 #ifndef __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIOBEESTONBARLOW__
2 #define __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIOBEESTONBARLOW__
13 #include "RooRealVar.h"
14 #include "RooRealConstant.h"
15 #include "RooParamHistFunc.h"
16 #include "RooHistConstraint.h"
17 #include "RooProdPdf.h"
18 #include "RooDataSet.h"
19 #include "RooDataHist.h"
20 #include "RooHistFunc.h"
21 #include "RooHistPdf.h"
22 #include "RooRealSumPdf.h"
23 #include "RooFitResult.h"
29 namespace JCOMPAREHISTOGRAMS {
47 JTest_t(
"Poissonian Log-Likelihood Ratio with Beeston-Barlow Correction",
"NLLR")
64 static const double penalty = -1e9;
66 const TH1* h1 =
dynamic_cast<const TH1*
>(o1);
67 const TH1* h2 =
dynamic_cast<const TH1*
>(o2);
69 if (h1 == NULL || h2 == NULL) {
70 THROW(JValueOutOfRange,
"JTestPoissonLogLikelihood::test(): Could not cast given TObjects to TH1.");
73 if(h1->GetNbinsX() != h2->GetNbinsX() ||
74 h1->GetNbinsY() != h2->GetNbinsY() ||
75 h1->GetNbinsZ() != h2->GetNbinsZ()) {
76 THROW(JValueOutOfRange,
"JTestPoissonLogLikelihood::test(): Histograms with different binning. The objects: " <<
77 h1->GetName() <<
" and " << h2->GetName() <<
" can not be compared." << endl);
80 RooRealVar
x(
"x",
"x", h1->GetBinLowEdge(1), h1->GetBinLowEdge(h1->GetNbinsX()) + h1->GetBinWidth(h1->GetNbinsX()));
82 RooDataHist dh1(
"dh1", h1->GetTitle(),
x, h1);
83 RooDataHist dh2(
"dh2", h2->GetTitle(),
x, h2);
85 RooParamHistFunc ph1(
"ph1",
"ph1", dh1);
86 RooHistConstraint hc1(
"hc1",
"hc1", ph1);
88 RooRealSumPdf model_tmp1(
"model_tmp1",
"model_tmp1", ph1, RooRealConstant::value(1.0));
89 RooProdPdf model1 (
"model1",
"model1", hc1, RooFit::Conditional(model_tmp1,x));
91 RooFitResult* result1 = (RooFitResult*) model1.fitTo(dh2,
92 RooFit::SumW2Error(
false),
93 RooFit::PrintLevel(-1),
94 RooFit::Verbose(
false),
97 if (result1 == NULL) {
98 THROW(JValueOutOfRange,
"JTestPoissonLogLikelihoodRatioBeestonBarlow::test(): Unable to retrieve fit results");
103 for (
int i=1 ;
i <= h1->GetNbinsX() ; ++
i) {
104 for (
int j=1 ;
j <= h1->GetNbinsY() ; ++
j) {
105 for (
int k=1 ;
k <= h1->GetNbinsZ() ; ++
k) {
107 const double n = h2->GetBinContent(
i,
j,
k);
110 nll2 += -n*
log(n) +
log(tgamma(n+1)) +
n;
118 const double nllratio = 2 * (result1->minNll() - nll2);
119 const bool passed = (nllratio <
threshold);
123 TH2* h3 = result1->correlationHist();
124 h3->SetTitle(title.
getTitle().c_str());
141 std::istream&
read(std::istream&
in)
override
148 THROW(JValueOutOfRange,
"JTestPoissonLogLikelihoodRatioBeestonBarlow::read(): Given chi-square threshold " <<
threshold <<
" is invalid");
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
Interface to read input and write output for TObject tests.
std::string getTitle() const
Returns a standard string to be used as title of a graphical root object.
then set_variable DIR else fatal Wrong number of arguments fi for INPUT_FILE in ls rt $DIR stage * log
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Auxiliary class to handle file name, ROOT directory and object name.
Implementation of the Poisson log-likelihood ratio test.
void test(const TObject *o1, const TObject *o2) override
Applies a Poissonian log-likelihood ratio test to the two given histograms.
#define MAKE_STRING(A)
Make string.
double threshold
threshold chi-square to decide if test is passed.
const std::string resultType
test result type
const std::string testName
test name
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Structure containing the result of the test.
std::istream & read(std::istream &in) override
Read test parameters from input.
JTestPoissonLogLikelihoodRatioBeestonBarlow()
Default constructor.