ROOT logo
/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id: RooHistError.cxx 25184 2008-08-20 13:59:55Z wouter $
 * 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)             *
 *****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// RooHistError is a 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.
// END_HTML
//

#include "RooFit.h"

#include "RooHistError.h"
#include "RooBrentRootFinder.h"
#include "RooMsgService.h"
#include "TMath.h"

#include "Riostream.h"
#include <assert.h>

ClassImp(RooHistError)
  ;



//_____________________________________________________________________________
const RooHistError &RooHistError::instance() 
{
  // 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.

  static RooHistError _theInstance;
  return _theInstance;
}


//_____________________________________________________________________________
RooHistError::RooHistError() 
{
  // Construct our singleton object.

  // Initialize lookup table ;
  Int_t i ;
  for (i=0 ; i<1000 ; i++) {
    getPoissonIntervalCalc(i,_poissonLoLUT[i],_poissonHiLUT[i],1.) ;
  }

}



//_____________________________________________________________________________
Bool_t RooHistError::getPoissonInterval(Int_t n, Double_t &mu1, Double_t &mu2, Double_t nSigma) const
{
  // Return 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. This method uses a lookup table to return precalculated results
  // for n<1000

  // Use lookup table for most common cases
  if (n<1000 && nSigma==1.) {
    mu1=_poissonLoLUT[n] ;
    mu2=_poissonHiLUT[n] ;
    return kTRUE ;
  }

  // Forward to calculation method 
  Bool_t ret =  getPoissonIntervalCalc(n,mu1,mu2,nSigma) ;
  return ret ;
}



//_____________________________________________________________________________
Bool_t RooHistError::getPoissonIntervalCalc(Int_t n, Double_t &mu1, Double_t &mu2, Double_t nSigma) const
{
  // 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.

  // sanity checks
  if(n < 0) {
    oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n << endl;
    return kFALSE;
  }

  // use assymptotic error if possible
  if(n > 100) {
    mu1= n - sqrt(n+0.25) + 0.5;
    mu2= n + sqrt(n+0.25) + 0.5;
    return kTRUE;
  }

  // create a function object to use
  PoissonSum upper(n);
  if(n > 0) {
    PoissonSum lower(n-1);
    return getInterval(&upper,&lower,(Double_t)n,1.0,mu1,mu2,nSigma);
  }

  // Backup solution for negative numbers
  return getInterval(&upper,0,(Double_t)n,1.0,mu1,mu2,nSigma);
}


//_____________________________________________________________________________
Bool_t RooHistError::getBinomialIntervalAsym(Int_t n, Int_t m,
					     Double_t &asym1, Double_t &asym2, Double_t nSigma) const
{
  // Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
  // If the return values is kFALSE and error occurred.

  // sanity checks
  if(n < 0 || m < 0) {
    oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
    return kFALSE;
  }

  // handle the special case of no events in either category
  if(n == 0 && m == 0) {
    asym1= -1;
    asym2= +1;
    return kTRUE;
  }

  // handle cases when n,m>100 (factorials in BinomialSum will overflow around 170)
  if ((n>100&&m>100)) {
    Double_t N = n ;
    Double_t M = m ;
    Double_t asym = 1.0*(N-M)/(N+M) ;
    Double_t approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;

    asym1 = asym-nSigma*approxErr ;
    asym2 = asym+nSigma*approxErr ;
    return kTRUE ;
  }

  // swap n and m to ensure that n <= m
  Bool_t swapped(kFALSE);
  if(n > m) {
    swapped= kTRUE;
    Int_t tmp(m);
    m= n;
    n= tmp;
  }

  // create the function objects to use
  Bool_t status(kFALSE);
  BinomialSumAsym upper(n,m);
  if(n > 0) {
    BinomialSumAsym lower(n-1,m+1);
    status= getInterval(&upper,&lower,(Double_t)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
  }
  else {
    status= getInterval(&upper,0,(Double_t)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
  }

  // undo the swap here
  if(swapped) {
    Double_t tmp(asym1);
    asym1= -asym2;
    asym2= -tmp;
  }

  return status;
}


//_____________________________________________________________________________
Bool_t RooHistError::getBinomialIntervalEff(Int_t n, Int_t m,
					     Double_t &asym1, Double_t &asym2, Double_t nSigma) const
{
  // Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
  // If the return values is kFALSE and error occurred.

  // sanity checks
  if(n < 0 || m < 0) {
    oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
    return kFALSE;
  }

  // handle the special case of no events in either category
  if(n == 0 && m == 0) {
    asym1= -1;
    asym2= +1;
    return kTRUE;
  }

  // handle cases when n,m>80 (factorials in BinomialSum will overflow around 170)
  if ((n>80&&m>80)) {
    Double_t N = n ;
    Double_t M = m ;
    Double_t asym = 1.0*(N)/(N+M) ;
    Double_t approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;

    asym1 = asym-nSigma*approxErr ;
    asym2 = asym+nSigma*approxErr ;
    return kTRUE ;
  }

  // swap n and m to ensure that n <= m
  Bool_t swapped(kFALSE);
  if(n > m) {
    swapped= kTRUE;
    Int_t tmp(m);
    m= n;
    n= tmp;
  }

  // create the function objects to use
  Bool_t status(kFALSE);
  BinomialSumEff upper(n,m);
  Double_t eff = (Double_t)(n)/(n+m) ;
  if(n > 0) {
    BinomialSumEff lower(n-1,m+1);
    status= getInterval(&upper,&lower,eff,0.1,asym1,asym2,nSigma);
  }
  else {
    status= getInterval(&upper,0,eff,0.1,asym1,asym2,nSigma);
  }

  // undo the swap here
  if(swapped) {
    Double_t tmp(asym1);
    asym1= 1-asym2;
    asym2= 1-tmp;
  }

  return status;
}



//_____________________________________________________________________________
Bool_t RooHistError::getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, Double_t pointEstimate,
				 Double_t stepSize, Double_t &lo, Double_t &hi, Double_t nSigma) const
{
  // 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.

  // sanity checks
  assert(0 != Qu || 0 != Ql);

  // convert number of sigma into a confidence level
  Double_t beta= TMath::Erf(nSigma/sqrt(2.));
  Double_t alpha= 0.5*(1-beta);

  // Does the central interval contain the point estimate?
  Bool_t ok(kTRUE);
  Double_t loProb(1),hiProb(0);
  if(0 != Ql) loProb= (*Ql)(&pointEstimate);
  if(0 != Qu) hiProb= (*Qu)(&pointEstimate);

  if(0 == Ql || loProb > alpha + beta)  {
    // Force the low edge to be at the pointEstimate
    lo= pointEstimate;
    Double_t target= loProb - beta;
    hi= seek(*Qu,lo,+stepSize,target);
    RooBrentRootFinder uFinder(*Qu);
    ok= uFinder.findRoot(hi,hi-stepSize,hi,target);
  }
  else if(0 == Qu || hiProb < alpha) {
    // Force the high edge to be at pointEstimate
    hi= pointEstimate;
    Double_t target= hiProb + beta;
    lo= seek(*Ql,hi,-stepSize,target);
    RooBrentRootFinder lFinder(*Ql);
    ok= lFinder.findRoot(lo,lo,lo+stepSize,target);
  }
  else {
    // use a central interval
    lo= seek(*Ql,pointEstimate,-stepSize,alpha+beta);
    hi= seek(*Qu,pointEstimate,+stepSize,alpha);
    RooBrentRootFinder lFinder(*Ql),uFinder(*Qu);
    ok= lFinder.findRoot(lo,lo,lo+stepSize,alpha+beta);
    ok|= uFinder.findRoot(hi,hi-stepSize,hi,alpha);
  }
  if(!ok) oocoutE((TObject*)0,Plotting) << "RooHistError::getInterval: failed to find root(s)" << endl;

  return ok;
}


//_____________________________________________________________________________
Double_t RooHistError::seek(const RooAbsFunc &f, Double_t startAt, Double_t step, Double_t value) const 
{
  // 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.
  
  Int_t steps(1000);
  Double_t min(f.getMinLimit(1)),max(f.getMaxLimit(1));
  Double_t x(startAt), 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;
}



//_____________________________________________________________________________
RooAbsFunc *RooHistError::createPoissonSum(Int_t n) 
{ 
  // Create and return a PoissonSum function binding

  return new PoissonSum(n); 
}


//_____________________________________________________________________________
RooAbsFunc *RooHistError::createBinomialSum(Int_t n, Int_t m, Bool_t eff) 
{ 
  // Create and return a BinomialSum function binding
  if (eff) {
    return new BinomialSumEff(n,m) ;
  } else {
    return new BinomialSumAsym(n,m) ;
  }
}
 RooHistError.cxx:1
 RooHistError.cxx:2
 RooHistError.cxx:3
 RooHistError.cxx:4
 RooHistError.cxx:5
 RooHistError.cxx:6
 RooHistError.cxx:7
 RooHistError.cxx:8
 RooHistError.cxx:9
 RooHistError.cxx:10
 RooHistError.cxx:11
 RooHistError.cxx:12
 RooHistError.cxx:13
 RooHistError.cxx:14
 RooHistError.cxx:15
 RooHistError.cxx:16
 RooHistError.cxx:17
 RooHistError.cxx:18
 RooHistError.cxx:19
 RooHistError.cxx:20
 RooHistError.cxx:21
 RooHistError.cxx:22
 RooHistError.cxx:23
 RooHistError.cxx:24
 RooHistError.cxx:25
 RooHistError.cxx:26
 RooHistError.cxx:27
 RooHistError.cxx:28
 RooHistError.cxx:29
 RooHistError.cxx:30
 RooHistError.cxx:31
 RooHistError.cxx:32
 RooHistError.cxx:33
 RooHistError.cxx:34
 RooHistError.cxx:35
 RooHistError.cxx:36
 RooHistError.cxx:37
 RooHistError.cxx:38
 RooHistError.cxx:39
 RooHistError.cxx:40
 RooHistError.cxx:41
 RooHistError.cxx:42
 RooHistError.cxx:43
 RooHistError.cxx:44
 RooHistError.cxx:45
 RooHistError.cxx:46
 RooHistError.cxx:47
 RooHistError.cxx:48
 RooHistError.cxx:49
 RooHistError.cxx:50
 RooHistError.cxx:51
 RooHistError.cxx:52
 RooHistError.cxx:53
 RooHistError.cxx:54
 RooHistError.cxx:55
 RooHistError.cxx:56
 RooHistError.cxx:57
 RooHistError.cxx:58
 RooHistError.cxx:59
 RooHistError.cxx:60
 RooHistError.cxx:61
 RooHistError.cxx:62
 RooHistError.cxx:63
 RooHistError.cxx:64
 RooHistError.cxx:65
 RooHistError.cxx:66
 RooHistError.cxx:67
 RooHistError.cxx:68
 RooHistError.cxx:69
 RooHistError.cxx:70
 RooHistError.cxx:71
 RooHistError.cxx:72
 RooHistError.cxx:73
 RooHistError.cxx:74
 RooHistError.cxx:75
 RooHistError.cxx:76
 RooHistError.cxx:77
 RooHistError.cxx:78
 RooHistError.cxx:79
 RooHistError.cxx:80
 RooHistError.cxx:81
 RooHistError.cxx:82
 RooHistError.cxx:83
 RooHistError.cxx:84
 RooHistError.cxx:85
 RooHistError.cxx:86
 RooHistError.cxx:87
 RooHistError.cxx:88
 RooHistError.cxx:89
 RooHistError.cxx:90
 RooHistError.cxx:91
 RooHistError.cxx:92
 RooHistError.cxx:93
 RooHistError.cxx:94
 RooHistError.cxx:95
 RooHistError.cxx:96
 RooHistError.cxx:97
 RooHistError.cxx:98
 RooHistError.cxx:99
 RooHistError.cxx:100
 RooHistError.cxx:101
 RooHistError.cxx:102
 RooHistError.cxx:103
 RooHistError.cxx:104
 RooHistError.cxx:105
 RooHistError.cxx:106
 RooHistError.cxx:107
 RooHistError.cxx:108
 RooHistError.cxx:109
 RooHistError.cxx:110
 RooHistError.cxx:111
 RooHistError.cxx:112
 RooHistError.cxx:113
 RooHistError.cxx:114
 RooHistError.cxx:115
 RooHistError.cxx:116
 RooHistError.cxx:117
 RooHistError.cxx:118
 RooHistError.cxx:119
 RooHistError.cxx:120
 RooHistError.cxx:121
 RooHistError.cxx:122
 RooHistError.cxx:123
 RooHistError.cxx:124
 RooHistError.cxx:125
 RooHistError.cxx:126
 RooHistError.cxx:127
 RooHistError.cxx:128
 RooHistError.cxx:129
 RooHistError.cxx:130
 RooHistError.cxx:131
 RooHistError.cxx:132
 RooHistError.cxx:133
 RooHistError.cxx:134
 RooHistError.cxx:135
 RooHistError.cxx:136
 RooHistError.cxx:137
 RooHistError.cxx:138
 RooHistError.cxx:139
 RooHistError.cxx:140
 RooHistError.cxx:141
 RooHistError.cxx:142
 RooHistError.cxx:143
 RooHistError.cxx:144
 RooHistError.cxx:145
 RooHistError.cxx:146
 RooHistError.cxx:147
 RooHistError.cxx:148
 RooHistError.cxx:149
 RooHistError.cxx:150
 RooHistError.cxx:151
 RooHistError.cxx:152
 RooHistError.cxx:153
 RooHistError.cxx:154
 RooHistError.cxx:155
 RooHistError.cxx:156
 RooHistError.cxx:157
 RooHistError.cxx:158
 RooHistError.cxx:159
 RooHistError.cxx:160
 RooHistError.cxx:161
 RooHistError.cxx:162
 RooHistError.cxx:163
 RooHistError.cxx:164
 RooHistError.cxx:165
 RooHistError.cxx:166
 RooHistError.cxx:167
 RooHistError.cxx:168
 RooHistError.cxx:169
 RooHistError.cxx:170
 RooHistError.cxx:171
 RooHistError.cxx:172
 RooHistError.cxx:173
 RooHistError.cxx:174
 RooHistError.cxx:175
 RooHistError.cxx:176
 RooHistError.cxx:177
 RooHistError.cxx:178
 RooHistError.cxx:179
 RooHistError.cxx:180
 RooHistError.cxx:181
 RooHistError.cxx:182
 RooHistError.cxx:183
 RooHistError.cxx:184
 RooHistError.cxx:185
 RooHistError.cxx:186
 RooHistError.cxx:187
 RooHistError.cxx:188
 RooHistError.cxx:189
 RooHistError.cxx:190
 RooHistError.cxx:191
 RooHistError.cxx:192
 RooHistError.cxx:193
 RooHistError.cxx:194
 RooHistError.cxx:195
 RooHistError.cxx:196
 RooHistError.cxx:197
 RooHistError.cxx:198
 RooHistError.cxx:199
 RooHistError.cxx:200
 RooHistError.cxx:201
 RooHistError.cxx:202
 RooHistError.cxx:203
 RooHistError.cxx:204
 RooHistError.cxx:205
 RooHistError.cxx:206
 RooHistError.cxx:207
 RooHistError.cxx:208
 RooHistError.cxx:209
 RooHistError.cxx:210
 RooHistError.cxx:211
 RooHistError.cxx:212
 RooHistError.cxx:213
 RooHistError.cxx:214
 RooHistError.cxx:215
 RooHistError.cxx:216
 RooHistError.cxx:217
 RooHistError.cxx:218
 RooHistError.cxx:219
 RooHistError.cxx:220
 RooHistError.cxx:221
 RooHistError.cxx:222
 RooHistError.cxx:223
 RooHistError.cxx:224
 RooHistError.cxx:225
 RooHistError.cxx:226
 RooHistError.cxx:227
 RooHistError.cxx:228
 RooHistError.cxx:229
 RooHistError.cxx:230
 RooHistError.cxx:231
 RooHistError.cxx:232
 RooHistError.cxx:233
 RooHistError.cxx:234
 RooHistError.cxx:235
 RooHistError.cxx:236
 RooHistError.cxx:237
 RooHistError.cxx:238
 RooHistError.cxx:239
 RooHistError.cxx:240
 RooHistError.cxx:241
 RooHistError.cxx:242
 RooHistError.cxx:243
 RooHistError.cxx:244
 RooHistError.cxx:245
 RooHistError.cxx:246
 RooHistError.cxx:247
 RooHistError.cxx:248
 RooHistError.cxx:249
 RooHistError.cxx:250
 RooHistError.cxx:251
 RooHistError.cxx:252
 RooHistError.cxx:253
 RooHistError.cxx:254
 RooHistError.cxx:255
 RooHistError.cxx:256
 RooHistError.cxx:257
 RooHistError.cxx:258
 RooHistError.cxx:259
 RooHistError.cxx:260
 RooHistError.cxx:261
 RooHistError.cxx:262
 RooHistError.cxx:263
 RooHistError.cxx:264
 RooHistError.cxx:265
 RooHistError.cxx:266
 RooHistError.cxx:267
 RooHistError.cxx:268
 RooHistError.cxx:269
 RooHistError.cxx:270
 RooHistError.cxx:271
 RooHistError.cxx:272
 RooHistError.cxx:273
 RooHistError.cxx:274
 RooHistError.cxx:275
 RooHistError.cxx:276
 RooHistError.cxx:277
 RooHistError.cxx:278
 RooHistError.cxx:279
 RooHistError.cxx:280
 RooHistError.cxx:281
 RooHistError.cxx:282
 RooHistError.cxx:283
 RooHistError.cxx:284
 RooHistError.cxx:285
 RooHistError.cxx:286
 RooHistError.cxx:287
 RooHistError.cxx:288
 RooHistError.cxx:289
 RooHistError.cxx:290
 RooHistError.cxx:291
 RooHistError.cxx:292
 RooHistError.cxx:293
 RooHistError.cxx:294
 RooHistError.cxx:295
 RooHistError.cxx:296
 RooHistError.cxx:297
 RooHistError.cxx:298
 RooHistError.cxx:299
 RooHistError.cxx:300
 RooHistError.cxx:301
 RooHistError.cxx:302
 RooHistError.cxx:303
 RooHistError.cxx:304
 RooHistError.cxx:305
 RooHistError.cxx:306
 RooHistError.cxx:307
 RooHistError.cxx:308
 RooHistError.cxx:309
 RooHistError.cxx:310
 RooHistError.cxx:311
 RooHistError.cxx:312
 RooHistError.cxx:313
 RooHistError.cxx:314
 RooHistError.cxx:315
 RooHistError.cxx:316
 RooHistError.cxx:317
 RooHistError.cxx:318
 RooHistError.cxx:319
 RooHistError.cxx:320
 RooHistError.cxx:321
 RooHistError.cxx:322
 RooHistError.cxx:323
 RooHistError.cxx:324
 RooHistError.cxx:325
 RooHistError.cxx:326
 RooHistError.cxx:327
 RooHistError.cxx:328
 RooHistError.cxx:329
 RooHistError.cxx:330
 RooHistError.cxx:331
 RooHistError.cxx:332
 RooHistError.cxx:333
 RooHistError.cxx:334
 RooHistError.cxx:335
 RooHistError.cxx:336
 RooHistError.cxx:337
 RooHistError.cxx:338
 RooHistError.cxx:339
 RooHistError.cxx:340
 RooHistError.cxx:341
 RooHistError.cxx:342
 RooHistError.cxx:343
 RooHistError.cxx:344
 RooHistError.cxx:345
 RooHistError.cxx:346