RooHistError.cxx

Go to the documentation of this file.

/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id$
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/
 
/**
\file RooHistError.cxx
\class RooHistError
\ingroup Roofitcore
 
Singleton class used to calculate the error bars
for each bin of a RooHist object. Errors are calculated by integrating
a specified area of a Poisson or Binomail error distribution.
**/
 
#include "RooHistError.h"
#include "RooBrentRootFinder.h"
#include "RooMsgService.h"
 
#include "Math/QuantFuncMathCore.h" // ROOT::Math::chisquared_quantile, chisquared_quantile_c
 
#include "Riostream.h"
 
#include <cassert>
#include <cmath>
 
using std::endl;
 
 
////////////////////////////////////////////////////////////////////////////////
/// Return a reference to a singleton object that is created the
/// first time this method is called. Only one object will be
/// constructed per ROOT session.
 
const RooHistError &RooHistError::instance()
{
  static RooHistError _theInstance;
  return _theInstance;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Calculate a confidence interval for the expected number of events given n
/// observed (unweighted) events. The interval will contain the same probability
/// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
/// the point estimate n (ie, for small n) in which case the interval is adjusted
/// to start at n.
 
bool RooHistError::getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma) const
{
   // sanity checks
   if (n < 0) {
      oocoutE(nullptr, Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n
                                 << std::endl;
      return false;
   }
   if (!(nSigma > 0.0)) {
      oocoutE(nullptr, Plotting) << "RooHistError::getPoissonInterval: nSigma must be > 0, got " << nSigma << std::endl;
      return false;
   }
 
   // Convert "number of sigmas" to central Gaussian probability beta, and
   // corresponding two-sided tail probability alpha.
   const double beta = std::erf(nSigma / std::sqrt(2.0));
   const double alpha = 1.0 - beta;
 
   // Special case n = 0 (boundary at mu >= 0).
   // Use a one-sided interval including the MLE mu=0:
   //   mu1 = 0
   //   P(N <= 0 | mu2) = alpha,
   // with alpha = 1 - erf(nSigma / sqrt(2)).
   if (n == 0) {
      mu1 = 0.0;
      mu2 = 0.5 * ROOT::Math::chisquared_quantile(1.0 - alpha, 2.0 * (n + 1));
      return true;
   }
 
   // For n>0, use the central (equal-tailed) Garwood interval, which
   // corresponds to allocating alpha/2 in each tail.
   const double a2 = 0.5 * alpha;
   mu1 = 0.5 * ROOT::Math::chisquared_quantile(a2, 2.0 * n);
   mu2 = 0.5 * ROOT::Math::chisquared_quantile_c(a2, 2.0 * (n + 1));
   return true;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
/// If the return values is false and error occurred.
 
bool RooHistError::getBinomialIntervalAsym(Int_t n, Int_t m,
                    double &asym1, double &asym2, double nSigma) const
{
  // sanity checks
  if(n < 0 || m < 0) {
    oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << std::endl;
    return false;
  }
 
  // handle the special case of no events in either category
  if(n == 0 && m == 0) {
    asym1= -1;
    asym2= +1;
    return true;
  }
 
  // handle cases when n,m>100 (factorials in BinomialSum will overflow around 170)
  if ((n>100&&m>100)) {
    double N = n ;
    double M = m ;
    double asym = 1.0*(N-M)/(N+M) ;
    double approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
 
    asym1 = asym-nSigma*approxErr ;
    asym2 = asym+nSigma*approxErr ;
    return true ;
  }
 
  // swap n and m to ensure that n <= m
  bool swapped(false);
  if(n > m) {
    swapped= true;
    Int_t tmp(m);
    m= n;
    n= tmp;
  }
 
  // create the function objects to use
  bool status(false);
  BinomialSumAsym upper(n,m);
  if(n > 0) {
    BinomialSumAsym lower(n-1,m+1);
    status= getInterval(&upper,&lower,(double)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
  }
  else {
    status= getInterval(&upper,nullptr,(double)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
  }
 
  // undo the swap here
  if(swapped) {
    double tmp(asym1);
    asym1= -asym2;
    asym2= -tmp;
  }
 
  return status;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
/// If the return values is false and error occurred.
 
bool RooHistError::getBinomialIntervalEff(Int_t n, Int_t m,
                    double &asym1, double &asym2, double nSigma) const
{
  // sanity checks
  if(n < 0 || m < 0) {
    oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << std::endl;
    return false;
  }
 
  // handle the special case of no events in either category
  if(n == 0 && m == 0) {
    asym1= -1;
    asym2= +1;
    return true;
  }
 
  // handle cases when n,m>80 (factorials in BinomialSum will overflow around 170)
  if ((n>80&&m>80)) {
    double N = n ;
    double M = m ;
    double asym = 1.0*(N)/(N+M) ;
    double approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
 
    asym1 = asym-nSigma*0.5*approxErr ;
    asym2 = asym+nSigma*0.5*approxErr ;
    return true ;
  }
 
  // swap n and m to ensure that n <= m
  bool swapped(false);
  if(n > m) {
    swapped= true;
    Int_t tmp(m);
    m= n;
    n= tmp;
  }
 
  // create the function objects to use
  bool status(false);
  BinomialSumEff upper(n,m);
  double eff = (double)(n)/(n+m) ;
  if(n > 0) {
    BinomialSumEff lower(n-1,m+1);
    status= getInterval(&upper,&lower,eff,0.1,asym1,asym2,nSigma*0.5);
  }
  else {
    status= getInterval(&upper,nullptr,eff,0.1,asym1,asym2,nSigma*0.5);
  }
 
  // undo the swap here
  if(swapped) {
    double tmp(asym1);
    asym1= 1-asym2;
    asym2= 1-tmp;
  }
 
  return status;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Calculate a confidence interval using the cumulative functions provided.
/// The interval will be "central" when both cumulative functions are provided,
/// unless this would exclude the pointEstimate, in which case a one-sided interval
/// pinned at the point estimate is returned instead.
 
bool RooHistError::getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, double pointEstimate,
             double stepSize, double &lo, double &hi, double nSigma) const
{
  // sanity checks
  assert(nullptr != Qu || nullptr != Ql);
 
  // convert number of sigma into a confidence level
  double beta= std::erf(nSigma/sqrt(2.));
  double alpha= 0.5*(1-beta);
 
  // Does the central interval contain the point estimate?
  bool ok(true);
  double loProb(1);
  double hiProb(0);
  if(nullptr != Ql) loProb= (*Ql)(&pointEstimate);
  if(nullptr != Qu) hiProb= (*Qu)(&pointEstimate);
 
  if (Qu && (nullptr == Ql || loProb > alpha + beta))  {
    // Force the low edge to be at the pointEstimate
    lo= pointEstimate;
    double target= loProb - beta;
    hi= seek(*Qu,lo,+stepSize,target);
    RooBrentRootFinder uFinder(*Qu);
    ok= uFinder.findRoot(hi,hi-stepSize,hi,target);
  }
  else if(Ql && (nullptr == Qu || hiProb < alpha)) {
    // Force the high edge to be at pointEstimate
    hi= pointEstimate;
    double target= hiProb + beta;
    lo= seek(*Ql,hi,-stepSize,target);
    RooBrentRootFinder lFinder(*Ql);
    ok= lFinder.findRoot(lo,lo,lo+stepSize,target);
  }
  else if (Qu && Ql) {
    // use a central interval
    lo= seek(*Ql,pointEstimate,-stepSize,alpha+beta);
    hi= seek(*Qu,pointEstimate,+stepSize,alpha);
    RooBrentRootFinder lFinder(*Ql);
    RooBrentRootFinder uFinder(*Qu);
    ok= lFinder.findRoot(lo,lo,lo+stepSize,alpha+beta);
    ok|= uFinder.findRoot(hi,hi-stepSize,hi,alpha);
  }
  if(!ok) oocoutE(nullptr,Plotting) << "RooHistError::getInterval: failed to find root(s)" << std::endl;
 
  return ok;
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Scan f(x)-value until it changes sign. Start at the specified point and take constant
/// steps of the specified size. Give up after 1000 steps.
 
double RooHistError::seek(const RooAbsFunc &f, double startAt, double step, double value) const
{
  Int_t steps(1000);
  double min(f.getMinLimit(1));
  double max(f.getMaxLimit(1));
  double x(startAt);
  double f0 = f(&startAt) - value;
  do {
    x+= step;
  }
  while(steps-- && (f0*(f(&x)-value) >= 0) && ((x-min)*(max-x) >= 0));
  assert(0 != steps);
  if(x < min) x= min;
  if(x > max) x= max;
 
  return x;
}
 
 
 
////////////////////////////////////////////////////////////////////////////////
/// Create and return a PoissonSum function binding
 
RooAbsFunc *RooHistError::createPoissonSum(Int_t n)
{
  return new PoissonSum(n);
}
 
 
////////////////////////////////////////////////////////////////////////////////
/// Create and return a BinomialSum function binding
 
RooAbsFunc *RooHistError::createBinomialSum(Int_t n, Int_t m, bool eff)
{
  if (eff) {
    return new BinomialSumEff(n,m) ;
  } else {
    return new BinomialSumAsym(n,m) ;
  }
}