JTestRuns_t.hh

Go to the documentation of this file.

#ifndef __JTESTRUNS_T__
#define __JTESTRUNS_T__

#include <iostream>
#include "JCompareHistograms/JTest_t.hh"

/**
 *  \author rgruiz
 */

/**
 * Implementation of the different Runs-related tests.
 */
class JTestRuns_t 
{ 
public:  
  /**
   * Default constructor.
   */
  JTestRuns_t(){}

  /**
   * Implements the Wald-Wolfowitx runs test: <https://en.wikipedia.org/wiki/Wald%E2%80%93Wolfowitz_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 ; i<h1->GetNbinsX() ; ++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 ; i<h1->GetNbinsX() ; ++i){
      
      bool b = ((T/R)*h1->GetBinContent(i) - h2->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 ? passed = false : passed = true);

    JResultTitle title(testName, parameterName , passed , fabs(d));
    
    h3->SetTitle(title.getTitle().c_str());    
    
    JTestResult r (testName,
                   string (h1->GetDirectory()->GetPath()).append(h1->GetName()),
                   string (h2->GetDirectory()->GetPath()).append(h2->GetName()),
                   h1->GetDirectory()->GetFile()->GetName(),
                   h2->GetDirectory()->GetFile()->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(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(h1->GetName() + to_string("_") + to_string(i));
             
        TH1D* s1 = h1->ProjectionY (sliceName.c_str(),i,i);
        TH1D* s2 = h2->ProjectionY (sliceName.c_str(),i,i);
                
        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 ; i<s1->GetNbinsX() ; ++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 ? passed = false : passed = true);

        if (!passed) nFailures++;
        
        h3->SetBinContent(i,fabs(d));
        
      }

      bool passed;

      (nFailures/nSlices1 > failuresThreshold ? passed = false : passed = true);
      
      JResultTitle title(testName, parameterName, passed , nFailures);
      
      h3->SetTitle(title.getTitle().c_str());
      
      r = JTestResult  (testName,
                        string (h1->GetDirectory()->GetPath()).append(h1->GetName()),
                        string (h2->GetDirectory()->GetPath()).append(h2->GetName()),
                        h1->GetDirectory()->GetFile()->GetName(),
                        h2->GetDirectory()->GetFile()->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(h1->GetName() + to_string("_") + to_string(i));
             
        TH1D* s1 = h1->ProjectionX (sliceName.c_str(),i,i);
        TH1D* s2 = h2->ProjectionX (sliceName.c_str(),i,i);
                
        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 ; i<s1->GetNbinsX() ; ++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 ? passed = false : passed = true);

        if (!passed) nFailures++;
        
        h3->SetBinContent(i,fabs(d));
        
      }

      bool passed;

      (nFailures/nSlices1 > failuresThreshold ? passed = false : passed = true);
      
      JResultTitle title(testName, parameterName, passed , nFailures);
      
      h3->SetTitle(title.getTitle().c_str());
      
      r = JTestResult  (testName,
                        string (h1->GetDirectory()->GetPath()).append(h1->GetName()),
                        string (h2->GetDirectory()->GetPath()).append(h2->GetName()),
                        h1->GetDirectory()->GetFile()->GetName(),
                        h2->GetDirectory()->GetFile()->GetName(),
                        parameterName, nFailures, failuresThreshold, h3, passed);

    }

    return r;
  };  
}; 
  
#endif