Logo ROOT   6.10/09
Reference Guide
RooHistError.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 /**
18 \file RooHistError.cxx
19 \class RooHistError
20 \ingroup Roofitcore
21 
22 RooHistError is a singleton class used to calculate the error bars
23 for each bin of a RooHist object. Errors are calculated by integrating
24 a specified area of a Poisson or Binomail error distribution.
25 **/
26 
27 #include "RooFit.h"
28 
29 #include "RooHistError.h"
30 #include "RooBrentRootFinder.h"
31 #include "RooMsgService.h"
32 #include "TMath.h"
33 
34 #include "Riostream.h"
35 #include <assert.h>
36 
37 using namespace std;
38 
40  ;
41 
42 
43 
44 ////////////////////////////////////////////////////////////////////////////////
45 /// Return a reference to a singleton object that is created the
46 /// first time this method is called. Only one object will be
47 /// constructed per ROOT session.
48 
50 {
51  static RooHistError _theInstance;
52  return _theInstance;
53 }
54 
55 
56 ////////////////////////////////////////////////////////////////////////////////
57 /// Construct our singleton object.
58 
60 {
61  // Initialize lookup table ;
62  Int_t i ;
63  for (i=0 ; i<1000 ; i++) {
64  getPoissonIntervalCalc(i,_poissonLoLUT[i],_poissonHiLUT[i],1.) ;
65  }
66 
67 }
68 
69 
70 
71 ////////////////////////////////////////////////////////////////////////////////
72 /// Return a confidence interval for the expected number of events given n
73 /// observed (unweighted) events. The interval will contain the same probability
74 /// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
75 /// the point estimate n (ie, for small n) in which case the interval is adjusted
76 /// to start at n. This method uses a lookup table to return precalculated results
77 /// for n<1000
78 
80 {
81  // Use lookup table for most common cases
82  if (n<1000 && nSigma==1.) {
83  mu1=_poissonLoLUT[n] ;
84  mu2=_poissonHiLUT[n] ;
85  return kTRUE ;
86  }
87 
88  // Forward to calculation method
89  Bool_t ret = getPoissonIntervalCalc(n,mu1,mu2,nSigma) ;
90  return ret ;
91 }
92 
93 
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Calculate a confidence interval for the expected number of events given n
97 /// observed (unweighted) events. The interval will contain the same probability
98 /// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
99 /// the point estimate n (ie, for small n) in which case the interval is adjusted
100 /// to start at n.
101 
103 {
104  // sanity checks
105  if(n < 0) {
106  oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n << endl;
107  return kFALSE;
108  }
109 
110  // use assymptotic error if possible
111  if(n > 100) {
112  mu1= n - sqrt(n+0.25) + 0.5;
113  mu2= n + sqrt(n+0.25) + 0.5;
114  return kTRUE;
115  }
116 
117  // create a function object to use
118  PoissonSum upper(n);
119  if(n > 0) {
120  PoissonSum lower(n-1);
121  return getInterval(&upper,&lower,(Double_t)n,1.0,mu1,mu2,nSigma);
122  }
123 
124  // Backup solution for negative numbers
125  return getInterval(&upper,0,(Double_t)n,1.0,mu1,mu2,nSigma);
126 }
127 
128 
129 ////////////////////////////////////////////////////////////////////////////////
130 /// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
131 /// If the return values is kFALSE and error occurred.
132 
134  Double_t &asym1, Double_t &asym2, Double_t nSigma) const
135 {
136  // sanity checks
137  if(n < 0 || m < 0) {
138  oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
139  return kFALSE;
140  }
141 
142  // handle the special case of no events in either category
143  if(n == 0 && m == 0) {
144  asym1= -1;
145  asym2= +1;
146  return kTRUE;
147  }
148 
149  // handle cases when n,m>100 (factorials in BinomialSum will overflow around 170)
150  if ((n>100&&m>100)) {
151  Double_t N = n ;
152  Double_t M = m ;
153  Double_t asym = 1.0*(N-M)/(N+M) ;
154  Double_t approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
155 
156  asym1 = asym-nSigma*approxErr ;
157  asym2 = asym+nSigma*approxErr ;
158  return kTRUE ;
159  }
160 
161  // swap n and m to ensure that n <= m
162  Bool_t swapped(kFALSE);
163  if(n > m) {
164  swapped= kTRUE;
165  Int_t tmp(m);
166  m= n;
167  n= tmp;
168  }
169 
170  // create the function objects to use
171  Bool_t status(kFALSE);
172  BinomialSumAsym upper(n,m);
173  if(n > 0) {
174  BinomialSumAsym lower(n-1,m+1);
175  status= getInterval(&upper,&lower,(Double_t)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
176  }
177  else {
178  status= getInterval(&upper,0,(Double_t)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
179  }
180 
181  // undo the swap here
182  if(swapped) {
183  Double_t tmp(asym1);
184  asym1= -asym2;
185  asym2= -tmp;
186  }
187 
188  return status;
189 }
190 
191 
192 ////////////////////////////////////////////////////////////////////////////////
193 /// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
194 /// If the return values is kFALSE and error occurred.
195 
197  Double_t &asym1, Double_t &asym2, Double_t nSigma) const
198 {
199  // sanity checks
200  if(n < 0 || m < 0) {
201  oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
202  return kFALSE;
203  }
204 
205  // handle the special case of no events in either category
206  if(n == 0 && m == 0) {
207  asym1= -1;
208  asym2= +1;
209  return kTRUE;
210  }
211 
212  // handle cases when n,m>80 (factorials in BinomialSum will overflow around 170)
213  if ((n>80&&m>80)) {
214  Double_t N = n ;
215  Double_t M = m ;
216  Double_t asym = 1.0*(N)/(N+M) ;
217  Double_t approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
218 
219  asym1 = asym-nSigma*0.5*approxErr ;
220  asym2 = asym+nSigma*0.5*approxErr ;
221  return kTRUE ;
222  }
223 
224  // swap n and m to ensure that n <= m
225  Bool_t swapped(kFALSE);
226  if(n > m) {
227  swapped= kTRUE;
228  Int_t tmp(m);
229  m= n;
230  n= tmp;
231  }
232 
233  // create the function objects to use
234  Bool_t status(kFALSE);
235  BinomialSumEff upper(n,m);
236  Double_t eff = (Double_t)(n)/(n+m) ;
237  if(n > 0) {
238  BinomialSumEff lower(n-1,m+1);
239  status= getInterval(&upper,&lower,eff,0.1,asym1,asym2,nSigma*0.5);
240  }
241  else {
242  status= getInterval(&upper,0,eff,0.1,asym1,asym2,nSigma*0.5);
243  }
244 
245  // undo the swap here
246  if(swapped) {
247  Double_t tmp(asym1);
248  asym1= 1-asym2;
249  asym2= 1-tmp;
250  }
251 
252  return status;
253 }
254 
255 
256 
257 ////////////////////////////////////////////////////////////////////////////////
258 /// Calculate a confidence interval using the cumulative functions provided.
259 /// The interval will be "central" when both cumulative functions are provided,
260 /// unless this would exclude the pointEstimate, in which case a one-sided interval
261 /// pinned at the point estimate is returned instead.
262 
263 Bool_t RooHistError::getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, Double_t pointEstimate,
264  Double_t stepSize, Double_t &lo, Double_t &hi, Double_t nSigma) const
265 {
266  // sanity checks
267  assert(0 != Qu || 0 != Ql);
268 
269  // convert number of sigma into a confidence level
270  Double_t beta= TMath::Erf(nSigma/sqrt(2.));
271  Double_t alpha= 0.5*(1-beta);
272 
273  // Does the central interval contain the point estimate?
274  Bool_t ok(kTRUE);
275  Double_t loProb(1),hiProb(0);
276  if(0 != Ql) loProb= (*Ql)(&pointEstimate);
277  if(0 != Qu) hiProb= (*Qu)(&pointEstimate);
278 
279  if (Qu && (0 == Ql || loProb > alpha + beta)) {
280  // Force the low edge to be at the pointEstimate
281  lo= pointEstimate;
282  Double_t target= loProb - beta;
283  hi= seek(*Qu,lo,+stepSize,target);
284  RooBrentRootFinder uFinder(*Qu);
285  ok= uFinder.findRoot(hi,hi-stepSize,hi,target);
286  }
287  else if(Ql && (0 == Qu || hiProb < alpha)) {
288  // Force the high edge to be at pointEstimate
289  hi= pointEstimate;
290  Double_t target= hiProb + beta;
291  lo= seek(*Ql,hi,-stepSize,target);
292  RooBrentRootFinder lFinder(*Ql);
293  ok= lFinder.findRoot(lo,lo,lo+stepSize,target);
294  }
295  else if (Qu && Ql) {
296  // use a central interval
297  lo= seek(*Ql,pointEstimate,-stepSize,alpha+beta);
298  hi= seek(*Qu,pointEstimate,+stepSize,alpha);
299  RooBrentRootFinder lFinder(*Ql),uFinder(*Qu);
300  ok= lFinder.findRoot(lo,lo,lo+stepSize,alpha+beta);
301  ok|= uFinder.findRoot(hi,hi-stepSize,hi,alpha);
302  }
303  if(!ok) oocoutE((TObject*)0,Plotting) << "RooHistError::getInterval: failed to find root(s)" << endl;
304 
305  return ok;
306 }
307 
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 /// Scan f(x)-value until it changes sign. Start at the specified point and take constant
311 /// steps of the specified size. Give up after 1000 steps.
312 
313 Double_t RooHistError::seek(const RooAbsFunc &f, Double_t startAt, Double_t step, Double_t value) const
314 {
315  Int_t steps(1000);
316  Double_t min(f.getMinLimit(1)),max(f.getMaxLimit(1));
317  Double_t x(startAt), f0= f(&startAt) - value;
318  do {
319  x+= step;
320  }
321  while(steps-- && (f0*(f(&x)-value) >= 0) && ((x-min)*(max-x) >= 0));
322  assert(0 != steps);
323  if(x < min) x= min;
324  if(x > max) x= max;
325 
326  return x;
327 }
328 
329 
330 
331 ////////////////////////////////////////////////////////////////////////////////
332 /// Create and return a PoissonSum function binding
333 
335 {
336  return new PoissonSum(n);
337 }
338 
339 
340 ////////////////////////////////////////////////////////////////////////////////
341 /// Create and return a BinomialSum function binding
342 
344 {
345  if (eff) {
346  return new BinomialSumEff(n,m) ;
347  } else {
348  return new BinomialSumAsym(n,m) ;
349  }
350 }
#define N
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
static RooAbsFunc * createPoissonSum(Int_t n)
Create and return a PoissonSum function binding.
STL namespace.
double beta(double x, double y)
Calculates the beta function.
double sqrt(double)
Double_t seek(const RooAbsFunc &f, Double_t startAt, Double_t step, Double_t value) const
Scan f(x)-value until it changes sign.
Double_t x[n]
Definition: legend1.C:17
#define oocoutE(o, a)
Definition: RooMsgService.h:47
Bool_t 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.
virtual Bool_t findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value=0) const
Do the root finding using the Brent-Decker method.
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:187
TMarker * m
Definition: textangle.C:8
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method...
virtual Double_t getMinLimit(UInt_t dimension) const =0
static const RooHistError & instance()
Return a reference to a singleton object that is created the first time this method is called...
RooHistError()
Construct our singleton object.
const Bool_t kFALSE
Definition: RtypesCore.h:92
#define ClassImp(name)
Definition: Rtypes.h:336
double f(double x)
virtual Double_t getMaxLimit(UInt_t dimension) const =0
Bool_t getBinomialIntervalAsym(Int_t n, Int_t m, Double_t &a1, Double_t &a2, Double_t nSigma=1) const
Return &#39;nSigma&#39; binomial confidence interval for (n,m).
double Double_t
Definition: RtypesCore.h:55
Bool_t getPoissonInterval(Int_t n, Double_t &mu1, Double_t &mu2, Double_t nSigma=1) const
Return a confidence interval for the expected number of events given n observed (unweighted) events...
Mother of all ROOT objects.
Definition: TObject.h:37
static RooAbsFunc * createBinomialSum(Int_t n, Int_t m, Bool_t eff)
Create and return a BinomialSum function binding.
Bool_t getBinomialIntervalEff(Int_t n, Int_t m, Double_t &a1, Double_t &a2, Double_t nSigma=1) const
Return &#39;nSigma&#39; binomial confidence interval for (n,m).
float type_of_call hi(const int &, const int &)
Bool_t getPoissonIntervalCalc(Int_t n, Double_t &mu1, Double_t &mu2, Double_t nSigma=1) const
Calculate a confidence interval for the expected number of events given n observed (unweighted) event...
const Bool_t kTRUE
Definition: RtypesCore.h:91
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
const Int_t n
Definition: legend1.C:16
RooHistError is a singleton class used to calculate the error bars for each bin of a RooHist object...
Definition: RooHistError.h:25