Jpp  15.0.1-rc.1-highQE
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JTestKolmogorov_t.hh
Go to the documentation of this file.
1 #ifndef __JCOMPAREHISTOGRAMS__JTESTKOLMOGOROV_T__
2 #define __JCOMPAREHISTOGRAMS__JTESTKOLMOGOROV_T__
3 
4 #include <istream>
5 #include <ostream>
6 
9 
10 
11 /**
12  * \author rgruiz
13  */
14 namespace JCOMPAREHISTOGRAMS {
15 
17 
18  /**
19  * Implementation of the different Kolmogorov-related tests.
20  */
22  {
23  public:
24  /**
25  * Default constructor.
26  */
28 
29  /**
30  * Kolmogorov test for 1D histograms.\n
31  * The input parameter threshold, is used to evaluate whether the test is passed or failed.\n
32  * The evaluation is done by comparing the threshold value with the result of the Kolmogorov test. The output of a Kolmogorov test is a p-value.\n
33  * The parameter threshold should therefore be a real value between 0 and 1.
34  *
35  * \param h1 First histogram
36  * \param h2 Second historgram
37  * \param threshold p-value
38  * \param parameterName Name of the parameter used to test the histograms
39  * \param testName Name of the test used to compare the histograms
40  *
41  * \return Test result
42  */
43  JTestResult JKolmogorovTest(TH1* h1, TH1* h2, double threshold, std::string testName, std::string parameterName) {
44 
45  using namespace std;
46  using namespace JPP;
47 
48  int n1 = h1 -> GetNbinsX();
49  int n2 = h2 -> GetNbinsX();
50 
51  if(n1 != n2)
52  ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl);
53 
54  TH2D* h3 = (TH2D*)h1->Clone(h1->GetName()==h2->GetName() ?
55  MAKE_CSTRING(to_string(h1->GetName())) :
56  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));
57 
58  bool afunc1 = false;
59  bool afunc2 = false;
60  double s1 = 1./h1->Integral();
61  double s2 = 1./h2->Integral();
62 
63  double ew1, ew2, w1 = 0, w2 = 0;
64 
65  for (int bin = 1; bin <= n1; ++bin){
66  ew1 = h1->GetBinError(bin);
67  ew2 = h2->GetBinError(bin);
68  w1 += ew1*ew1;
69  w2 += ew2*ew2;
70  }
71 
72  double esum1 = 0, esum2 = 0;
73  if (w1 > 0)
74  esum1 = 1./s1/s1/w1;
75  else
76  afunc1 = true;
77  if (w2 > 0)
78  esum2 = 1./s2/s2/w2;
79  else
80  afunc2 = true;
81  if (afunc2 && afunc1) {
82  ERROR("Errors are zero for both histograms\n");
83  }
84 
85  double c1 = 0, c2 = 0;
86 
87  double dmax = 0;
88  for (int bin=1 ; bin<=n1 ; ++bin){
89  c1 += s1*h1->GetBinContent(bin);
90  c2 += s2*h2->GetBinContent(bin);
91  double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
92  double p = TMath::KolmogorovProb(d);
93  h3->SetBinContent(bin,p);
94  dmax = TMath::Max(dmax,TMath::Abs(c1-c2));
95  }
96 
97  double z;
98 
99  if (afunc1)
100  z = dmax*TMath::Sqrt(esum2);
101  else if (afunc2)
102  z = dmax*TMath::Sqrt(esum1);
103  else
104  z = dmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
105 
106  double pValue = TMath::KolmogorovProb(z);
107 
108  bool passed;
109 
110  (pValue < threshold ? passed = false : passed = true);
111 
112  JResultTitle title(testName, parameterName, passed , pValue);
113 
114  h3->SetTitle(title.getTitle().c_str());
115 
116  JTestResult r (testName,
117  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
118  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
119  parameterName, pValue, threshold, h3, passed);
120 
121  return r;
122  };
123 
124 
125  /**
126  * Kolmogorov test for sliced 2D histograms.\n
127  * The histograms are sliced along the axis specified by the slice parameter. A slice per bin is made.\n
128  * For each of the slices, the input parameter threshold is used to evaluate whether the test is passed or failed.\n
129  * The evaluation is done by comparing the threshold value with the result of the Kolmogorov test. The output of the Kolmogorov test is a p-value.\n
130  * The parameter threshold should therefore be a real value between 0 and 1.
131  * The fraction of failed tests is compared to the input parameter failuresThreshold. If this fraction is larger than failuresThreshold, the test fails.
132  *
133  * \param h1 First histogram
134  * \param h2 Second histogram
135  * \param threshold Threshold value for the test result
136  * \param failuresThreshold Threshold value for the test result
137  * \param parameterName Name of the parameter used to test the histograms
138  * \param testName Name of the test used to compare the histograms
139  * \param slice The axis along which the histogram is sliced.
140  *
141  * \return Test result.
142  */
143  JTestResult JKolmogorovTestSlice(TH2* h1, TH2* h2, double threshold, double failuresThreshold, std::string testName, std::string parameterName, char slice) {
144 
145  using namespace std;
146  using namespace JPP;
147 
148  int nFailures = 0;
149 
150  JTestResult r;
151 
152  if(slice == 'x' || slice == 'X'){
153 
154  int nSlices1 = h1->GetNbinsX();
155  int nSlices2 = h2->GetNbinsX();
156 
157  TH1* h3 = h1->ProjectionX(h1->GetName()==h2->GetName() ?
158  MAKE_CSTRING(to_string(h1->GetName())) :
159  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));
160 
161  if(nSlices1 != nSlices2)
162  ERROR("Histograms with different binning. The objects: " << h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
163 
164  for (int i=1 ; i<=nSlices1 ; ++i){
165 
166  std::string sliceName = MAKE_STRING(h1->GetName() + to_string("_") + to_string(i));
167 
168  TH1D* s1 = h1->ProjectionY (sliceName.c_str(),i,i);
169  TH1D* s2 = h2->ProjectionY (sliceName.c_str(),i,i);
170 
171  double p = s1 -> KolmogorovTest (s2);
172 
173  bool passed;
174 
175  (p < threshold ? passed = false : passed = true);
176 
177  if (!passed) nFailures++;
178 
179  h3->SetBinContent(i,p);
180 
181  }
182 
183  bool passed;
184 
185  (nFailures/nSlices1 > failuresThreshold ? passed = false : passed = true);
186 
187  JResultTitle title(testName, parameterName, passed , nFailures);
188 
189  h3->SetTitle(title.getTitle().c_str());
190 
191  r = JTestResult (testName,
192  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
193  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
194  parameterName, nFailures, failuresThreshold, h3, passed);
195 
196  }else if (slice == 'y' || slice == 'Y'){
197 
198  int nSlices1 = h1->GetNbinsY();
199  int nSlices2 = h2->GetNbinsY();
200 
201  TH1* h3 = h1->ProjectionY(h1->GetName()==h2->GetName() ?
202  MAKE_CSTRING(to_string(h1->GetName())) :
203  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));
204 
205  if(nSlices1 != nSlices2)
206  ERROR("Histograms with different binning. The objects: " << h1->GetName() << " can not be compared." << endl);
207 
208  for (int i=1 ; i<=nSlices1 ; ++i){
209 
210  std::string sliceName = MAKE_STRING(h1->GetName() + to_string("_") + to_string(i));
211 
212  TH1D* s1 = h1->ProjectionX (sliceName.c_str(),i,i);
213  TH1D* s2 = h2->ProjectionX (sliceName.c_str(),i,i);
214 
215  double p = s1 -> KolmogorovTest (s2);
216 
217  bool passed;
218 
219  (p < threshold ? passed = false : passed = true);
220 
221  if (!passed) nFailures++;
222 
223  h3->SetBinContent(i,p);
224  }
225 
226  bool passed;
227 
228  (nFailures/nSlices1 > failuresThreshold ? passed = false : passed = true);
229 
230  JResultTitle title(testName, parameterName, passed , nFailures);
231 
232  h3->SetTitle(title.getTitle().c_str());
233 
234  r = JTestResult (testName,
235  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
236  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
237  parameterName, nFailures, failuresThreshold, h3, passed);
238 
239  }
240 
241  return r;
242  };
243 
244 
245 
246  /**
247  * Kolmogorov test for 2D histograms.\n
248  * The input parameter threshold, is used to evaluate whether the test is passed or failed.\n
249  * The evaluation is done by comparing the threshold value with the result of the Kolmogorov test. The output of a Kolmogorov test is a p-value.\n
250  * The parameter threshold should therefore be a real value between 0 and 1.
251  *
252  * \param h1 First histogram
253  * \param h2 Second historgram
254  * \param threshold p-value
255  * \param parameterName Name of the parameter used to test the histograms
256  * \param testName Name of the test used to compare the histograms
257  *
258  * \return Test result
259  */
260  JTestResult JKolmogorovTest2D(TH2* h1, TH2* h2, double threshold, std::string testName, std::string parameterName) {
261 
262  using namespace std;
263  using namespace JPP;
264 
265  int n1x = h1 -> GetNbinsX();
266  int n2x = h2 -> GetNbinsX();
267  int n1y = h1 -> GetNbinsY();
268  int n2y = h2 -> GetNbinsY();
269 
270  if(n1x != n2x || n1y != n2y)
271  ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl);
272 
273  if(h1->Integral()==0 || h2->Integral()==0)
274  ERROR("Empty histogram: " << h1 -> GetName() << " can not be compared." << endl);
275 
276  double s1 = 1./h1->Integral();
277  double s2 = 1./h2->Integral();
278 
279  TH2D* h3 = (TH2D*)h1->Clone(h1->GetName()==h2->GetName() ?
280  MAKE_CSTRING(to_string(h1->GetName())) :
281  MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName())));
282 
283  bool afunc1 = false;
284  bool afunc2 = false;
285 
286  double ew1, ew2, w1 = 0, w2 = 0;
287 
288  for (int i = 1; i <= n1x; ++i){
289  for (int j = 1; j <= n1y; ++j){
290  ew1 = h1->GetBinError(i,j);
291  ew2 = h2->GetBinError(i,j);
292  w1 += ew1*ew1;
293  w2 += ew2*ew2;
294  }
295  }
296 
297  double esum1 = 0, esum2 = 0;
298  if (w1 > 0)
299  esum1 = 1./s1/s1/w1;
300  else
301  afunc1 = true;
302  if (w2 > 0)
303  esum2 = 1./s2/s2/w2;
304  else
305  afunc2 = true;
306  if (afunc2 && afunc1) {
307  ERROR("Errors are zero for both histograms\n");
308  }
309 
310  double c1 = 0, c2 = 0;
311 
312  double dmax1 = 0;
313  for (int i=1 ; i<=n1x ; ++i){
314  for (int j=1 ; j<=n1y ; ++j){
315  c1 += s1*h1->GetBinContent(i,j);
316  c2 += s2*h2->GetBinContent(i,j);
317  double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
318  h3->Fill(i,j,d);
319  dmax1 = TMath::Max(dmax1,TMath::Abs(c1-c2));
320  }
321  }
322 
323  c1 = 0, c2 = 0;
324 
325  double dmax2 = 0;
326  for (int j=1 ; j<=n1y ; ++j){
327  for (int i=1 ; i<=n1x ; ++i){
328  c1 += s1*h1->GetBinContent(i,j);
329  c2 += s2*h2->GetBinContent(i,j);
330  double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
331  h3->Fill(i,j,d);
332  dmax1 = TMath::Max(dmax2,TMath::Abs(c1-c2));
333  }
334  }
335 
336  double dmax = 0.5*(dmax1+dmax2);
337  double z;
338 
339  if (afunc1)
340  z = dmax*TMath::Sqrt(esum2);
341  else if (afunc2)
342  z = dmax*TMath::Sqrt(esum1);
343  else
344  z = dmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
345 
346  double pValue = TMath::KolmogorovProb(z);
347 
348  for (int i=1 ; i<=n1x ; ++i){
349  for (int j=1 ; j<=n1y ; ++j){
350  h3->SetBinContent(i,j,TMath::KolmogorovProb(0.5 * h3->GetBinContent(i,j)));
351  }
352  }
353 
354  bool passed;
355 
356  (pValue < threshold ? passed = false : passed = true);
357 
358  JResultTitle title(testName, parameterName, passed , pValue);
359 
360  h3->SetTitle(title.getTitle().c_str());
361 
362  JTestResult r (testName,
363  JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
364  JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
365  parameterName, pValue, threshold, h3, passed);
366 
367  return r;
368  };
369  };
370 }
371 
372 #endif
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
Definition: JResultTitle.hh:22
Implementation of the different Kolmogorov-related tests.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:71
Auxiliary class to handle file name, ROOT directory and object name.
data_type r[M+1]
Definition: JPolint.hh:742
JTestResult JKolmogorovTest(TH1 *h1, TH1 *h2, double threshold, std::string testName, std::string parameterName)
Kolmogorov test for 1D histograms.
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:142
#define ERROR(A)
Definition: JMessage.hh:66
JTestResult JKolmogorovTestSlice(TH2 *h1, TH2 *h2, double threshold, double failuresThreshold, std::string testName, std::string parameterName, char slice)
Kolmogorov test for sliced 2D histograms.
TCanvas * c1
Global variables to handle mouse events.
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 JKolmogorovTest2D(TH2 *h1, TH2 *h2, double threshold, std::string testName, std::string parameterName)
Kolmogorov test for 2D histograms.
int j
Definition: JPolint.hh:666
std::string getTitle()
Returns a standard string to be used as title of a graphical root object.
Definition: JResultTitle.hh:64