JTestRuns_t.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTRUNS_T__
#define __JCOMPAREHISTOGRAMS__JTESTRUNS_T__

#include <iostream>

#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: <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;
        }
      }

      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 ; 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 ? 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 ; 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 ? 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