Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JTestSignificance.hh
Go to the documentation of this file.
1#ifndef __JCOMPAREHISTOGRAMS__JTESTSIGNIFICANCE__
2#define __JCOMPAREHISTOGRAMS__JTESTSIGNIFICANCE__
3
4#include <istream>
5#include <ostream>
6
7#include "JLang/JException.hh"
8
11
12#include "TH1.h"
13
14
15/**
16 * \author rgruiz, bofearreigh, bjung
17 */
18namespace JCOMPAREHISTOGRAMS {}
19namespace JPP { using namespace JCOMPAREHISTOGRAMS; }
20
21namespace JCOMPAREHISTOGRAMS {
22
23 /**
24 * Significance test
25 *
26 * Compares two histograms \f$ H_a \f$ and \f$ H_b \f$ by calculating the average bin significance \n
27 *
28 * \f[
29 * S = \frac{1}{N} \cdot \sum_{i=0}^{N} \frac{|a_i - K \cdot b_i|}{\sqrt{\sigma ^2(a_i) + K^2 \cdot \sigma ^2(b_i)}}
30 * \f]
31 *
32 * where:
33 * - \f$ N \f$ is the number of bins in \f$ H_a \f$ and \f$ H_b \f$.
34 * - \f$ a_i \f$ and \f$ b_i \f$ are the contents of bin \f$ i \f$ in \f$ H_a \f$ and \f$ H_b \f$.
35 * - \f$ K \f$ is a normalisation parameter which can be specified by the user.
36 * - It should be interpreted as the integral of \f$ H_a \f$ divided by the integral of \f$ H_b \f$.
37 * - Note that if the specified value of K is negative, the value will be reset to the integral of \f$ H_a \f$ divided by the integral of \f$ H_b \f$.
38 * - \f$ \sigma \f$ is calculated assuming that the bin contents follow Poisson distributions.
39 *
40 * A threshold value for \f$ S \f$ is used to tag the test as passed or failed.
41 */
43 public JTest_t
44 {
45 public:
46
47 /**
48 * Default constructor.
49 */
51 JTest_t("Significance", "Significance")
52 {}
53
54
55 /**
56 * Applies Significance test for two ROOT TH1 histograms.
57 *
58 * \param o1 First histogram
59 * \param o2 Second histogram
60 */
61 void test(const TObject* o1, const TObject* o2) override
62 {
63 using namespace std;
64 using namespace JPP;
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, "JTestSignificance::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, "JTestSignificance::test(): Histograms with different bining. The objects: " <<
77 h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
78 }
79
80 if (K < 0) {
81 K = h2->GetEntries() / h1->GetEntries();
82 }
83
84 TH1* h3 = (TH1*) h1->Clone(h1->GetName() == h2->GetName() ?
85 MAKE_CSTRING(h1->GetName() << "_" << testName) :
86 MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_" << testName));
87
88 h3->Add(h2,-1*K);
89
90 double S = 0;
91
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){
95
96 const double a = h1->GetBinContent(i,j,k);
97 const double b = h2->GetBinContent(i,j,k);
98
99 if (a!=0 || b!=0) {
100 S += fabs((a - K*b)/sqrt(a + K*K*b));
101 }
102 }
103 }
104 }
105
106 S /= (h1->GetNbinsX()*h1->GetNbinsY()*h1->GetNbinsZ());
107
108 const bool passed = (S > threshold);
109
110 const JResultTitle title(testName, resultType, passed, S);
111
112 h3->SetTitle(title.getTitle().c_str());
113
114 const int Ndims = h3->GetDimension();
115
116 if (Ndims == 1) {
117 h3->GetYaxis()->SetTitle(resultType.c_str());
118 } else if (Ndims == 2) {
119 h3->GetZaxis()->SetTitle(resultType.c_str());
120 }
121
122 const JTestResult r(testName,
123 JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
124 JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
125 resultType, S, threshold, h3, passed);
126
127 this->push_back(r);
128 }
129
130
131 /**
132 * Read test parameters from input.
133 *
134 * \param in input stream
135 * \return input stream
136 */
137 std::istream& read(std::istream& in) override
138 {
139 using namespace JPP;
140
141 in >> threshold >> K;
142
143 if (threshold < 0.0) {
144 THROW(JValueOutOfRange, "JTestSignificance::read(): Invalid threshold value " << threshold);
145 }
146
147 return in;
148 }
149
150 private:
151
152 double threshold; //!< threshold p-value to decide if test is passed.
153 double K; //!< normalization factor between histograms.
154 };
155}
156
157#endif
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
#define MAKE_CSTRING(A)
Make C-string.
Definition JPrint.hh:72
#define MAKE_STRING(A)
Make string.
Definition JPrint.hh:63
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
std::string getTitle() const
Returns a standard string to be used as title of a graphical root object.
void test(const TObject *o1, const TObject *o2) override
Applies Significance test for two ROOT TH1 histograms.
double K
normalization factor between histograms.
std::istream & read(std::istream &in) override
Read test parameters from input.
double threshold
threshold p-value to decide if test is passed.
Interface to read input and write output for TObject tests.
Definition JTest_t.hh:42
const std::string resultType
test result type
Definition JTest_t.hh:181
const std::string testName
test name
Definition JTest_t.hh:180
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.