Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
22RooHistError is a singleton class used to calculate the error bars
23for each bin of a RooHist object. Errors are calculated by integrating
24a 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
37using 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++) {
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
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
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 f(i)
Definition RSha256.hxx:104
#define oocoutE(o, a)
const Bool_t kFALSE
Definition RtypesCore.h:101
const Bool_t kTRUE
Definition RtypesCore.h:100
#define ClassImp(name)
Definition Rtypes.h:364
#define N
#define hi
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method.
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.
RooHistError is a singleton class used to calculate the error bars for each bin of a RooHist object.
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.
Bool_t getBinomialIntervalEff(Int_t n, Int_t m, Double_t &a1, Double_t &a2, Double_t nSigma=1) const
Return 'nSigma' binomial confidence interval for (n,m).
static const RooHistError & instance()
Return a reference to a singleton object that is created the first time this method is called.
static RooAbsFunc * createBinomialSum(Int_t n, Int_t m, Bool_t eff)
Create and return a BinomialSum function binding.
Bool_t getBinomialIntervalAsym(Int_t n, Int_t m, Double_t &a1, Double_t &a2, Double_t nSigma=1) const
Return 'nSigma' binomial confidence interval for (n,m).
Double_t _poissonLoLUT[1000]
Double_t seek(const RooAbsFunc &f, Double_t startAt, Double_t step, Double_t value) const
Scan f(x)-value until it changes sign.
static RooAbsFunc * createPoissonSum(Int_t n)
Create and return a PoissonSum function binding.
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.
RooHistError()
Construct our singleton object.
Double_t _poissonHiLUT[1000]
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...
Mother of all ROOT objects.
Definition TObject.h:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition TMath.cxx:184
auto * m
Definition textangle.C:8