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/JResultTitle.hh"
#include "JCompareHistograms/JTest_t.hh"
 
#include "TH1.h"
 
 
/**
 *  \author rgruiz, bjung
 */
namespace JCOMPAREHISTOGRAMS {}
namespace JPP { using namespace JCOMPAREHISTOGRAMS; }
 
namespace JCOMPAREHISTOGRAMS {
 
  /**
   * Implementation of the Wald-Wolfowitz 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;
      int p = (a ? 1 : 0);
      int 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 int N = 1 + 2*p*q/(p+q);
      
      const double S1 = 2*p*q*(2*p*q-p-q);
      const double S2 = (p+q)*(p+q)*(p+q-1);
      const double s  = sqrt(S1/S2);
      
      const double d = (N == n ? 0.0 : (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