Applies Kolmogorov test for two ROOT TH2 histograms.
50 {
53
54 const TH2* h1 = dynamic_cast<const TH2*>(o1);
55 const TH2* h2 = dynamic_cast<const TH2*>(o2);
56
57 if (h1 == NULL || h2 == NULL) {
59 }
60
61 const int n1x = h1->GetNbinsX();
62 const int n2x = h2->GetNbinsX();
63 const int n1y = h1->GetNbinsY();
64 const int n2y = h2->GetNbinsY();
65
66 if(n1x != n2x || n1y != n2y)
67 THROW(
JValueOutOfRange,
"JTestKolmogorov_2D::test(): Histograms with different bining. The objects: " <<
68 h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
69
70 if(h1->Integral() == 0 || h2->Integral() == 0) {
72 h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
73 }
74
75 const double s1 = 1./h1->Integral();
76 const double s2 = 1./h2->Integral();
77
78 TH2* h3 = (TH2*) h1->Clone(h1->GetName() == h2->GetName() ?
81
82 double ew1, ew2, w1 = 0, w2 = 0;
83
84 for (int i = 1; i <= n1x; ++i) {
85 for (
int j = 1;
j <= n1y; ++
j) {
86 ew1 = h1->GetBinError(i,j);
87 ew2 = h2->GetBinError(i,j);
88 w1 += ew1*ew1;
89 w2 += ew2*ew2;
90 }
91 }
92
93 bool afunc1 = false;
94 bool afunc2 = false;
95
96 double esum1 = 0, esum2 = 0;
97
98 if (w1 > 0) {
99 esum1 = 1./s1/s1/w1;
100 } else {
101 afunc1 = true;
102 }
103
104 if (w2 > 0) {
105 esum2 = 1./s2/s2/w2;
106 } else {
107 afunc2 = true;
108 }
109
110 if (afunc2 && afunc1) {
112 }
113
114 double c1 = 0, c2 = 0;
115
116 double dmax1 = 0;
117
118 for (int i=1 ; i<=n1x ; ++i){
119 for (
int j=1 ;
j<=n1y ; ++
j){
120
121 c1 += s1*h1->GetBinContent(i,j);
122 c2 += s2*h2->GetBinContent(i,j);
123
124 double d = TMath::Abs(
c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
125
126 dmax1 = TMath::Max(dmax1,TMath::Abs(
c1-c2));
127
128 h3->Fill(i,j,d);
129 }
130 }
131
133
134 double dmax2 = 0;
135
136 for (
int j=1 ;
j<=n1y ; ++
j){
137 for (int i=1 ; i<=n1x ; ++i){
138
139 c1 += s1*h1->GetBinContent(i,j);
140 c2 += s2*h2->GetBinContent(i,j);
141
142 double d = TMath::Abs(
c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
143
144 dmax2 = TMath::Max(dmax2,TMath::Abs(
c1-c2));
145
146 h3->Fill(i,j,d);
147 }
148 }
149
150 double dmax = 0.5*(dmax1+dmax2);
151
152 double z;
153
154 if (afunc1) {
155 z = dmax*TMath::Sqrt(esum2);
156 } else if (afunc2) {
157 z = dmax*TMath::Sqrt(esum1);
158 } else {
159 z = dmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
160 }
161
162 const double pValue = TMath::KolmogorovProb(z);
163
164 for (int i=1 ; i<=n1x ; ++i) {
165 for (
int j=1 ;
j<=n1y ; ++
j) {
166 h3->SetBinContent(i,j,TMath::KolmogorovProb(0.5 * h3->GetBinContent(i,j)));
167 }
168 }
169
170 const bool passed = (pValue >
threshold);
171
173
174 h3->SetTitle(title.getTitle().c_str());
176
181
182 this->push_back(r);
183 }
TCanvas * c1
Global variables to handle mouse events.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
#define MAKE_CSTRING(A)
Make C-string.
#define MAKE_STRING(A)
Make string.
Class dedicated to standardize the title of the graphical objects produced by the JTest_t() derived c...
double threshold
threshold p-value to decide if test is passed.
const std::string resultType
test result type
const std::string testName
test name
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.