JTestPoissonLogLikelihoodRatio.hh

Go to the documentation of this file.

#ifndef __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIO__
#define __JCOMPAREHISTOGRAMS__JTESTPOISSONLOGLIKELIHOODRATIO__
 
#include <istream>
#include <ostream>
 
#include "JLang/JException.hh"
 
#include "JCompareHistograms/JResultTitle.hh"
#include "JCompareHistograms/JTest_t.hh"
 
#include "TH1.h"
 
 
/**
 * \author bjung
 */
namespace JCOMPAREHISTOGRAMS {
  
  /**
   * Implementation of the Poisson log-likelihood ratio test.
   * The first histogram is treated as an Asimov dataset corresponding to a given null hypothesis,
   * which is compared to an alternative hypothesis given by the second histogram.
   * The uncertainty on the expectation is taken into account as an additional Poissonian constraint term in the likelihood of each bin (\f$\Lambda_i\f$):
   *
   * \f[
   *   \Lambda_i = \frac{(\gamma_i m_i)^{n_i}}{n_i!} e^{-\gamma_i m_i} \cdot \frac{(\gamma_i M_i)^{M_i}}{M_i!} e^{-\gamma_i M_i},
   * \f]
   *
   * where:
   *
   * - \f$\gamma_i\f$ corresponds to a scaling parameter,
   * - \f$m_i\f$ corresponds to the expectation in bin \f$i\f$,
   * - \f$n_i\f$ corresponds to the observed number of events in bin \f$i\f$ of the Asimov dataset and
   * - \f$M_i = \frac{m_i^2}{\sigma_{m_i}^2}\f$ corresponds to the effective number of exected events in bin \f$i\f$,\n where \f$\sigma_{m_i}\f$ corresponds to the standard deviation of the expectation \f$m_i\f$.
   *
   * This is inspired on the Beeston-Barlow method, described in the following publication:
   * Computer Physics Communications, Volume 77, Issue 2, 1993
   * "Fitting using finite Monte Carlo samples",
   * Roger Barlow and Christine Beeston.
   * DOI: 10.1016/0010-4655(93)90005-W
   */
  class JTestPoissonLogLikelihoodRatio: 
    public JTest_t
  { 
  public:  
 
    /**
     * Default constructor.
     */
    JTestPoissonLogLikelihoodRatio() :
      JTest_t("Poisson_NLLR", "NLLR"),
      threshold(0.0)
    {}
 
    
    /**
     * 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(JValueOutOfRange, "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);
      }
 
      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 LnLratio = 0.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) {
            
            double contribution = 0.0;
            
            const double n1 = h1->GetBinContent(i,j,k);
            const double n2 = h2->GetBinContent(i,j,k);
              
            const double e1 = h1->GetBinError(i,j,k);
            const double e2 = h2->GetBinError(i,j,k);
 
            if (n2 > 0.0) {
 
              if (!(e2 > 0.0)) {
                THROW(JDivisionByZero, "JTestPoissonLogLikelihoodRatio::test(): Invalid effective bin content for histogram " << h2->GetName() << " at bin (" << i << ',' << j << ',' << k << ')');
              }
 
              const double N2 = n2*n2/e2/e2; // Effective bin content of alternative model
                
              if (n1 > 0.0) {
 
                if (!(e1 > 0.0)) {
                  THROW(JDivisionByZero, "JTestPoissonLogLikelihoodRatio::test(): Invalid effective bin content for histogram " << h1->GetName() << " at bin (" << i << ',' << j << ',' << k << ')');
                }
 
                const double f1 = (n1 + N2) / (n2 + N2); // Optimal scale factor for alternative model bin
 
                const double N1 = n1*n1/e1/e1; // Effective bin content of Asimov dataset
                    
                contribution = 2 * (N2               - N1             +
                                    lgamma(N2 + 1)   - lgamma(N1 + 1) +
                                    n1*log(n1/n2/f1) - N2*log(f1*N2)  +
                                    N1*log(N1));
                    
              } else {
 
                const double f1 = N2 / (n2 + N2); // Optimal scale factor for alternative model bin
 
                contribution = 2 * (N2 + lgamma(N2 + 1) - N2*log(f1*N2));
              }
            }
            
            h3->SetBinContent(i,j,k, contribution);
            
            LnLratio += contribution;
          }
        }
      }
 
      const bool passed = (LnLratio < threshold);
 
      const JResultTitle title(testName, resultType, passed, LnLratio);
 
      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, LnLratio, 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, "JTestPoissonLogLikelihoodRatio::read(): Given chi-square threshold " << threshold << " is invalid");
      }
       
      return in;
    }
  
  private:
    
    double threshold;  //!< threshold chi-square to decide if test is passed.
  };
}
  
#endif