JTestKolmogorov_t.hh

Go to the documentation of this file.

#ifndef __JTESTKOLMOGOROV_T__
#define __JTESTKOLMOGOROV_T__

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

/**
 * \author rgruiz
 */

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

  /**
   * Kolmogorov 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 Kolmogorov test. The output of a Kolmogorov 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 historgram
   * \param threshold 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 JKolmogorovTest(TH1* h1, TH1* h2, double threshold, std::string testName, std::string parameterName) {    

    using namespace std;
    using namespace JPP;
      
    int n1 = h1 -> GetNbinsX();
    int n2 = h2 -> GetNbinsX();

    if(n1 != n2)
      ERROR("Histograms with different bining. 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())));

    bool afunc1 = false;
    bool afunc2 = false;
    double s1 = 1./h1->Integral();
    double s2 = 1./h2->Integral();

    double ew1, ew2, w1 = 0, w2 = 0;

    for (int bin = 1; bin <= n1; ++bin){
      ew1 = h1->GetBinError(bin);
      ew2 = h2->GetBinError(bin);
      w1 += ew1*ew1;
      w2 += ew2*ew2;
    }
    
    double esum1 = 0, esum2 = 0;
    if (w1 > 0)
       esum1 = 1./s1/s1/w1;
    else
       afunc1 = true;
    if (w2 > 0)
       esum2 = 1./s2/s2/w2;
    else
       afunc2 = true;
    if (afunc2 && afunc1) {
       ERROR("Errors are zero for both histograms\n");
    }
 
    double c1 = 0, c2 = 0;
 
    double dmax = 0;
    for (int bin=1 ; bin<=n1 ; ++bin){
      c1 += s1*h1->GetBinContent(bin);
      c2 += s2*h2->GetBinContent(bin);
      double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
      double p = TMath::KolmogorovProb(d);
      h3->SetBinContent(bin,p);
      dmax = TMath::Max(dmax,TMath::Abs(c1-c2));
    }
 
    double z;
    
    if  (afunc1)
      z = dmax*TMath::Sqrt(esum2);
    else if (afunc2)
      z = dmax*TMath::Sqrt(esum1);
    else
      z = dmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
    
    double pValue = TMath::KolmogorovProb(z);

    bool passed;
    
    (pValue < threshold ? passed = false : passed = true);
    
    JResultTitle title(testName, parameterName, passed , pValue);
    
    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, pValue, threshold, h3, passed);

    return r;
  };


  /**
   * Kolmogorov 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 Kolmogorov test. The output of the Kolmogorov test is a p-value.\n 
   * The parameter threshold should therefore be a real value between 0 and 1.
   * 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 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 slice The axis along which the histogram is sliced.
   *
   * \return Test result.
   */
  JTestResult JKolmogorovTestSlice(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 p = s1 -> KolmogorovTest (s2);

        bool passed;
        
        (p < threshold ? passed = false : passed = true);

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

      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 p = s1 -> KolmogorovTest (s2);

        bool passed;
        
        (p < threshold ? passed = false : passed = true);

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


  
  /**
   * Kolmogorov test for 2D 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 Kolmogorov test. The output of a Kolmogorov 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 historgram
   * \param threshold 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 JKolmogorovTest2D(TH2* h1, TH2* h2, double threshold, std::string testName, std::string parameterName) {    

    using namespace std;
    using namespace JPP;
    
    int n1x = h1 -> GetNbinsX();
    int n2x = h2 -> GetNbinsX();
    int n1y = h1 -> GetNbinsY();
    int n2y = h2 -> GetNbinsY();

    if(n1x != n2x || n1y != n2y)
      ERROR("Histograms with different bining. The objects: " << h1 -> GetName() << " can not be compared." << endl);
     
    if(h1->Integral()==0 || h2->Integral()==0)
      ERROR("Empty histogram: " << h1 -> GetName() << " can not be compared." << endl);

    double s1 = 1./h1->Integral();
    double s2 = 1./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())));

    bool afunc1 = false;
    bool afunc2 = false;

    double ew1, ew2, w1 = 0, w2 = 0;

    for (int i = 1; i <= n1x; ++i){
      for (int j = 1; j <= n1y; ++j){
        ew1 = h1->GetBinError(i,j);
        ew2 = h2->GetBinError(i,j);
        w1 += ew1*ew1;
        w2 += ew2*ew2;
      }
    }
    
    double esum1 = 0, esum2 = 0;
    if (w1 > 0)
       esum1 = 1./s1/s1/w1;
    else
       afunc1 = true;
    if (w2 > 0)
       esum2 = 1./s2/s2/w2;
    else
       afunc2 = true;
    if (afunc2 && afunc1) {
       ERROR("Errors are zero for both histograms\n");
    }
 
    double c1 = 0, c2 = 0;

    double dmax1 = 0;
    for (int i=1 ; i<=n1x ; ++i){
      for (int j=1 ; j<=n1y ; ++j){
        c1 += s1*h1->GetBinContent(i,j);
        c2 += s2*h2->GetBinContent(i,j);
        double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
        h3->Fill(i,j,d);
        dmax1 = TMath::Max(dmax1,TMath::Abs(c1-c2));
      }
    }

    c1 = 0, c2 = 0;
    
    double dmax2 = 0;
    for (int j=1 ; j<=n1y ; ++j){
      for (int i=1 ; i<=n1x ; ++i){
        c1 += s1*h1->GetBinContent(i,j);
        c2 += s2*h2->GetBinContent(i,j);
        double d = TMath::Abs(c1-c2)*TMath::Sqrt(esum1*esum2/(esum1+esum2));
        h3->Fill(i,j,d);
        dmax1 = TMath::Max(dmax2,TMath::Abs(c1-c2));
      }
    }
        
    double dmax = 0.5*(dmax1+dmax2);
    double z;

    if  (afunc1)
      z = dmax*TMath::Sqrt(esum2);
    else if (afunc2)
      z = dmax*TMath::Sqrt(esum1);
    else
      z = dmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
    
    double pValue = TMath::KolmogorovProb(z);

    for (int i=1 ; i<=n1x ; ++i){
      for (int j=1 ; j<=n1y ; ++j){
        h3->SetBinContent(i,j,TMath::KolmogorovProb(0.5 * h3->GetBinContent(i,j)));
      }
    }

    bool passed;
    
    (pValue < threshold ? passed = false : passed = true);
    
    JResultTitle title(testName, parameterName, passed , pValue);
    
    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, pValue, threshold, h3, passed);

    return r;
  };
}; 
  
#endif