JTestChi2_t.hh

Go to the documentation of this file.

#ifndef __JTESTCHI2_T__
#define __JTESTCHI2_T__

#include <istream>
#include <ostream>
#include "JCompareHistograms/JTest_t.hh"

/**
 * \author rgruiz
 */

/**
 * 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,
                   string (h1->GetDirectory()->GetPath()).append(h1->GetName()),
                   string (h2->GetDirectory()->GetPath()).append(h2->GetName()),
                   h1->GetDirectory()->GetFile()->GetName(),
                   h2->GetDirectory()->GetFile()->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(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 chi2 = s1 -> Chi2Test (s2 , options.c_str());

        bool passed;
        
        if (options.find("CHI2") != std::string::npos) {
          (chi2 > threshold ? passed = false : passed = true);
        }else{
          (chi2 < threshold ? passed = false : passed = true);
        }

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

      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->GetNbinsY();
      int nSlices2 = h2->GetNbinsY();

      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() << " 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 chi2 = s1 -> Chi2Test (s2 , options.c_str());

        bool passed;
        
        if (options.find("CHI2") != std::string::npos) {
          (chi2 > threshold ? passed = false : passed = true);
        }else{
          (chi2 < threshold ? passed = false : passed = true);
        }

        if (!passed) nFailures++;
        
        h3->SetBinContent(i,chi2);
      }
      
      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;
  };

  
  /**
   * 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,
                   string (h1->GetDirectory()->GetPath()).append(h1->GetName()),
                   string (h2->GetDirectory()->GetPath()).append(h2->GetName()),
                   h1->GetDirectory()->GetFile()->GetName(),
                   h2->GetDirectory()->GetFile()->GetName(),
                   parameterName, 100*outliers, 100*outliersThreshold, h3, passed);
    
    return r;      
  };
}; 
  
#endif