JTestChi2.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTCHI2__
#define __JCOMPAREHISTOGRAMS__JTESTCHI2__

#include <istream>
#include <ostream>

#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 the Chi2 test for ROOT 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
   * 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 produced by JChi2Test(). The output of a Chi2 test is a p-value.\n 
   * The parameter threshold() should therefore be a real value between 0 and 1.
   */
  class JTestChi2: 
    public JTest_t
  { 
  public:  

    /**
     * Default constructor.
     */
    JTestChi2() :
      JTest_t("Chi2", "p-Value(chi2)")
    {}

    
    /**
     * Applies Chi2 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 == NULL || h2 == NULL) {
        THROW(JValueOutOfRange, "JTestChi2::test(): Could not cast given TObjects to TH1.");
      }

      vector<Double_t> residuals(h1->GetNcells());
      
      const double chi2 = h1->Chi2Test(h2 , options.c_str(), &residuals[0]);

      TH1* h3 = (TH1*) h1->Clone(h1->GetName() == h2->GetName() ?
                                 MAKE_CSTRING(h1->GetName()                            << "_" << testName) :
                                 MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_" << testName));

      h3->Reset();
      
      for (int i = 0; i < h1->GetNcells(); ++i) {
        h3->SetBinContent(i+1, residuals[i]);
      }

      const bool passed = ( options.find("CHI2") != string::npos ?
                            (chi2 < threshold) :
                            (chi2 > threshold) );

      const JResultTitle title(testName, resultType, passed, chi2);
      
      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, chi2, 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 >> options;

      if (threshold < 0.0) {
        THROW(JValueOutOfRange, "JTestChi2::read(): Invalid threshold value " << threshold);    
      }

      return in;
    }    
  
  private:
    
    double      threshold;  //!< Threshold p-value or chi2-value to decide if test is passed.
    std::string options;    //!< Options for the ROOT chi2 test.
  };
}
  
#endif