JTestChi2_Bin.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTCHI2_BIN__
#define __JCOMPAREHISTOGRAMS__JTESTCHI2_BIN__
 
#include <istream>
#include <ostream>
 
#include "JLang/JException.hh"
 
#include "JCompareHistograms/JTest_t.hh"
#include "JCompareHistograms/JResultTitle.hh"
 
#include "TH1.h"
 
 
/**
 * \author rgruiz, bjung
 */
namespace JCOMPAREHISTOGRAMS {}
namespace JPP { using namespace JCOMPAREHISTOGRAMS; }
 
namespace JCOMPAREHISTOGRAMS {
 
  /**
   * Implementation of a bin-by-bin Chi2 test for 2D 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.
   */
  class JTestChi2_Bin: 
    public JTest_t
  { 
  public:
    
    /**
     * Default constructor.
     */
    JTestChi2_Bin() :
      JTest_t("Chi2_Bin", "Outliers[%]")
    {}
 
    
    /**
     * Tests the statistical compatibility of two ROOT 2D 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 == NULL || h2 == NULL) {
        THROW(JValueOutOfRange, "JTestChi2_Bin::test(): Could not cast given TObjects to TH1.");
      }
    
      const int nx1 = h1->GetNbinsX();
      const int nx2 = h2->GetNbinsX();      
      const int ny1 = h1->GetNbinsY();
      const int ny2 = h2->GetNbinsY();
      const int nz1 = h1->GetNbinsZ();
      const int nz2 = h2->GetNbinsZ();
      
      const double M = h1->GetSumOfWeights();
      const double N = h2->GetSumOfWeights();
 
      if(nx1 != nx2 || ny1 != ny2 || nz1 != nz2 || M == 0 || N == 0) {
        THROW(JValueOutOfRange, "JTestChi2_t::test(): Histograms with different binning or with zero entries.");
      }
 
      TH1* h3 = (TH1*) h1->Clone(h1->GetName() == h2->GetName() ?
                                 MAKE_CSTRING(h1->GetName()                            << "_" << testName) :
                                 MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_" << testName));
 
      h3->Reset();
      
      double outliers = 0;
    
      for (int i=1 ; i<nx1 ; ++i) {
        for (int j=1 ; j<ny1 ; ++j) {
          for (int k=1 ; k<nz1 ; ++k) {
              
            const double m = h1->GetBinContent(i,j,k);
            const double n = h2->GetBinContent(i,j,k);
 
            const double em = h1->GetBinError(i,j,k);
            const double en = h2->GetBinError(i,j,k);
 
            if (n > 0.0 && m > 0.0 && em > 0.0 && en > 0.0) {
 
              const double delta = N*m - M*n;
              const double sigma = M*M*en*en + N*N*em*em;
              const double chi2  = delta*delta / sigma;
 
              if (fabs(chi2) > chi2Threshold) {
                outliers += 1./(nx1*ny1);
              }
            
              h3->SetBinContent(i,j,k,chi2);
            }
          }
        }
      }
 
      const bool passed = (outliers < outliersThreshold);
 
      const JResultTitle title(testName, resultType, passed, 100.*outliers);
      
      h3->SetTitle(title.getTitle().c_str());
 
      const int Ndims = h3->GetDimension();
 
      if        (Ndims == 1) {
        h3->GetYaxis()->SetTitle(resultType.c_str());
      } else if (Ndims == 2) {
        h3->GetZaxis()->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, 100.*outliers, 100.*outliersThreshold, 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 >> chi2Threshold >> outliersThreshold;
 
      if (chi2Threshold < 0.0) {
        THROW(JValueOutOfRange, "JTestChi2_Bin::read(): Invalid chi2-threshold value " << chi2Threshold);       
      }
 
      if (outliersThreshold < 0.0 || outliersThreshold > 1.0) {
        THROW(JValueOutOfRange, "JTestChi2_Bin::read(): Invalid outliers threshold value " << outliersThreshold);       
      }
 
      return in;
    }
 
    
  private:
    
    double chi2Threshold;                        //!< threshold chi2-value to decide if test is passed for a bin.
    double outliersThreshold;                    //!< fraction of bins allowed to fail the test.
  };
}
  
#endif