Jpp  18.0.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JTestRuns_t.hh
Go to the documentation of this file.
1 #ifndef __JCOMPAREHISTOGRAMS__JTESTRUNS_T__
2 #define __JCOMPAREHISTOGRAMS__JTESTRUNS_T__
3 
4 #include <iostream>
5 
8 
9 
10 /**
11  * \author rgruiz
12  */
13 namespace JCOMPAREHISTOGRAMS {
14 
16 
17  /**
18  * Implementation of the different Runs-related tests.
19  */
20  class JTestRuns_t
21  {
22  public:
23  /**
24  * Default constructor.
25  */
27 
28  /**
29  * Implements the Wald-Wolfowitx runs test: <https://en.wikipedia.org/wiki/Wald%E2%80%93Wolfowitz_runs_test>
30  * In this, an expected number of runs and a standard deviation are computed from the number of bins and the number of "aboves" and "belows".
31  * The test returns the difference between the observed number of runs and the expected number of runs, expressed in standard deviations.
32  * This is compared to the threshold input parameter.
33  *
34  * \param h1 First histogram
35  * \param h2 Second histogram
36  * \param threshold
37  * \param parameterName Name of the parameter used to test the histograms
38  * \param testName Name of the test used to compare the histograms
39  */
40  JTestResult JRunsTest(TH1* h1, TH1* h2, double threshold, std::string testName, std::string parameterName) {
41 
42  using namespace std;
43  using namespace JPP;
44 
45  if(h1 -> GetNbinsX() != h2 -> GetNbinsX())
46  ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl);
47 
48  double R = h1->Integral();
49  double T = h2->Integral();
50 
51  TH2D* h3 = (TH2D*)h1->Clone(h1->GetName()==h2->GetName() ?
52  MAKE_CSTRING(to_string(h1->GetName())) :
53  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));
54 
55  for (int i=1 ; i<h1->GetNbinsX() ; ++i) {
56  h3->SetBinContent(i , (T/R)*h1->GetBinContent(i) - h2->GetBinContent(i));
57  }
58 
59  int n = 1;
60  double p = 0;
61  double q = 0;
62 
63  bool a = ((T/R)*h1->GetBinContent(1) - h2->GetBinContent(1)) < 0;
64 
65  (a ? p++ : q++);
66 
67  for (int i = 2 ; i<h1->GetNbinsX() ; ++i){
68 
69  bool b = ((T/R)*h1->GetBinContent(i) - h2->GetBinContent(i)) < 0;
70 
71  (b ? p++ : q++);
72 
73  if (b != a){
74  n++;
75  a=b;
76  }
77  }
78 
79  const double N = 1 + 2*p*q/(p+q) ;
80  const double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) );
81  const double d = (n-N)/s;
82 
83  const bool passed = (fabs(d) > threshold ? false : true);
84 
85  JResultTitle title(testName, parameterName , passed , fabs(d));
86 
87  h3->SetTitle(title.getTitle().c_str());
88 
89  JTestResult r (testName,
90  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
91  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
92  parameterName, fabs(d), threshold, h3, passed);
93 
94  return r;
95  };
96 
97  /**
98  * Runs test for sliced 2D histograms.\n
99  * The histograms are sliced along the axis specified by the slice parameter. A slice per bin is made.\n
100  * For each of the slices, the input parameter threshold, is used to evaluate whether the test is passed or failed.\n
101  * The evaluation is done by comparing the threshold value with the result of the Runs test.\n
102  * The fraction of failed tests is compared to the input parameter failuresThreshold. If this fraction is larger than failuresThreshold, the test fails.
103  *
104  * \param h1 First histogram
105  * \param h2 Second histogram
106  * \param threshold Threshold value for the test result
107  * \param failuresThreshold Threshold value for the fraction of failed tests.
108  * \param testName Name of the test used to compare the histograms
109  * \param parameterName Name of the parameter used to test the histograms
110  * \param slice The axis along which the histogram is sliced.
111  *
112  * \return Test result.
113  */
114  JTestResult JRunsTestSlice(TH2* h1, TH2* h2, double threshold, double failuresThreshold, std::string testName, std::string parameterName, char slice) {
115 
116  using namespace std;
117  using namespace JPP;
118 
119  int nFailures = 0;
120 
121  JTestResult r;
122 
123  if(slice == 'x' || slice == 'X'){
124 
125  int nSlices1 = h1->GetNbinsX();
126  int nSlices2 = h2->GetNbinsX();
127 
128  TH1* h3 = h1->ProjectionX(MAKE_CSTRING(h1->GetName() << "_VS_" <<
129  h2->GetName() << "_RunsTestSliceX"));
130 
131  if (nSlices1 != nSlices2) {
132  ERROR("Histograms with different binning. The objects: " << h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
133  }
134 
135  for (int i=1 ; i<=nSlices1 ; ++i){
136 
137  std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i));
138 
139  TH1D* s1 = h1->ProjectionY (sliceName.c_str(),i,i);
140  TH1D* s2 = h2->ProjectionY (sliceName.c_str(),i,i);
141 
142  if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; }
143 
144  double R = s1->Integral();
145  double T = s2->Integral();
146 
147  int n = 1;
148  double p = 0;
149  double q = 0;
150 
151  bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0;
152 
153  (a ? p++ : q++);
154 
155  for (int i = 2 ; i<s1->GetNbinsX() ; ++i){
156 
157  bool b = ((T/R)*s1->GetBinContent(i) - s2->GetBinContent(i)) < 0;
158 
159  (b ? p++ : q++);
160 
161  if (b != a){
162  n++;
163  a=b;
164  }
165  }
166 
167  double N = 1 + 2*p*q/(p+q) ;
168  double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) );
169  double d = (n-N)/s;
170 
171  bool passed = (fabs(d) > threshold ? false : true);
172 
173  if (!passed) nFailures++;
174 
175  h3->SetBinContent(i,fabs(d));
176 
177  }
178 
179  bool passed = (nFailures/nSlices1 > failuresThreshold ? false : true);
180 
181  JResultTitle title(testName, parameterName, passed , nFailures);
182 
183  h3->SetTitle(title.getTitle().c_str());
184 
185  r = JTestResult(testName,
186  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
187  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
188  parameterName, nFailures, failuresThreshold, h3, passed);
189 
190  }else if (slice == 'y' || slice == 'Y'){
191 
192  int nSlices1 = h1->GetNbinsX();
193  int nSlices2 = h2->GetNbinsX();
194 
195  TH1* h3 = h1->ProjectionY(h1->GetName()==h2->GetName() ?
196  MAKE_CSTRING(to_string(h1->GetName())) :
197  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));
198 
199  if(nSlices1 != nSlices2) {
200  ERROR("Histograms with different binning. The objects: " << h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
201  }
202 
203  for (int i=1 ; i<=nSlices1 ; ++i){
204 
205  std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i));
206 
207  TH1D* s1 = h1->ProjectionX (sliceName.c_str(),i,i);
208  TH1D* s2 = h2->ProjectionX (sliceName.c_str(),i,i);
209 
210  if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; }
211 
212  double R = s1->Integral();
213  double T = s2->Integral();
214 
215  int n = 1;
216  double p = 0;
217  double q = 0;
218 
219  bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0;
220 
221  (a ? p++ : q++);
222 
223  for (int i = 2 ; i<s1->GetNbinsX() ; ++i){
224 
225  bool b = ((T/R)*s1->GetBinContent(i) - s2->GetBinContent(i)) < 0;
226 
227  (b ? p++ : q++);
228 
229  if (b != a){
230  n++;
231  a=b;
232  }
233  }
234 
235  double N = 1 + 2*p*q/(p+q) ;
236  double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) );
237  double d = (n-N)/s;
238 
239  bool passed = (fabs(d) > threshold ? false : true);
240 
241  if (!passed) nFailures++;
242 
243  h3->SetBinContent(i,fabs(d));
244 
245  }
246 
247  bool passed = (nFailures/nSlices1 > failuresThreshold ? false : true);
248 
249  JResultTitle title(testName, parameterName, passed , nFailures);
250 
251  h3->SetTitle(title.getTitle().c_str());
252 
253  r = JTestResult (testName,
254  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
255  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
256  parameterName, nFailures, failuresThreshold, h3, passed);
257 
258  }
259 
260  return r;
261  };
262  };
263 }
264 
265 #endif
JTestResult JRunsTest(TH1 *h1, TH1 *h2, double threshold, std::string testName, std::string parameterName)
Implements the Wald-Wolfowitx runs test: https://en.wikipedia.org/wiki/Wald%E2%80%93Wolfowitz_runs_te...
Definition: JTestRuns_t.hh:40
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
Definition: JResultTitle.hh:22
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:136
Auxiliary class to handle file name, ROOT directory and object name.
data_type r[M+1]
Definition: JPolint.hh:779
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:127
JTestRuns_t()
Default constructor.
Definition: JTestRuns_t.hh:26
const int n
Definition: JPolint.hh:697
then JCalibrateToT a
Definition: JTuneHV.sh:116
do set_variable OUTPUT_DIRECTORY $WORKDIR T
#define ERROR(A)
Definition: JMessage.hh:66
then awk string
Implementation of the different Runs-related tests.
Definition: JTestRuns_t.hh:20
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
std::string to_string(const T &value)
Convert value to string.
Structure containing the result of the test.
Definition: JTestResult.hh:27
JTestResult JRunsTestSlice(TH2 *h1, TH2 *h2, double threshold, double failuresThreshold, std::string testName, std::string parameterName, char slice)
Runs test for sliced 2D histograms.
Definition: JTestRuns_t.hh:114
std::string getTitle()
Returns a standard string to be used as title of a graphical root object.
Definition: JResultTitle.hh:64