Logo ROOT  
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
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 "RooHistError.h"
28#include "RooBrentRootFinder.h"
29#include "RooMsgService.h"
30#include "TMath.h"
31
32#include "Riostream.h"
33#include <assert.h>
34
35using namespace std;
36
38 ;
39
40
41
42////////////////////////////////////////////////////////////////////////////////
43/// Return a reference to a singleton object that is created the
44/// first time this method is called. Only one object will be
45/// constructed per ROOT session.
46
48{
49 static RooHistError _theInstance;
50 return _theInstance;
51}
52
53
54////////////////////////////////////////////////////////////////////////////////
55/// Construct our singleton object.
56
58{
59 // Initialize lookup table ;
60 Int_t i ;
61 for (i=0 ; i<1000 ; i++) {
63 }
64
65}
66
67
68
69////////////////////////////////////////////////////////////////////////////////
70/// Return a confidence interval for the expected number of events given n
71/// observed (unweighted) events. The interval will contain the same probability
72/// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
73/// the point estimate n (ie, for small n) in which case the interval is adjusted
74/// to start at n. This method uses a lookup table to return precalculated results
75/// for n<1000
76
78{
79 // Use lookup table for most common cases
80 if (n<1000 && nSigma==1.) {
81 mu1=_poissonLoLUT[n] ;
82 mu2=_poissonHiLUT[n] ;
83 return kTRUE ;
84 }
85
86 // Forward to calculation method
87 Bool_t ret = getPoissonIntervalCalc(n,mu1,mu2,nSigma) ;
88 return ret ;
89}
90
91
92
93////////////////////////////////////////////////////////////////////////////////
94/// Calculate a confidence interval for the expected number of events given n
95/// observed (unweighted) events. The interval will contain the same probability
96/// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
97/// the point estimate n (ie, for small n) in which case the interval is adjusted
98/// to start at n.
99
101{
102 // sanity checks
103 if(n < 0) {
104 oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n << endl;
105 return kFALSE;
106 }
107
108 // use assymptotic error if possible
109 if(n > 100) {
110 mu1= n - sqrt(n+0.25) + 0.5;
111 mu2= n + sqrt(n+0.25) + 0.5;
112 return kTRUE;
113 }
114
115 // create a function object to use
116 PoissonSum upper(n);
117 if(n > 0) {
118 PoissonSum lower(n-1);
119 return getInterval(&upper,&lower,(Double_t)n,1.0,mu1,mu2,nSigma);
120 }
121
122 // Backup solution for negative numbers
123 return getInterval(&upper,0,(Double_t)n,1.0,mu1,mu2,nSigma);
124}
125
126
127////////////////////////////////////////////////////////////////////////////////
128/// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
129/// If the return values is kFALSE and error occurred.
130
132 Double_t &asym1, Double_t &asym2, Double_t nSigma) const
133{
134 // sanity checks
135 if(n < 0 || m < 0) {
136 oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
137 return kFALSE;
138 }
139
140 // handle the special case of no events in either category
141 if(n == 0 && m == 0) {
142 asym1= -1;
143 asym2= +1;
144 return kTRUE;
145 }
146
147 // handle cases when n,m>100 (factorials in BinomialSum will overflow around 170)
148 if ((n>100&&m>100)) {
149 Double_t N = n ;
150 Double_t M = m ;
151 Double_t asym = 1.0*(N-M)/(N+M) ;
152 Double_t approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
153
154 asym1 = asym-nSigma*approxErr ;
155 asym2 = asym+nSigma*approxErr ;
156 return kTRUE ;
157 }
158
159 // swap n and m to ensure that n <= m
160 Bool_t swapped(kFALSE);
161 if(n > m) {
162 swapped= kTRUE;
163 Int_t tmp(m);
164 m= n;
165 n= tmp;
166 }
167
168 // create the function objects to use
169 Bool_t status(kFALSE);
170 BinomialSumAsym upper(n,m);
171 if(n > 0) {
172 BinomialSumAsym lower(n-1,m+1);
173 status= getInterval(&upper,&lower,(Double_t)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
174 }
175 else {
176 status= getInterval(&upper,0,(Double_t)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
177 }
178
179 // undo the swap here
180 if(swapped) {
181 Double_t tmp(asym1);
182 asym1= -asym2;
183 asym2= -tmp;
184 }
185
186 return status;
187}
188
189
190////////////////////////////////////////////////////////////////////////////////
191/// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
192/// If the return values is kFALSE and error occurred.
193
195 Double_t &asym1, Double_t &asym2, Double_t nSigma) const
196{
197 // sanity checks
198 if(n < 0 || m < 0) {
199 oocoutE((TObject*)0,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
200 return kFALSE;
201 }
202
203 // handle the special case of no events in either category
204 if(n == 0 && m == 0) {
205 asym1= -1;
206 asym2= +1;
207 return kTRUE;
208 }
209
210 // handle cases when n,m>80 (factorials in BinomialSum will overflow around 170)
211 if ((n>80&&m>80)) {
212 Double_t N = n ;
213 Double_t M = m ;
214 Double_t asym = 1.0*(N)/(N+M) ;
215 Double_t approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
216
217 asym1 = asym-nSigma*0.5*approxErr ;
218 asym2 = asym+nSigma*0.5*approxErr ;
219 return kTRUE ;
220 }
221
222 // swap n and m to ensure that n <= m
223 Bool_t swapped(kFALSE);
224 if(n > m) {
225 swapped= kTRUE;
226 Int_t tmp(m);
227 m= n;
228 n= tmp;
229 }
230
231 // create the function objects to use
232 Bool_t status(kFALSE);
233 BinomialSumEff upper(n,m);
234 Double_t eff = (Double_t)(n)/(n+m) ;
235 if(n > 0) {
236 BinomialSumEff lower(n-1,m+1);
237 status= getInterval(&upper,&lower,eff,0.1,asym1,asym2,nSigma*0.5);
238 }
239 else {
240 status= getInterval(&upper,0,eff,0.1,asym1,asym2,nSigma*0.5);
241 }
242
243 // undo the swap here
244 if(swapped) {
245 Double_t tmp(asym1);
246 asym1= 1-asym2;
247 asym2= 1-tmp;
248 }
249
250 return status;
251}
252
253
254
255////////////////////////////////////////////////////////////////////////////////
256/// Calculate a confidence interval using the cumulative functions provided.
257/// The interval will be "central" when both cumulative functions are provided,
258/// unless this would exclude the pointEstimate, in which case a one-sided interval
259/// pinned at the point estimate is returned instead.
260
262 Double_t stepSize, Double_t &lo, Double_t &hi, Double_t nSigma) const
263{
264 // sanity checks
265 assert(0 != Qu || 0 != Ql);
266
267 // convert number of sigma into a confidence level
268 Double_t beta= TMath::Erf(nSigma/sqrt(2.));
269 Double_t alpha= 0.5*(1-beta);
270
271 // Does the central interval contain the point estimate?
272 Bool_t ok(kTRUE);
273 Double_t loProb(1),hiProb(0);
274 if(0 != Ql) loProb= (*Ql)(&pointEstimate);
275 if(0 != Qu) hiProb= (*Qu)(&pointEstimate);
276
277 if (Qu && (0 == Ql || loProb > alpha + beta)) {
278 // Force the low edge to be at the pointEstimate
279 lo= pointEstimate;
280 Double_t target= loProb - beta;
281 hi= seek(*Qu,lo,+stepSize,target);
282 RooBrentRootFinder uFinder(*Qu);
283 ok= uFinder.findRoot(hi,hi-stepSize,hi,target);
284 }
285 else if(Ql && (0 == Qu || hiProb < alpha)) {
286 // Force the high edge to be at pointEstimate
287 hi= pointEstimate;
288 Double_t target= hiProb + beta;
289 lo= seek(*Ql,hi,-stepSize,target);
290 RooBrentRootFinder lFinder(*Ql);
291 ok= lFinder.findRoot(lo,lo,lo+stepSize,target);
292 }
293 else if (Qu && Ql) {
294 // use a central interval
295 lo= seek(*Ql,pointEstimate,-stepSize,alpha+beta);
296 hi= seek(*Qu,pointEstimate,+stepSize,alpha);
297 RooBrentRootFinder lFinder(*Ql),uFinder(*Qu);
298 ok= lFinder.findRoot(lo,lo,lo+stepSize,alpha+beta);
299 ok|= uFinder.findRoot(hi,hi-stepSize,hi,alpha);
300 }
301 if(!ok) oocoutE((TObject*)0,Plotting) << "RooHistError::getInterval: failed to find root(s)" << endl;
302
303 return ok;
304}
305
306
307////////////////////////////////////////////////////////////////////////////////
308/// Scan f(x)-value until it changes sign. Start at the specified point and take constant
309/// steps of the specified size. Give up after 1000 steps.
310
312{
313 Int_t steps(1000);
314 Double_t min(f.getMinLimit(1)),max(f.getMaxLimit(1));
315 Double_t x(startAt), f0= f(&startAt) - value;
316 do {
317 x+= step;
318 }
319 while(steps-- && (f0*(f(&x)-value) >= 0) && ((x-min)*(max-x) >= 0));
320 assert(0 != steps);
321 if(x < min) x= min;
322 if(x > max) x= max;
323
324 return x;
325}
326
327
328
329////////////////////////////////////////////////////////////////////////////////
330/// Create and return a PoissonSum function binding
331
333{
334 return new PoissonSum(n);
335}
336
337
338////////////////////////////////////////////////////////////////////////////////
339/// Create and return a BinomialSum function binding
340
342{
343 if (eff) {
344 return new BinomialSumEff(n,m) ;
345 } else {
346 return new BinomialSumAsym(n,m) ;
347 }
348}
#define f(i)
Definition: RSha256.hxx:104
#define oocoutE(o, a)
Definition: RooMsgService.h:48
const Bool_t kFALSE
Definition: RtypesCore.h:101
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:375
#define N
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t target
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
#define hi
Definition: THbookFile.cxx:128
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.
Bool_t findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value=0) const override
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.
Definition: RooHistError.h:25
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]
Definition: RooHistError.h:43
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]
Definition: RooHistError.h:44
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:37
double beta(double x, double y)
Calculates the beta function.
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:184
auto * m
Definition: textangle.C:8