#ifndef __JCOMPAREHISTOGRAMS__JTESTCHI2_T__ #define __JCOMPAREHISTOGRAMS__JTESTCHI2_T__ #include #include #include "JCompareHistograms/JTest_t.hh" #include "JCompareHistograms/JResultTitle.hh" /** * \author rgruiz */ namespace JCOMPAREHISTOGRAMS { using JGIZMO::JRootObjectID; /** * Implementation of the different Chi2-related tests. */ class JTestChi2_t { public: /** * Default constructor. */ JTestChi2_t(){} /** * Chi2 test for 1D histograms.\n * The input parameter threshold, is used to evaluate whether the test is passed or failed.\n * The evaluation is done by comparing the threshold value with the result of the Chi2 test. The output of a Chi2 test is a p-value.\n * The parameter threshold should therefore be a real value between 0 and 1. * * \param h1 First histogram * \param h2 Second histogram * \param threshold Threshold value for the test result * \param parameterName Name of the parameter used to test the histograms * \param testName Name of the test used to compare the histograms * \param options ROOT options for the test. * * \return Test result. */ JTestResult JChi2Test(TH1* h1, TH1* h2, double threshold, std::string testName, std::string parameterName, std::string options) { using namespace std; using namespace JPP; if(h1 -> GetNbinsX() != h2 -> GetNbinsX()) ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl); double chi2 = h1 -> Chi2Test (h2 , options.c_str()); double M = h1->Integral(); double N = h2->Integral(); TH2D* h3 = (TH2D*)h1->Clone(h1->GetName()==h2->GetName() ? MAKE_CSTRING(to_string(h1->GetName())) : MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName()))); h3->Reset(); for (int i=1 ; i < h1->GetNbinsX() ; ++i){ double m = h1->GetBinContent(i); double n = h2->GetBinContent(i); if(n!=0 && m!=0){ double c = (M*n - N*m)/sqrt((n+m)*(N*M)); h3->SetBinContent(i,c); } } bool passed; if (options.find("CHI2") != std::string::npos) { (chi2 > threshold ? passed = false : passed = true); }else{ (chi2 < threshold ? passed = false : passed = true); } JResultTitle title(testName, parameterName, passed , chi2); h3->SetTitle(title.getTitle().c_str()); JTestResult r (testName, JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())), JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())), parameterName, chi2, threshold, h3, passed); return r; }; /** * Chi2 test for sliced 2D histograms.\n * The histograms are sliced along the axis specified by the slice parameter. A slice per bin is made.\n * For each of the slices, the input parameter threshold, is used to evaluate whether the test is passed or failed.\n * The evaluation is done by comparing the threshold value with the result of the Chi2 test. The output of a Chi2 test is a p-value.\n * The parameter threshold should therefore be a real value between 0 and 1.\n * The fraction of failed tests is compared to the input parameter failuresThreshold. If this fraction is larger than failuresThreshold, the test fails. * * \param h1 First histogram * \param h2 Second histogram * \param threshold Threshold value for the test result * \param failuresThreshold Threshold value for the fraction of failed tests. * \param parameterName Name of the parameter used to test the histograms * \param testName Name of the test used to compare the histograms * \param options ROOT options for the test. * \param slice The axis along which the histogram is sliced. * * \return Test result. */ JTestResult JChi2TestSlice(TH2* h1, TH2* h2, double threshold, double failuresThreshold, std::string testName, std::string parameterName, std::string options, char slice) { using namespace std; using namespace JPP; int nFailures = 0; JTestResult r; if(slice == 'x' || slice == 'X'){ int nSlices1 = h1->GetNbinsX(); int nSlices2 = h2->GetNbinsX(); TH1* h3 = h1->ProjectionX(MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_Chi2TestSliceX")); if (nSlices1 != nSlices2) { ERROR("Histograms with different binning. The objects: " << h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl); } for (int i=1 ; i<=nSlices1 ; ++i){ std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i)); TH1D* s1 = h1->ProjectionY (sliceName.c_str(),i,i); TH1D* s2 = h2->ProjectionY (sliceName.c_str(),i,i); if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; } double chi2 = s1->Chi2Test(s2 , options.c_str()); bool passed = (options.find("CHI2") != std::string::npos ? (chi2 > threshold ? false : true) : (chi2 < threshold ? false : true)); if (!passed) nFailures++; h3->SetBinContent(i,chi2); } bool passed = (nFailures/nSlices1 > failuresThreshold ? false : true); JResultTitle title(testName, parameterName, passed , nFailures); h3->SetTitle(title.getTitle().c_str()); r = JTestResult(testName, JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())), JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())), parameterName, nFailures, failuresThreshold, h3, passed); } else if (slice == 'y' || slice == 'Y') { int nSlices1 = h1->GetNbinsY(); int nSlices2 = h2->GetNbinsY(); TH1* h3 = h1->ProjectionY(MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_Chi2TestSliceY")); if (nSlices1 != nSlices2) { ERROR("Histograms with different binning. The objects: " << h1->GetName() << " can not be compared." << endl); } for (int i=1 ; i<=nSlices1 ; ++i){ std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i)); TH1D* s1 = h1->ProjectionX (sliceName.c_str(),i,i); TH1D* s2 = h2->ProjectionX (sliceName.c_str(),i,i); if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; } double chi2 = s1 -> Chi2Test (s2 , options.c_str()); bool passed = (options.find("CHI2") != std::string::npos ? (chi2 > threshold ? false : true) : (chi2 < threshold ? false : true));; if (!passed) nFailures++; h3->SetBinContent(i,chi2); } bool passed = (nFailures/nSlices1 > failuresThreshold ? false : true); JResultTitle title(testName, parameterName, passed , nFailures); h3->SetTitle(title.getTitle().c_str()); r = JTestResult (testName, JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())), JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())), parameterName, nFailures, failuresThreshold, h3, passed); } return r; }; /** * Bin-by-Bin Chi2 comparison of 2D histograms.\n * The Chi distance between h1 and h2 is calculated for each bin, and compared to the chi2Threshold() parameter.\n * If the calculated Chi distance is above this threshold, the test is passed for that bin.\n * If the fraction of failures is above the input parameter outliersThreshold(), the test is failed. * * \param h1 First object * \param h2 Second object * \param outliersThreshold fraction of bins allowed to fail the test * \param chi2Threshold p-value * \param parameterName Name of the parameter used to test the histograms * \param testName Name of the test used to compare the histograms * * \return Test result. */ JTestResult JChi2TestBin_2D(TH2* h1, TH2* h2, double outliersThreshold, double chi2Threshold, std::string testName, std::string parameterName) { using namespace std; using namespace JPP; int nx1 = h1->GetNbinsX(); int nx2 = h2->GetNbinsX(); int ny1 = h1->GetNbinsY(); int ny2 = h2->GetNbinsY(); double M = h1->Integral(); double N = h2->Integral(); if(nx1 != nx2 || ny1 != ny2 || M == 0 || N == 0) ERROR("Histograms with different binning. The objects: " << h1->GetName() << " can not be compared." << endl); TH2D* h3 = (TH2D*)h1->Clone(h1->GetName()==h2->GetName() ? MAKE_CSTRING(to_string(h1->GetName())) : MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName()))); h3->Reset(); double outliers = 0; for (int i=1 ; i GetBinContent(i,j); double n = h2 -> GetBinContent(i,j); double chi2 = (n-m*N/M)/sqrt(m*N/M); (fabs(chi2) > chi2Threshold ? outliers+=1./(nx1*ny1) : outliers+=0 ); h3->SetBinContent(i,j,chi2); } } bool passed; (outliers > outliersThreshold ? passed = false : passed = true); JResultTitle title(testName, parameterName , passed , 100*outliers); h3->SetTitle(title.getTitle().c_str()); JTestResult r (testName, JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())), JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())), parameterName, 100*outliers, 100*outliersThreshold, h3, passed); return r; }; }; } #endif