JTestKolmogorov_1D.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTKOLMOGOROV_1D__
#define __JCOMPAREHISTOGRAMS__JTESTKOLMOGOROV_1D__
 
#include <istream>
#include <ostream>
 
#include "JLang/JException.hh"
 
#include "JCompareHistograms/JTest_t.hh"
#include "JCompareHistograms/JResultTitle.hh"
 
#include "TH1.h"
#include "TMath.h"
 
 
/**
 * \author rgruiz, bjung
 */
namespace JCOMPAREHISTOGRAMS {}
namespace JPP { using namespace JCOMPAREHISTOGRAMS; }
 
namespace JCOMPAREHISTOGRAMS {
 
  /**
   * 
   * Implementation of the Kolmogorov test for 1D histograms.\n
   * This class is derived from the abstract class JTest_t(). For a general description of the implementation of this and other tests derived from JTest_t(), see its documentation.\n
   * This test uses the input parameter threshold to evaluate whether the test is passed or failed.\n 
   * The evaluation is done by comparing the threshold value with the result of the JTestKolmogorov_1D() test. This is a p-value.\n 
   * The parameter threshold should therefore be a real value between 0 and 1.
   */
  class JTestKolmogorov_1D:
    public JTest_t
  { 
  public:  
 
    /**
     * Default constructor.
     */ 
    JTestKolmogorov_1D() :
      JTest_t("KS_1D", "p-Value(KS)")
    {}
 
    
    /**
     * Applies Kolmogorov test for two ROOT TH1 histograms.
     *
     * \param o1 First histogram
     * \param o2 Second histogram
     */
    void test(const TObject* o1, const TObject* o2) override
    {    
      using namespace std;
      using namespace JPP;
 
      const TH1* h1 =  dynamic_cast<const TH1*>(o1);
      const TH1* h2 =  dynamic_cast<const TH1*>(o2);
    
      if (h1->GetDimension() != 1 || h2->GetDimension() != 1) {
        THROW(JValueOutOfRange, "JTestKolmogorov_1D::test(): Given histograms must be 1-D.");
      }
 
      const int n1 = h1 -> GetNbinsX();
      const int n2 = h2 -> GetNbinsX();
 
      if(n1 != n2) {
        THROW(JValueOutOfRange, "JTestKolmogorov_1D::test(): Histograms with different bining. The objects: " <<
              h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
      }
 
      TH1* h3 = (TH1*) h1->Clone(h1->GetName() == h2->GetName() ?
                                 MAKE_CSTRING(h1->GetName()                            << "_" << testName) :
                                 MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_" << testName));
      
      const double s1 = 1./h1->Integral();
      const 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;
      }
 
      bool afunc1 = false;
      bool afunc2 = false;
    
      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) {
        THROW(JValueOutOfRange, "JTestKolmogorov_1D::test(): Errors are zero for both histograms.");
      }
 
      double dmax = 0, c1 = 0, c2 = 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);
        
        dmax = TMath::Max(dmax,TMath::Abs(c1-c2));
        
        h3->SetBinContent(bin,p);
      }
 
      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));
      }
    
      const double pValue = TMath::KolmogorovProb(z);
 
      const bool passed = (pValue > threshold);
 
      const JResultTitle title(testName, resultType, passed , pValue);
      
      h3->SetTitle(title.getTitle().c_str());
      h3->GetYaxis()->SetTitle(resultType.c_str());      
    
      const JTestResult r (testName,
                           JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
                           JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
                           resultType, pValue, threshold, h3, passed);
  
      this->push_back(r);
    }
 
    
    /**
     * Read test parameters from input.
     *
     * \param  in            input stream
     * \return               input stream
     */
    std::istream& read(std::istream& in) override
    {
      using namespace JPP;
      
      in >> threshold;
 
      if (threshold < 0.0 || threshold > 1.0) {
        THROW(JValueOutOfRange, "JTestKolmogorov_1D::read(): Invalid threshold value " << threshold);
      }
      
      return in;
    }
  
  private:
    
    double threshold;  //!< threshold p-value to decide if test is passed.
  };
}
  
#endif