JTestRuns_1D.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTRUNS_1D__
#define __JCOMPAREHISTOGRAMS__JTESTRUNS_1D__

#include <istream>
#include <ostream>

#include "JLang/JException.hh"

#include "JCompareHistograms/JTest_t.hh"

#include "TH1.h"


/**
 *  \author rgruiz, bjung
 */
namespace JCOMPAREHISTOGRAMS {}
namespace JPP { using namespace JCOMPAREHISTOGRAMS; }

namespace JCOMPAREHISTOGRAMS {

  /**
   * Implementation of the Runs 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 runs test. 
   */
  class JTestRuns_1D: 
    public JTest_t
  { 
  public:  

    /**
     * Default constructor.
     */
    JTestRuns_1D() :
      JTest_t("Runs_1D", "#sigma")
    {}

    
    /**
     * Tests the statistical compatibility of two ROOT TObjects
     * 
     * \param o1 First object
     * \param o2 Second object
     */
    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, "JTestRuns_1D::test(): Could not cast given TObjects to TH1.");
      }

      if (h1->GetDimension() != 1 || h2->GetDimension() != 1) {
        THROW(JValueOutOfRange, "JTestRuns_1D::test(): Given histograms must be 1-D.");
      }

      if(h1->GetNbinsX() != h2->GetNbinsX()) {
        THROW(JValueOutOfRange, "JTestIdentical::test(): Histograms with different bining. The objects: " <<
              h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
      }

      const double R = h1->Integral();
      const double T = h2->Integral();

      TH1* h3 = (TH1*) h1->Clone(h1->GetName() == h2->GetName() ?
                                 MAKE_CSTRING(h1->GetName()                            << "_" << testName) :
                                 MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_" << testName));
    
      for (int i=1 ; i<h1->GetNbinsX() ; ++i) {
        h3->SetBinContent(i, (T/R)*h1->GetBinContent(i) - h2->GetBinContent(i));
      }
        
      bool a = ((T/R)*h1->GetBinContent(1) - h2->GetBinContent(1)) < 0;

      int    n = 1;
      double p = (a ? 1 : 0);
      double q = (a ? 0 : 1);

      for (int i = 2 ; i<h1->GetNbinsX() ; ++i){
      
        const bool b = ((T/R)*h1->GetBinContent(i) - h2->GetBinContent(i)) < 0;

        (b ? ++p : ++q);

        if (b != a){
          ++n;
          a=b;
        }
      }

      const double N = 1 + 2*p*q/(p+q) ;
      const double s = sqrt( 2*p*q*(2*p*q-p-q)/(p+q)/(p+q)/(p+q-1) );
      const double d = (n-N)/s;
    
      const bool passed = (fabs(d) < threshold);

      const JResultTitle title(testName, resultType , passed , fabs(d));
      
      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, fabs(d), 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) {
        THROW(JValueOutOfRange, "JTestRuns_1D::read(): Invalid threshold value " << threshold);
      }
      
      return in;
    }
    
  private:
    
    double threshold;  //!< threshold value to decide if test is passed.
  };
}
  
#endif