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
22Singleton 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 <cassert>
34
35using namespace std;
36
37
38////////////////////////////////////////////////////////////////////////////////
39/// Return a reference to a singleton object that is created the
40/// first time this method is called. Only one object will be
41/// constructed per ROOT session.
42
44{
45 static RooHistError _theInstance;
46 return _theInstance;
47}
48
49
50////////////////////////////////////////////////////////////////////////////////
51/// Construct our singleton object.
52
54{
55 // Initialize lookup table ;
56 Int_t i ;
57 for (i=0 ; i<1000 ; i++) {
59 }
60
61}
62
63
64
65////////////////////////////////////////////////////////////////////////////////
66/// Return a confidence interval for the expected number of events given n
67/// observed (unweighted) events. The interval will contain the same probability
68/// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
69/// the point estimate n (ie, for small n) in which case the interval is adjusted
70/// to start at n. This method uses a lookup table to return precalculated results
71/// for n<1000
72
73bool RooHistError::getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma) const
74{
75 // Use lookup table for most common cases
76 if (n<1000 && nSigma==1.) {
77 mu1=_poissonLoLUT[n] ;
78 mu2=_poissonHiLUT[n] ;
79 return true ;
80 }
81
82 // Forward to calculation method
83 bool ret = getPoissonIntervalCalc(n,mu1,mu2,nSigma) ;
84 return ret ;
85}
86
87
88
89////////////////////////////////////////////////////////////////////////////////
90/// Calculate a confidence interval for the expected number of events given n
91/// observed (unweighted) events. The interval will contain the same probability
92/// as nSigma of a Gaussian. Uses a central interval unless this does not enclose
93/// the point estimate n (ie, for small n) in which case the interval is adjusted
94/// to start at n.
95
96bool RooHistError::getPoissonIntervalCalc(Int_t n, double &mu1, double &mu2, double nSigma) const
97{
98 // sanity checks
99 if(n < 0) {
100 oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n = " << n << endl;
101 return false;
102 }
103
104 // use asymptotic error if possible
105 if(n > 100) {
106 mu1= n - sqrt(n+0.25) + 0.5;
107 mu2= n + sqrt(n+0.25) + 0.5;
108 return true;
109 }
110
111 // create a function object to use
112 PoissonSum upper(n);
113 if(n > 0) {
114 PoissonSum lower(n-1);
115 return getInterval(&upper,&lower,(double)n,1.0,mu1,mu2,nSigma);
116 }
117
118 // Backup solution for negative numbers
119 return getInterval(&upper,nullptr,(double)n,1.0,mu1,mu2,nSigma);
120}
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
125/// If the return values is false and error occurred.
126
128 double &asym1, double &asym2, double nSigma) const
129{
130 // sanity checks
131 if(n < 0 || m < 0) {
132 oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
133 return false;
134 }
135
136 // handle the special case of no events in either category
137 if(n == 0 && m == 0) {
138 asym1= -1;
139 asym2= +1;
140 return true;
141 }
142
143 // handle cases when n,m>100 (factorials in BinomialSum will overflow around 170)
144 if ((n>100&&m>100)) {
145 double N = n ;
146 double M = m ;
147 double asym = 1.0*(N-M)/(N+M) ;
148 double approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
149
150 asym1 = asym-nSigma*approxErr ;
151 asym2 = asym+nSigma*approxErr ;
152 return true ;
153 }
154
155 // swap n and m to ensure that n <= m
156 bool swapped(false);
157 if(n > m) {
158 swapped= true;
159 Int_t tmp(m);
160 m= n;
161 n= tmp;
162 }
163
164 // create the function objects to use
165 bool status(false);
166 BinomialSumAsym upper(n,m);
167 if(n > 0) {
168 BinomialSumAsym lower(n-1,m+1);
169 status= getInterval(&upper,&lower,(double)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
170 }
171 else {
172 status= getInterval(&upper,nullptr,(double)(n-m)/(n+m),0.1,asym1,asym2,nSigma);
173 }
174
175 // undo the swap here
176 if(swapped) {
177 double tmp(asym1);
178 asym1= -asym2;
179 asym2= -tmp;
180 }
181
182 return status;
183}
184
185
186////////////////////////////////////////////////////////////////////////////////
187/// Return 'nSigma' binomial confidence interval for (n,m). The result is return in asym1 and asym2.
188/// If the return values is false and error occurred.
189
191 double &asym1, double &asym2, double nSigma) const
192{
193 // sanity checks
194 if(n < 0 || m < 0) {
195 oocoutE(nullptr,Plotting) << "RooHistError::getPoissonInterval: cannot calculate interval for n,m = " << n << "," << m << endl;
196 return false;
197 }
198
199 // handle the special case of no events in either category
200 if(n == 0 && m == 0) {
201 asym1= -1;
202 asym2= +1;
203 return true;
204 }
205
206 // handle cases when n,m>80 (factorials in BinomialSum will overflow around 170)
207 if ((n>80&&m>80)) {
208 double N = n ;
209 double M = m ;
210 double asym = 1.0*(N)/(N+M) ;
211 double approxErr = sqrt(4.0*n/(N+M)*(1-N/(N+M))/(N+M)) ;
212
213 asym1 = asym-nSigma*0.5*approxErr ;
214 asym2 = asym+nSigma*0.5*approxErr ;
215 return true ;
216 }
217
218 // swap n and m to ensure that n <= m
219 bool swapped(false);
220 if(n > m) {
221 swapped= true;
222 Int_t tmp(m);
223 m= n;
224 n= tmp;
225 }
226
227 // create the function objects to use
228 bool status(false);
229 BinomialSumEff upper(n,m);
230 double eff = (double)(n)/(n+m) ;
231 if(n > 0) {
232 BinomialSumEff lower(n-1,m+1);
233 status= getInterval(&upper,&lower,eff,0.1,asym1,asym2,nSigma*0.5);
234 }
235 else {
236 status= getInterval(&upper,nullptr,eff,0.1,asym1,asym2,nSigma*0.5);
237 }
238
239 // undo the swap here
240 if(swapped) {
241 double tmp(asym1);
242 asym1= 1-asym2;
243 asym2= 1-tmp;
244 }
245
246 return status;
247}
248
249
250
251////////////////////////////////////////////////////////////////////////////////
252/// Calculate a confidence interval using the cumulative functions provided.
253/// The interval will be "central" when both cumulative functions are provided,
254/// unless this would exclude the pointEstimate, in which case a one-sided interval
255/// pinned at the point estimate is returned instead.
256
257bool RooHistError::getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, double pointEstimate,
258 double stepSize, double &lo, double &hi, double nSigma) const
259{
260 // sanity checks
261 assert(nullptr != Qu || nullptr != Ql);
262
263 // convert number of sigma into a confidence level
264 double beta= TMath::Erf(nSigma/sqrt(2.));
265 double alpha= 0.5*(1-beta);
266
267 // Does the central interval contain the point estimate?
268 bool ok(true);
269 double loProb(1);
270 double hiProb(0);
271 if(nullptr != Ql) loProb= (*Ql)(&pointEstimate);
272 if(nullptr != Qu) hiProb= (*Qu)(&pointEstimate);
273
274 if (Qu && (nullptr == Ql || loProb > alpha + beta)) {
275 // Force the low edge to be at the pointEstimate
276 lo= pointEstimate;
277 double target= loProb - beta;
278 hi= seek(*Qu,lo,+stepSize,target);
279 RooBrentRootFinder uFinder(*Qu);
280 ok= uFinder.findRoot(hi,hi-stepSize,hi,target);
281 }
282 else if(Ql && (nullptr == Qu || hiProb < alpha)) {
283 // Force the high edge to be at pointEstimate
284 hi= pointEstimate;
285 double target= hiProb + beta;
286 lo= seek(*Ql,hi,-stepSize,target);
287 RooBrentRootFinder lFinder(*Ql);
288 ok= lFinder.findRoot(lo,lo,lo+stepSize,target);
289 }
290 else if (Qu && Ql) {
291 // use a central interval
292 lo= seek(*Ql,pointEstimate,-stepSize,alpha+beta);
293 hi= seek(*Qu,pointEstimate,+stepSize,alpha);
294 RooBrentRootFinder lFinder(*Ql);
295 RooBrentRootFinder uFinder(*Qu);
296 ok= lFinder.findRoot(lo,lo,lo+stepSize,alpha+beta);
297 ok|= uFinder.findRoot(hi,hi-stepSize,hi,alpha);
298 }
299 if(!ok) oocoutE(nullptr,Plotting) << "RooHistError::getInterval: failed to find root(s)" << endl;
300
301 return ok;
302}
303
304
305////////////////////////////////////////////////////////////////////////////////
306/// Scan f(x)-value until it changes sign. Start at the specified point and take constant
307/// steps of the specified size. Give up after 1000 steps.
308
309double RooHistError::seek(const RooAbsFunc &f, double startAt, double step, double value) const
310{
311 Int_t steps(1000);
312 double min(f.getMinLimit(1));
313 double max(f.getMaxLimit(1));
314 double x(startAt);
315 double 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)
#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
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
Do the root finding using the Brent-Decker method.
Singleton class used to calculate the error bars for each bin of a RooHist object.
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]
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]
static RooAbsFunc * createBinomialSum(Int_t n, Int_t m, bool eff)
Create and return a BinomialSum function binding.
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:190
TMarker m
Definition textangle.C:8