JTestRuns_Slice2D.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTRUNS_2D__
#define __JCOMPAREHISTOGRAMS__JTESTRUNS_2D__

#include <istream>
#include <ostream>

#include "JLang/JException.hh"

#include "JCompareHistograms/JResultTitle.hh"
#include "JCompareHistograms/JTest_t.hh"

#include "TH1.h"
#include "TH2.h"


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

namespace JCOMPAREHISTOGRAMS {

  /**
   * Implementation of the 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.\n
   * The parameter `slice` can have the values x, X, y or Y. The histograms are sliced along the corresponding axis, and the runs test is applied to each slice.\n 
   * This test uses the input parameter `threshold` to evaluate whether the test is passed or failed for each slice.\n 
   * The evaluation is done by comparing the threshold value with the result of the runs test.\n 
   * The results for each slice are stored in the results buffer (see JTest_t())
   */
  class JTestRuns_Slice2D: 
    public JTest_t
  { 
  public:  

    /**
     * Default constructor.
     */
    JTestRuns_Slice2D() :
      JTest_t("Runs_2D", "#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 TH2* h1 = dynamic_cast<const TH2*>(o1);
      const TH2* h2 = dynamic_cast<const TH2*>(o2);

      if (h1 == NULL || h2 == NULL) {
        THROW(JValueOutOfRange, "JTestKolmogorov_Slice2D::test(): Could not cast given TObjects to TH2.");
      }

      TH1* h3 = NULL;

      const char* const h3name = (h1->GetName() == h2->GetName() ?
                                  MAKE_CSTRING(h1->GetName()                            << "_" << testName << "_" << slice) :
                                  MAKE_CSTRING(h1->GetName() << "_VS_" << h2->GetName() << "_" << testName << "_" << slice));
      
      int  nFailures = 0;
      int  nSlices   = 0;
    
      if (slice == 'x' || slice == 'X') {
    
        const int nSlices1 = h1->GetNbinsX();
        const int nSlices2 = h2->GetNbinsX();
        
        if (nSlices1 != nSlices2) {
          THROW(JValueOutOfRange, "JTestRuns_Slice2D::test(): Histograms with different binning. The objects: " << h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
        }

        nSlices = nSlices1;

        h3 = h1->ProjectionX(h3name);

        for (int i=1 ; i<=nSlices1 ; ++i){

          const std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i));
             
          const TH1* s1 = h1->ProjectionY (sliceName.c_str(),i,i);
          const TH1* s2 = h2->ProjectionY (sliceName.c_str(),i,i);

          if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; }
                
          const double R = s1->Integral();
          const double T = s2->Integral();

          bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0;
          
          int    n = 1;
          double p = (a ? 1 : 0);
          double q = (a ? 0 : 1);

          for (int i = 2 ; i<s1->GetNbinsX() ; ++i){
      
            const bool b = ((T/R)*s1->GetBinContent(i) - s2->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;

          if (fabs(d) > threshold) nFailures++;
        
          h3->SetBinContent(i,fabs(d));
        }
            
      } else if (slice == 'y' || slice == 'Y') {
      
        const int nSlices1 = h1->GetNbinsX();
        const int nSlices2 = h2->GetNbinsX();
      
        if (nSlices1 != nSlices2) {
          THROW(JValueOutOfRange, "JTestRuns_Slice2D::test(): Histograms with different binning. The objects: " << h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
        }

        nSlices = nSlices1;

        h3 = h1->ProjectionY(h3name);

        for (int i=1 ; i<=nSlices1 ; ++i){

          const std::string sliceName = MAKE_STRING(h3->GetName() << "_" << to_string(i));
             
          const TH1* s1 = h1->ProjectionX (sliceName.c_str(),i,i);
          const TH1* s2 = h2->ProjectionX (sliceName.c_str(),i,i);

          if (!(s1->GetSumOfWeights() > 0 && s2->GetSumOfWeights() > 0)) { continue; }
                
          const double R = s1->Integral();
          const double T = s2->Integral();

          bool a = ((T/R)*s1->GetBinContent(1) - s2->GetBinContent(1)) < 0;
          
          int    n = 1;
          double p = (a ? 1 : 0);
          double q = (a ? 0 : 1);

          for (int i = 2 ; i<s1->GetNbinsX() ; ++i){
      
            const bool b = ((T/R)*s1->GetBinContent(i) - s2->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;

          if (fabs(d) > threshold) nFailures++;
        
          h3->SetBinContent(i,fabs(d));
        }
      }

      const bool passed = (nFailures/nSlices < failuresThreshold);
      
      const JResultTitle title(testName, resultType, passed , nFailures);

      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, nFailures, failuresThreshold, 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 >> failuresThreshold >> slice;

      if (threshold < 0.0) {
        THROW(JValueOutOfRange, "JTestRuns_Slice2D::read(): Invalid failuresThreshold value " << failuresThreshold);
      }
      
      if (failuresThreshold < 0.0 || failuresThreshold > 1.0) {
        THROW(JValueOutOfRange, "JTestRuns_Slice2D::read(): Invalid failuresThreshold value " << failuresThreshold);
      }

      if (slice != 'x' && slice != 'X' && slice != 'y' && slice != 'Y') {
        THROW(JValueOutOfRange, "JTestRuns_Slice2D::read(): Invalid slice option \'" << slice << "\'");
      }
                                                                           
      return in;
    };

  private:
    
    double threshold;                       //!< threshold value to decide if test is passed.
    double failuresThreshold;               //!< threshold value to decide if test is passed.
    char   slice;                           //!< Axis to slice. x or X for x-axis, y or Y for y-axis, n or N for None.
  };
}
  
#endif