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
77bool RooHistError::getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma) const
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 true ;
84 }
85
86 // Forward to calculation method
87 bool 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
100bool RooHistError::getPoissonIntervalCalc(Int_t n, double &mu1, double &mu2, double nSigma) const
101{
102 // sanity checks
103 if(n < 0) {
104 oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n << endl;
105 return false;
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 true;
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)n,1.0,mu1,mu2,nSigma);
120 }
121
122 // Backup solution for negative numbers
123 return getInterval(&upper,0,(double)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 false and error occurred.
130
132 double &asym1, double &asym2, double nSigma) const
133{
134 // sanity checks
135 if(n < 0 || m < 0) {
136 oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
137 return false;
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 true;
145 }
146
147 // handle cases when n,m>100 (factorials in BinomialSum will overflow around 170)
148 if ((n>100&&m>100)) {
149 double N = n ;
150 double M = m ;
151 double asym = 1.0*(N-M)/(N+M) ;
152 double 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 true ;
157 }
158
159 // swap n and m to ensure that n <= m
160 bool swapped(false);
161 if(n > m) {
162 swapped= true;
163 Int_t tmp(m);
164 m= n;
165 n= tmp;
166 }
167
168 // create the function objects to use
169 bool status(false);
170 BinomialSumAsym upper(n,m);
171 if(n > 0) {
172 BinomialSumAsym lower(n-1,m+1);
173 status= getInterval(&upper,&lower,(double)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
174 }
175 else {
176 status= getInterval(&upper,0,(double)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
177 }
178
179 // undo the swap here
180 if(swapped) {
181 double 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 false and error occurred.
193
195 double &asym1, double &asym2, double nSigma) const
196{
197 // sanity checks
198 if(n < 0 || m < 0) {
199 oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
200 return false;
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 true;
208 }
209
210 // handle cases when n,m>80 (factorials in BinomialSum will overflow around 170)
211 if ((n>80&&m>80)) {
212 double N = n ;
213 double M = m ;
214 double asym = 1.0*(N)/(N+M) ;
215 double 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 true ;
220 }
221
222 // swap n and m to ensure that n <= m
223 bool swapped(false);
224 if(n > m) {
225 swapped= true;
226 Int_t tmp(m);
227 m= n;
228 n= tmp;
229 }
230
231 // create the function objects to use
232 bool status(false);
233 BinomialSumEff upper(n,m);
234 double eff = (double)(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 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
261bool RooHistError::getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, double pointEstimate,
262 double stepSize, double &lo, double &hi, double nSigma) const
263{
264 // sanity checks
265 assert(0 != Qu || 0 != Ql);
266
267 // convert number of sigma into a confidence level
268 double beta= TMath::Erf(nSigma/sqrt(2.));
269 double alpha= 0.5*(1-beta);
270
271 // Does the central interval contain the point estimate?
272 bool ok(true);
273 double 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 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 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(nullptr,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
311double RooHistError::seek(const RooAbsFunc &f, double startAt, double step, double value) const
312{
313 Int_t steps(1000);
314 double min(f.getMinLimit(1)),max(f.getMaxLimit(1));
315 double 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:52
#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 findRoot(double &result, double xlo, double xhi, double 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 getBinomialIntervalAsym(Int_t n, Int_t m, double &a1, double &a2, double nSigma=1) const
Return 'nSigma' binomial confidence interval for (n,m).
double _poissonLoLUT[1000]
Definition: RooHistError.h:43
static const RooHistError & instance()
Return a reference to a singleton object that is created the first time this method is called.
bool getBinomialIntervalEff(Int_t n, Int_t m, double &a1, double &a2, double nSigma=1) const
Return 'nSigma' binomial confidence interval for (n,m).
bool getPoissonIntervalCalc(Int_t n, double &mu1, double &mu2, double nSigma=1) const
Calculate a confidence interval for the expected number of events given n observed (unweighted) event...
bool getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, double pointEstimate, double stepSize, double &lo, double &hi, double nSigma) const
Calculate a confidence interval using the cumulative functions provided.
double seek(const RooAbsFunc &f, double startAt, double step, double value) const
Scan f(x)-value until it changes sign.
bool getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma=1) const
Return a confidence interval for the expected number of events given n observed (unweighted) events.
static RooAbsFunc * createPoissonSum(Int_t n)
Create and return a PoissonSum function binding.
RooHistError()
Construct our singleton object.
double _poissonHiLUT[1000]
Definition: RooHistError.h:44
static RooAbsFunc * createBinomialSum(Int_t n, Int_t m, bool eff)
Create and return a BinomialSum function binding.
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:190
TMarker m
Definition: textangle.C:8