Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
10 
11 #include "TH1.h"
12 
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"
24 
25 
26 /**
27  * \author bjung
28  */
29 namespace JCOMPAREHISTOGRAMS {
30 
31  /**
32  * Implementation of the Poisson log-likelihood ratio test.\n
33  * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis,\n
34  * which is compared to an alternative hypothesis given by the second histogram.\n\n
35  *
36  * 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
37  */
39  public JTest_t
40  {
41  public:
42 
43  /**
44  * Default constructor.
45  */
47  JTest_t("Poissonian Log-Likelihood Ratio with Beeston-Barlow Correction", "NLLR")
48  {}
49 
50 
51  /**
52  * Applies a Poissonian log-likelihood ratio test to the two given histograms.\n
53  * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis.\n
54  * The second histogram is treated as the expectation for the alternative hypothesis, to which the null hypothesis is compared.
55  *
56  * \param o1 First histogram
57  * \param o2 Second histogram
58  */
59  void test(const TObject* o1, const TObject* o2) override
60  {
61  using namespace std;
62  using namespace JPP;
63 
64  static const double penalty = -1e9; //!< Penalty factor
65 
66  const TH1* h1 = dynamic_cast<const TH1*>(o1);
67  const TH1* h2 = dynamic_cast<const TH1*>(o2);
68 
69  if (h1 == NULL || h2 == NULL) {
70  THROW(JValueOutOfRange, "JTestPoissonLogLikelihood::test(): Could not cast given TObjects to TH1.");
71  }
72 
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);
78  }
79 
80  RooRealVar x("x", "x", h1->GetBinLowEdge(1), h1->GetBinLowEdge(h1->GetNbinsX()) + h1->GetBinWidth(h1->GetNbinsX()));
81 
82  RooDataHist dh1("dh1", h1->GetTitle(), x, h1);
83  RooDataHist dh2("dh2", h2->GetTitle(), x, h2);
84 
85  RooParamHistFunc ph1("ph1", "ph1", dh1);
86  RooHistConstraint hc1("hc1", "hc1", ph1);
87 
88  RooRealSumPdf model_tmp1("model_tmp1", "model_tmp1", ph1, RooRealConstant::value(1.0));
89  RooProdPdf model1 ("model1", "model1", hc1, RooFit::Conditional(model_tmp1,x));
90 
91  RooFitResult* result1 = (RooFitResult*) model1.fitTo(dh2,
92  RooFit::SumW2Error(false),
93  RooFit::PrintLevel(-1),
94  RooFit::Verbose(false),
95  RooFit::Save());
96 
97  if (result1 == NULL) {
98  THROW(JValueOutOfRange, "JTestPoissonLogLikelihoodRatioBeestonBarlow::test(): Unable to retrieve fit results");
99  }
100 
101  double nll2 = 0.0; // Poissonian negative log-likelihood of Asimov data-set onto itself
102 
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) {
106 
107  const double n = h2->GetBinContent(i,j,k);
108 
109  if (n > 0.0) {
110  nll2 += -n*log(n) + log(tgamma(n+1)) + n;
111  } else {
112  nll2 -= penalty;
113  }
114  }
115  }
116  }
117 
118  const double nllratio = 2 * (result1->minNll() - nll2);
119  const bool passed = (nllratio < threshold);
120 
121  const JResultTitle title(testName, resultType, passed, nllratio);
122 
123  TH2* h3 = result1->correlationHist();
124  h3->SetTitle(title.getTitle().c_str());
125 
126  const JTestResult r (testName,
127  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
128  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
129  resultType, nllratio, threshold, h3, passed);
130 
131  this->push_back(r);
132  }
133 
134 
135  /**
136  * Read test parameters from input.
137  *
138  * \param in input stream
139  * \return input stream
140  */
141  std::istream& read(std::istream& in) override
142  {
143  using namespace JPP;
144 
145  in >> threshold;
146 
147  if (!(threshold > 0.0)) {
148  THROW(JValueOutOfRange, "JTestPoissonLogLikelihoodRatioBeestonBarlow::read(): Given chi-square threshold " << threshold << " is invalid");
149  }
150 
151  return in;
152  }
153 
154  private:
155 
156  double threshold; //!< threshold chi-square to decide if test is passed.
157  };
158 }
159 
160 #endif
Exceptions.
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
Interface to read input and write output for TObject tests.
Definition: JTest_t.hh:40
std::string getTitle() const
Returns a standard string to be used as title of a graphical root object.
Definition: JResultTitle.hh:57
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
Definition: JResultTitle.hh:25
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
Definition: JRoot.hh:19
Auxiliary class to handle file name, ROOT directory and object name.
data_type r[M+1]
Definition: JPolint.hh:868
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.
Definition: JPrint.hh:127
const int n
Definition: JPolint.hh:786
const std::string resultType
test result type
Definition: JTest_t.hh:181
const std::string testName
test name
Definition: JTest_t.hh:180
then fatal The output file must have the wildcard in the e g root fi 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
Definition: JCanberra.sh:48
Structure containing the result of the test.
Definition: JTestResult.hh:28
std::istream & read(std::istream &in) override
Read test parameters from input.
int j
Definition: JPolint.hh:792