JTestChi2_t.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTCHI2_T__
#define __JCOMPAREHISTOGRAMS__JTESTCHI2_T__

#include <istream>
#include <ostream>

#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<nx1 ; ++i){
        for (int j=1 ; j<ny1 ; ++j){

          double m = h1 -> 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