JTestSignificance.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTSIGNIFICANCE__
#define __JCOMPAREHISTOGRAMS__JTESTSIGNIFICANCE__

#include <istream>
#include <ostream>

#include "JLang/JException.hh"

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

#include "TH1.h"


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

namespace JCOMPAREHISTOGRAMS {

  /**
   * Significance test
   *
   * Compares two histograms \f$ H_a \f$ and \f$ H_b \f$ by calculating the average bin significance \n
   *
   * \f[
   * S = \frac{1}{N} \cdot \sum_{i=0}^{N} \frac{|a_i - K \cdot b_i|}{\sqrt{\sigma ^2(a_i) + K^2 \cdot \sigma ^2(b_i)}}
   * \f]
   *
   * where:
   *  - \f$ N \f$ is the number of bins in \f$ H_a \f$ and \f$ H_b \f$.
   *  - \f$ a_i \f$ and \f$ b_i \f$ are the contents of bin \f$ i \f$ in \f$ H_a \f$ and \f$ H_b \f$.
   *  - \f$ K \f$ is a normalisation parameter which can be specified by the user. 
   *    - It should be interpreted as the integral of \f$ H_a \f$ divided by the integral of \f$ H_b \f$.
   *    - Note that if the specified value of K is negative, the value will be reset to the integral of \f$ H_a \f$ divided by the integral of \f$ H_b \f$.
   *  - \f$ \sigma \f$  is calculated assuming that the bin contents follow Poisson distributions.
   *
   * A threshold value for \f$ S \f$ is used to tag the test as passed or failed.
   */
  class JTestSignificance: 
    public JTest_t
  { 
  public:  

    /**
     * Default constructor.
     */
    JTestSignificance() :
      JTest_t("Significance", "Significance")
    {}

    
    /**
     * Applies Significance 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, "JTestSignificance::test(): Could not cast given TObjects to TH1.");
      }

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

      if (K < 0) {
        K = h2->GetEntries() / h1->GetEntries();
      }

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

      h3->Add(h2,-1*K);
    
      double S = 0;

      for (int i=1 ; i <= h1->GetNbinsX() ; ++i){
        for (int j=1 ; j <= h1->GetNbinsY() ; ++j){
          for (int k=1 ; k <= h1->GetNbinsZ() ; ++k){
          
            const double a = h1->GetBinContent(i,j,k);
            const double b = h2->GetBinContent(i,j,k);

            if (a!=0 || b!=0) {  
              S += fabs((a - K*b)/sqrt(a + K*K*b));
            }
          }
        }
      }

      S /= (h1->GetNbinsX()*h1->GetNbinsY()*h1->GetNbinsZ());

      const bool passed = (S > threshold);

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

      if (threshold < 0.0) {
        THROW(JValueOutOfRange, "JTestSignificance::read(): Invalid threshold value " << threshold);
      }
      
      return in;
    }
  
  private:
    
    double threshold;  //!< threshold p-value to decide if test is passed.
    double K;          //!< normalization factor between histograms.
  };
}

#endif