JTestPoissonLogLikelihoodRatioBeestonBarlow.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIOBEESTONBARLOW__
#define __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIOBEESTONBARLOW__
 
#include <istream>
#include <ostream>
 
#include "JLang/JException.hh"
 
#include "JCompareHistograms/JResultTitle.hh"
#include "JCompareHistograms/JTest_t.hh"
 
#include "TH1.h"
#include "TH2.h"
 
#include "RooRealVar.h"
#include "RooRealConstant.h"
#include "RooParamHistFunc.h"
#include "RooHistConstraint.h"
#include "RooProdPdf.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooHistFunc.h"
#include "RooHistPdf.h"
#include "RooRealSumPdf.h"
#include "RooFitResult.h"
 
 
/**
 * \author bjung
 */
namespace JCOMPAREHISTOGRAMS {
  
  /**
   * Implementation of the Poisson log-likelihood ratio test.\n
   * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis,\n
   * which is compared to an alternative hypothesis given by the second histogram.\n\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
   */
  class JTestPoissonLogLikelihoodRatioBeestonBarlow: 
    public JTest_t
  { 
  public:  
 
    /**
     * Default constructor.
     */
    JTestPoissonLogLikelihoodRatioBeestonBarlow() :
      JTest_t("Poisson_NLLR_Beeston-Barlow", "NLLR")
    {}
 
    
    /**
     * Applies a Poissonian log-likelihood ratio test to the two given histograms.\n
     * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis.\n
     * The second histogram is treated as the expectation for the alternative hypothesis, to which the null hypothesis is compared.
     * 
     * \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(JCastException, "JTestPoissonLogLikelihood::test(): Could not cast given TObjects to TH1.");
      }
 
      if(h1->GetNbinsX() != h2->GetNbinsX() ||
         h1->GetNbinsY() != h2->GetNbinsY() ||
         h1->GetNbinsZ() != h2->GetNbinsZ()) {
        THROW(JValueOutOfRange, "JTestPoissonLogLikelihood::test(): Histograms with different binning. The objects: " <<
              h1->GetName() << " and " << h2->GetName() << " can not be compared." << endl);
      }
 
      // Flatten histograms
      
      const Int_t N = h1->GetNbinsX() * h1->GetNbinsY() * h1->GetNbinsZ();
 
      TH1I H1("H1", NULL, N, 0.5, N+0.5);
      TH1I H2("H2", NULL, N, 0.5, N+0.5);
 
      for (Int_t i=1 ; i <= h1->GetNbinsX() ; ++i) {
        for (Int_t j=1 ; j <= h1->GetNbinsY() ; ++j) {
          for (Int_t k=1 ; k <= h1->GetNbinsZ() ; ++k) {
 
            const Int_t n = h1->GetBin(i,j,k);
 
            const double y1 = h1->GetBinContent(i,j,k);
            const double y2 = h2->GetBinContent(i,j,k);
 
            const double e1 = h1->GetBinError(i,j,k);
            const double e2 = h2->GetBinError(i,j,k);
 
            H1.SetBinContent(n, y1);
            H2.SetBinContent(n, y2);
 
            H1.SetBinError(n, e1);
            H2.SetBinError(n, e2);
          }
        }
      }
 
      RooRealVar index("index", "index", 0.5, N+0.5);
 
      RooDataHist dh1("dh1", h1->GetTitle(), index, &H1);
      RooDataHist dh2("dh2", h2->GetTitle(), index, &H2);
 
      // Compute NLL of Asimov data, given the altnerative model with Beeston-Barlow correction      
      
#if ROOT_VERSION_CODE < ROOT_VERSION(6,32,0)
      RooParamHistFunc ph2 ("ph2", "ph2", dh2);
#else
      RooParamHistFunc ph2 ("ph2", "ph2", dh2, index);
#endif
      RooHistConstraint hc2("hc2", "hc2", ph2);
 
      RooRealSumPdf model_tmp2("model_tmp2", "model_tmp2", ph2, RooRealConstant::value(1.0));
      
      RooProdPdf model2("model2", "model2", hc2, RooFit::Conditional(model_tmp2,index));
 
      RooFitResult* result2 = (RooFitResult*) model2.fitTo(dh1,
                                                           RooFit::Offset("bin"),
                                                           RooFit::SumW2Error(false),
                                                           RooFit::PrintLevel(-1),
                                                           RooFit::Verbose(false),
                                                           RooFit::Save(true));
      
      if (result2 == NULL) {
        THROW(JNullPointerException, "JTestPoissonLogLikelihoodRatioBeestonBarlow::test(): Unable to retrieve Beeston-Barlow fit results");
      }
      
      const double nllratio = 2 * (result2->minNll() + hc2.getLogVal());
      const bool   passed   = (nllratio < threshold);
      
      const JResultTitle title(testName, resultType, passed, nllratio);
      
      const JTestResult r (testName,
                           JRootObjectID(MAKE_STRING(h1->GetDirectory()->GetPath() << h1->GetName())),
                           JRootObjectID(MAKE_STRING(h2->GetDirectory()->GetPath() << h1->GetName())),
                           resultType, nllratio, threshold, result2, 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, "JTestPoissonLogLikelihoodRatioBeestonBarlow::read(): Given chi-square threshold " << threshold << " is invalid");
      }
      
      return in;
    }
  
  private:
    
    double threshold;  //!< threshold chi-square to decide if test is passed.
  };
}
  
#endif