#ifndef __JCOMPAREHISTOGRAMS__JTESTRUNS_T__ #define __JCOMPAREHISTOGRAMS__JTESTRUNS_T__ #include #include "JCompareHistograms/JTest_t.hh" #include "JCompareHistograms/JResultTitle.hh" /** * \author rgruiz */ namespace JCOMPAREHISTOGRAMS { using JGIZMO::JRootObjectID; /** * Implementation of the different Runs-related tests. */ class JTestRuns_t { public: /** * Default constructor. */ JTestRuns_t(){} /** * Implements the Wald-Wolfowitx runs test: * 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". * The test returns the difference between the observed number of runs and the expected number of runs, expressed in standard deviations. * This is compared to the threshold input parameter. * * \param h1 First histogram * \param h2 Second histogram * \param threshold * \param parameterName Name of the parameter used to test the histograms * \param testName Name of the test used to compare the histograms */ JTestResult JRunsTest(TH1* h1, TH1* h2, double threshold, std::string testName, std::string parameterName) { 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 R = h1->Integral(); double T = 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()))); for (int i=1 ; iGetNbinsX() ; ++i) { h3->SetBinContent(i , (T/R)*h1->GetBinContent(i) - h2->GetBinContent(i)); } int n = 1; double p = 0; double q = 0; bool a = ((T/R)*h1->GetBinContent(1) - h2->GetBinContent(1)) < 0; (a ? p++ : q++); for (int i = 2 ; iGetNbinsX() ; ++i){ bool b = ((T/R)*h1->GetBinContent(i) - h2->GetBinContent(i)) < 0; (b ? p++ : q++); if (b != a){ n++; a=b; } } const double N = 1 + 2*p*q/(p+q) ; const double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) ); const double d = (n-N)/s; const bool passed = (fabs(d) > threshold ? false : true); JResultTitle title(testName, parameterName , passed , fabs(d)); 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, fabs(d), threshold, h3, passed); return r; }; /** * Runs 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 Runs test.\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 testName Name of the test used to compare the histograms * \param parameterName Name of the parameter used to test the histograms * \param slice The axis along which the histogram is sliced. * * \return Test result. */ JTestResult JRunsTestSlice(TH2* h1, TH2* h2, double threshold, double failuresThreshold, std::string testName, std::string parameterName, 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() << "_RunsTestSliceX")); 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 R = s1->Integral(); double T = s2->Integral(); int n = 1; double p = 0; double q = 0; bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0; (a ? p++ : q++); for (int i = 2 ; iGetNbinsX() ; ++i){ bool b = ((T/R)*s1->GetBinContent(i) - s2->GetBinContent(i)) < 0; (b ? p++ : q++); if (b != a){ n++; a=b; } } double N = 1 + 2*p*q/(p+q) ; double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) ); double d = (n-N)/s; bool passed = (fabs(d) > threshold ? false : true); if (!passed) nFailures++; h3->SetBinContent(i,fabs(d)); } 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->GetNbinsX(); int nSlices2 = h2->GetNbinsX(); TH1* h3 = h1->ProjectionY(h1->GetName()==h2->GetName() ? MAKE_CSTRING(to_string(h1->GetName())) : MAKE_CSTRING(to_string(h1->GetName()) + "_VS_" + to_string(h2->GetName()))); 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->ProjectionX (sliceName.c_str(),i,i); TH1D* s2 = h2->ProjectionX (sliceName.c_str(),i,i); if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; } double R = s1->Integral(); double T = s2->Integral(); int n = 1; double p = 0; double q = 0; bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0; (a ? p++ : q++); for (int i = 2 ; iGetNbinsX() ; ++i){ bool b = ((T/R)*s1->GetBinContent(i) - s2->GetBinContent(i)) < 0; (b ? p++ : q++); if (b != a){ n++; a=b; } } double N = 1 + 2*p*q/(p+q) ; double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) ); double d = (n-N)/s; bool passed = (fabs(d) > threshold ? false : true); if (!passed) nFailures++; h3->SetBinContent(i,fabs(d)); } 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; }; }; } #endif