JTestZero.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTZERO__
#define __JCOMPAREHISTOGRAMS__JTESTZERO__

#include <istream>
#include <ostream>

#include "Math/Math.h"
#include "Math/Error.h"
#include "Math/ProbFuncMathCore.h"
#include "Math/SpecFuncMathCore.h"

#include "JLang/JException.hh"

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

#include "JTools/JQuantile.hh"
#include "JTools/JRange.hh"

#include "TObject.h"
#include "TMath.h"
#include "TH1.h"


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

namespace JCOMPAREHISTOGRAMS {

  /**
   * Implementation of a bin-by-bin compatibility test for histograms with low bin contents.\n
   * The test loops over all the bins of both histograms and performs the following operations.\n 
   * -Calculates the probability that the observed bin content for histogram A is obtained from a Poisson of parameter 1.\n
   * -Compares the previous result with the `threshold` value given as an input parameter. The result is true if the probability is higher than the `threshold`.\n
   * -Calculates the probability that the observed bin content for histogram B is obtained from a Poisson of parameter 1.\n
   * -Compares the previous result with the `threshold` value given as an input parameter. The result is true if the probability is higher than the `threshold`.\n
   * -Compares the results from both bins: If both are true, or both are false, the test is passed. If one is true and the other is false, the test is failed.\n
   * At the end of the loop, a failure fraction is computed and compared to the `outliersThreshold` parameter. If the failure fraction is above the `threshold`, the test is failed.\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 JTestZero: 
    public JTest_t
  { 
  public:  

    /**
     * Default constructor.
     */
    JTestZero() :
      JTest_t("Zero", "failure_fraction")
    {}

      
    /**
     * Bin-by-bin comparison for ROOT histograms, of compatibility with a Poisson pdf of parameter 1.
     * 
     * \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 TH2.");
      }
    
      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();

      if(nx1 != nx2 || ny1 != ny2 || nz1 != nz2) {
        THROW(JValueOutOfRange, "JTestZero::test(): Histograms with different binning. The objects: " <<
              h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
      }

      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 failures = 0;
    
      for (int i=1 ; i<nx1+1 ; ++i) {
        for (int j=1 ; j<ny1+1 ; ++j) {
          for (int k=1 ; k<nz1+1 ; ++k) {

            const int Nbin = h3->GetBin(i,j,k);
        
            const double m = h1->GetBinContent(i,j,k);
            const double n = h2->GetBinContent(i,j,k);
        
            const double p1 = 1 - ROOT::Math::poisson_cdf(m,1);
            const double p2 = 1 - ROOT::Math::poisson_cdf(n,1);         

            if ((p1 > threshold && p2 < threshold) || 
                (p1 < threshold && p2 > threshold)) {
              
              failures+=1./(nx1*ny1);
              h3->Fill(Nbin);
            }
          }
        }
      }
      
      const bool passed = (failures < outliersThreshold);

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

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

      if (outliersThreshold < 0.0 || outliersThreshold > 1.0) {
        THROW(JValueOutOfRange, "JTestZero::read(): Invalid outliersThreshold value " << outliersThreshold);
      }

      return in;
    }
      
  private:
    
    double threshold;                                  //!< threshold p-value to decide if test is passed.
    double outliersThreshold;                          //!< Fraction of bins allowed to fail.
  };
}
  
#endif