Logo ROOT  
Reference Guide
RooHist.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 RooHist.cxx
19\class RooHist
20\ingroup Roofitcore
21
22A RooHist is a graphical representation of binned data based on the
23TGraphAsymmErrors class. Error bars are calculated using either Poisson
24or Binomial statistics. A RooHist is used to represent histograms in
25a RooPlot.
26**/
27
28#include "RooFit.h"
29
30#include "RooHist.h"
31#include "RooHistError.h"
32#include "RooCurve.h"
33#include "RooScaledFunc.h"
34#include "RooMsgService.h"
35
36#include "TH1.h"
37#include "TClass.h"
38#include "Riostream.h"
39#include <iomanip>
40
41using namespace std;
42
44 ;
45
46
47////////////////////////////////////////////////////////////////////////////////
48/// Default constructor
49
51 _nominalBinWidth(1),
52 _nSigma(1),
53 _entries(0),
54 _rawEntries(0)
55{
56}
57
58
59
60////////////////////////////////////////////////////////////////////////////////
61/// Create an empty histogram that can be filled with the addBin()
62/// and addAsymmetryBin() methods. Use the optional parameter to
63/// specify the confidence level in units of sigma to use for
64/// calculating error bars. The nominal bin width specifies the
65/// default used by addBin(), and is used to set the relative
66/// normalization of bins with different widths.
67
68 RooHist::RooHist(Double_t nominalBinWidth, Double_t nSigma, Double_t /*xErrorFrac*/, Double_t /*scaleFactor*/) :
69 TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
70{
71 initialize();
72}
73
74
75////////////////////////////////////////////////////////////////////////////////
76/// Create a histogram from the contents of the specified TH1 object
77/// which may have fixed or variable bin widths. Error bars are
78/// calculated using Poisson statistics. Prints a warning and rounds
79/// any bins with non-integer contents. Use the optional parameter to
80/// specify the confidence level in units of sigma to use for
81/// calculating error bars. The nominal bin width specifies the
82/// default used by addBin(), and is used to set the relative
83/// normalization of bins with different widths. If not set, the
84/// nominal bin width is calculated as range/nbins.
85
86RooHist::RooHist(const TH1 &data, Double_t nominalBinWidth, Double_t nSigma, RooAbsData::ErrorType etype, Double_t xErrorFrac,
87 Bool_t correctForBinWidth, Double_t scaleFactor) :
88 TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
89{
90 initialize();
91 // copy the input histogram's name and title
92 SetName(data.GetName());
93 SetTitle(data.GetTitle());
94 // calculate our nominal bin width if necessary
95 if(_nominalBinWidth == 0) {
96 const TAxis *axis= ((TH1&)data).GetXaxis();
97 if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
98 }
100
101 // initialize our contents from the input histogram's contents
102 Int_t nbin= data.GetNbinsX();
103 for(Int_t bin= 1; bin <= nbin; bin++) {
104 Axis_t x= data.GetBinCenter(bin);
105 Stat_t y= data.GetBinContent(bin);
106 Stat_t dy = data.GetBinError(bin) ;
107 if (etype==RooAbsData::Poisson) {
108 addBin(x,y,data.GetBinWidth(bin),xErrorFrac,scaleFactor);
109 } else if (etype==RooAbsData::SumW2) {
110 addBinWithError(x,y,dy,dy,data.GetBinWidth(bin),xErrorFrac,correctForBinWidth,scaleFactor);
111 } else {
112 addBinWithError(x,y,0,0,data.GetBinWidth(bin),xErrorFrac,correctForBinWidth,scaleFactor);
113 }
114 }
115 // add over/underflow bins to our event count
116 _entries+= data.GetBinContent(0) + data.GetBinContent(nbin+1);
117}
118
119
120
121////////////////////////////////////////////////////////////////////////////////
122/// Create a histogram from the asymmetry between the specified TH1 objects
123/// which may have fixed or variable bin widths, but which must both have
124/// the same binning. The asymmetry is calculated as (1-2)/(1+2). Error bars are
125/// calculated using Binomial statistics. Prints a warning and rounds
126/// any bins with non-integer contents. Use the optional parameter to
127/// specify the confidence level in units of sigma to use for
128/// calculating error bars. The nominal bin width specifies the
129/// default used by addAsymmetryBin(), and is used to set the relative
130/// normalization of bins with different widths. If not set, the
131/// nominal bin width is calculated as range/nbins.
132
133RooHist::RooHist(const TH1 &data1, const TH1 &data2, Double_t nominalBinWidth, Double_t nSigma,
134 RooAbsData::ErrorType etype, Double_t xErrorFrac, Bool_t efficiency, Double_t scaleFactor) :
135 TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
136{
137 initialize();
138 // copy the first input histogram's name and title
139 SetName(data1.GetName());
140 SetTitle(data1.GetTitle());
141 // calculate our nominal bin width if necessary
142 if(_nominalBinWidth == 0) {
143 const TAxis *axis= ((TH1&)data1).GetXaxis();
144 if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
145 }
146
147 if (!efficiency) {
148 setYAxisLabel(Form("Asymmetry (%s - %s)/(%s + %s)",
149 data1.GetName(),data2.GetName(),data1.GetName(),data2.GetName()));
150 } else {
151 setYAxisLabel(Form("Efficiency (%s)/(%s + %s)",
152 data1.GetName(),data1.GetName(),data2.GetName()));
153 }
154 // initialize our contents from the input histogram contents
155 Int_t nbin= data1.GetNbinsX();
156 if(data2.GetNbinsX() != nbin) {
157 coutE(InputArguments) << "RooHist::RooHist: histograms have different number of bins" << endl;
158 return;
159 }
160 for(Int_t bin= 1; bin <= nbin; bin++) {
161 Axis_t x= data1.GetBinCenter(bin);
162 if(fabs(data2.GetBinCenter(bin)-x)>1e-10) {
163 coutW(InputArguments) << "RooHist::RooHist: histograms have different centers for bin " << bin << endl;
164 }
165 Stat_t y1= data1.GetBinContent(bin);
166 Stat_t y2= data2.GetBinContent(bin);
167 if (!efficiency) {
168
169 if (etype==RooAbsData::Poisson) {
170 addAsymmetryBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
171 } else if (etype==RooAbsData::SumW2) {
172 Stat_t dy1= data1.GetBinError(bin);
173 Stat_t dy2= data2.GetBinError(bin);
174 addAsymmetryBinWithError(x,y1,y2,dy1,dy2,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
175 } else {
176 addAsymmetryBinWithError(x,y1,y2,0,0,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
177 }
178
179 } else {
180
181 if (etype==RooAbsData::Poisson) {
182 addEfficiencyBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
183 } else if (etype==RooAbsData::SumW2) {
184 Stat_t dy1= data1.GetBinError(bin);
185 Stat_t dy2= data2.GetBinError(bin);
186 addEfficiencyBinWithError(x,y1,y2,dy1,dy2,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
187 } else {
188 addEfficiencyBinWithError(x,y1,y2,0,0,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
189 }
190
191 }
192
193 }
194 // we do not have a meaningful number of entries
195 _entries= -1;
196}
197
198
199
200////////////////////////////////////////////////////////////////////////////////
201/// Create histogram as sum of two existing histograms. If Poisson errors are selected the histograms are
202/// added and Poisson confidence intervals are calculated for the summed content. If wgt1 and wgt2 are not
203/// 1 in this mode, a warning message is printed. If SumW2 errors are selected the histograms are added
204/// and the histograms errors are added in quadrature, taking the weights into account.
205
206RooHist::RooHist(const RooHist& hist1, const RooHist& hist2, Double_t wgt1, Double_t wgt2,
207 RooAbsData::ErrorType etype, Double_t xErrorFrac) : _rawEntries(-1)
208{
209 // Initialize the histogram
210 initialize() ;
211
212 // Copy all non-content properties from hist1
213 SetName(hist1.GetName()) ;
214 SetTitle(hist1.GetTitle()) ;
216 _nSigma=hist1._nSigma ;
218
219 if (!hist1.hasIdenticalBinning(hist2)) {
220 coutE(InputArguments) << "RooHist::RooHist input histograms have incompatible binning, combined histogram will remain empty" << endl ;
221 return ;
222 }
223
224 if (etype==RooAbsData::Poisson) {
225 // Add histograms with Poisson errors
226
227 // Issue warning if weights are not 1
228 if (wgt1!=1.0 || wgt2 != 1.0) {
229 coutW(InputArguments) << "RooHist::RooHist: WARNING: Poisson errors of weighted sum of two histograms is not well defined! " << endl
230 << " Summed histogram bins will rounded to nearest integer for Poisson confidence interval calculation" << endl ;
231 }
232
233 // Add histograms, calculate Poisson confidence interval on sum value
234 Int_t i,n=hist1.GetN() ;
235 for(i=0 ; i<n ; i++) {
236 Double_t x1,y1,x2,y2,dx1 ;
237 hist1.GetPoint(i,x1,y1) ;
238 dx1 = hist1.GetErrorX(i) ;
239 hist2.GetPoint(i,x2,y2) ;
240 addBin(x1,roundBin(wgt1*y1+wgt2*y2),2*dx1/xErrorFrac,xErrorFrac) ;
241 }
242
243 } else {
244 // Add histograms with SumW2 errors
245
246 // Add histograms, calculate combined sum-of-weights error
247 Int_t i,n=hist1.GetN() ;
248 for(i=0 ; i<n ; i++) {
249 Double_t x1,y1,x2,y2,dx1,dy1,dy2 ;
250 hist1.GetPoint(i,x1,y1) ;
251 dx1 = hist1.GetErrorX(i) ;
252 dy1 = hist1.GetErrorY(i) ;
253 dy2 = hist2.GetErrorY(i) ;
254 hist2.GetPoint(i,x2,y2) ;
255 Double_t dy = sqrt(wgt1*wgt1*dy1*dy1+wgt2*wgt2*dy2*dy2) ;
256 addBinWithError(x1,wgt1*y1+wgt2*y2,dy,dy,2*dx1/xErrorFrac,xErrorFrac) ;
257 }
258 }
259
260}
261
262
263////////////////////////////////////////////////////////////////////////////////
264/// Create histogram from a pdf or function. Errors are computed based on the fit result provided.
265///
266/// This signature is intended for unfolding/deconvolution scenarios,
267/// where a pdf is constructed as "data minus background" and is thus
268/// intended to be displayed as "data" (or at least data-like).
269/// Usage of this signature is triggered by the draw style "P" in RooAbsReal::plotOn.
270///
271/// More details.
272/// \param[in] f The function to be plotted.
273/// \param[in] x The variable on the x-axis
274/// \param[in] xErrorFrac Size of the errror in x as a fraction of the bin width
275/// \param[in] scaleFactor arbitrary scaling of the y-values
276/// \param[in] normVars variables over which to normalize
277RooHist::RooHist(const RooAbsReal &f, RooAbsRealLValue &x, Double_t xErrorFrac, Double_t scaleFactor, const RooArgSet *normVars, const RooFitResult* fr) :
278 TGraphAsymmErrors(), _nSigma(1), _rawEntries(-1)
279{
280 // grab the function's name and title
281 TString name(f.GetName());
282 SetName(name.Data());
283 TString title(f.GetTitle());
284 SetTitle(title.Data());
285 // append " ( [<funit> ][/ <xunit> ])" to our y-axis label if necessary
286 if(0 != strlen(f.getUnit()) || 0 != strlen(x.getUnit())) {
287 title.Append(" ( ");
288 if(0 != strlen(f.getUnit())) {
289 title.Append(f.getUnit());
290 title.Append(" ");
291 }
292 if(0 != strlen(x.getUnit())) {
293 title.Append("/ ");
294 title.Append(x.getUnit());
295 title.Append(" ");
296 }
297 title.Append(")");
298 }
299 setYAxisLabel(title.Data());
300
301 RooAbsFunc *funcPtr = nullptr;
302 RooAbsFunc *rawPtr = nullptr;
303 funcPtr= f.bindVars(x,normVars,kTRUE);
304
305 // apply a scale factor if necessary
306 if(scaleFactor != 1) {
307 rawPtr= funcPtr;
308 funcPtr= new RooScaledFunc(*rawPtr,scaleFactor);
309 }
310
311 // apply a scale factor if necessary
312 assert(funcPtr);
313
314 // calculate the points to add to our curve
315 int xbins = x.numBins();
316 RooArgSet nset;
317 if(normVars) nset.add(*normVars);
318 for(int i=0; i<xbins; ++i){
319 double xval = x.getBinning().binCenter(i);
320 double xwidth = x.getBinning().binWidth(i);
321 Axis_t xval_ax = xval;
322 double yval = (*funcPtr)(&xval);
323 double yerr = sqrt(yval);
324 if(fr) yerr = f.getPropagatedError(*fr,nset);
325 addBinWithError(xval_ax,yval,yerr,yerr,xwidth,xErrorFrac,false,scaleFactor) ;
326 _entries += yval;
327 }
328 _nominalBinWidth = 1.;
329
330 // cleanup
331 delete funcPtr;
332 if(rawPtr) delete rawPtr;
333}
334
335
336////////////////////////////////////////////////////////////////////////////////
337/// Perform common initialization for all constructors.
338
340{
342 _entries= 0;
343}
344
345
346////////////////////////////////////////////////////////////////////////////////
347/// Return the number of events of the dataset associated with this RooHist.
348/// This is the number of events in the RooHist itself, unless a different
349/// value was specified through setRawEntries()
350
352{
353 return (_rawEntries==-1 ? _entries : _rawEntries) ;
354}
355
356
357////////////////////////////////////////////////////////////////////////////////
358/// Calculate integral of histogram in given range
359
361{
362 Double_t sum(0) ;
363 for (int i=0 ; i<GetN() ; i++) {
364 Double_t x,y ;
365
366 GetPoint(i,x,y) ;
367
368 if (x>=xlo && x<=xhi) {
369 sum += y ;
370 }
371 }
372
373 if (_rawEntries!=-1) {
374 coutW(Plotting) << "RooHist::getFitRangeNEvt() WARNING: The number of normalisation events associated to histogram " << GetName() << " is not equal to number of events in this histogram."
375 << "\n\t\t This is due a cut being applied while plotting the data. Automatic normalisation over a sub-range of a plot variable assumes"
376 << "\n\t\t that the effect of that cut is uniform across the plot, which may be an incorrect assumption. To obtain a correct normalisation, it needs to be passed explicitly:"
377 << "\n\t\t\t data->plotOn(frame01,CutRange(\"SB1\"));"
378 << "\n\t\t\t const double nData = data->sumEntries(\"\", \"SB1\"); //or the cut string such as sumEntries(\"x > 0.\");"
379 << "\n\t\t\t model.plotOn(frame01, RooFit::Normalization(nData, RooAbsReal::NumEvent), ProjectionRange(\"SB1\"));" << endl ;
381 }
382
383 return sum ;
384}
385
386
387
388////////////////////////////////////////////////////////////////////////////////
389/// Return (average) bin width of this RooHist
390
392{
393 return _nominalBinWidth ;
394}
395
396
397
398////////////////////////////////////////////////////////////////////////////////
399/// Return the nearest positive integer to the input value
400/// and print a warning if an adjustment is required.
401
403{
404 if(y < 0) {
405 coutW(Plotting) << fName << "::roundBin: rounding negative bin contents to zero: " << y << endl;
406 return 0;
407 }
408 Int_t n= (Int_t)(y+0.5);
409 if(fabs(y-n)>1e-6) {
410 coutW(Plotting) << fName << "::roundBin: rounding non-integer bin contents: " << y << endl;
411 }
412 return n;
413}
414
415
416
417////////////////////////////////////////////////////////////////////////////////
418/// Add a bin to this histogram with the specified integer bin contents
419/// and using an error bar calculated with Poisson statistics. The bin width
420/// is used to set the relative scale of bins with different widths.
421
422void RooHist::addBin(Axis_t binCenter, Double_t n, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
423{
424 if (n<0) {
425 coutW(Plotting) << "RooHist::addBin(" << GetName() << ") WARNING: negative entry set to zero when Poisson error bars are requested" << endl ;
426 }
427
428 Double_t scale= 1;
429 if(binWidth > 0) {
430 scale= _nominalBinWidth/binWidth;
431 }
432 _entries+= n;
433 Int_t index= GetN();
434
435 // calculate Poisson errors for this bin
436 Double_t ym,yp,dx(0.5*binWidth);
437
438 if (fabs((double)((n-Int_t(n))>1e-5))) {
439 // need interpolation
440 Double_t ym1(0),yp1(0),ym2(0),yp2(0) ;
441 Int_t n1 = Int_t(n) ;
442 Int_t n2 = n1+1 ;
443 if(!RooHistError::instance().getPoissonInterval(n1,ym1,yp1,_nSigma) ||
444 !RooHistError::instance().getPoissonInterval(n2,ym2,yp2,_nSigma)) {
445 coutE(Plotting) << "RooHist::addBin: unable to add bin with " << n << " events" << endl;
446 }
447 ym = ym1 + (n-n1)*(ym2-ym1) ;
448 yp = yp1 + (n-n1)*(yp2-yp1) ;
449 coutW(Plotting) << "RooHist::addBin(" << GetName()
450 << ") WARNING: non-integer bin entry " << n << " with Poisson errors, interpolating between Poisson errors of adjacent integer" << endl ;
451 } else {
452 // integer case
453 if(!RooHistError::instance().getPoissonInterval(Int_t(n),ym,yp,_nSigma)) {
454 coutE(Plotting) << "RooHist::addBin: unable to add bin with " << n << " events" << endl;
455 return;
456 }
457 }
458
459 SetPoint(index,binCenter,n*scale*scaleFactor);
460 SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,scale*(n-ym)*scaleFactor,scale*(yp-n)*scaleFactor);
461 updateYAxisLimits(scale*yp);
462 updateYAxisLimits(scale*ym);
463}
464
465
466
467////////////////////////////////////////////////////////////////////////////////
468/// Add a bin to this histogram with the specified bin contents
469/// and error. The bin width is used to set the relative scale of
470/// bins with different widths.
471
472void RooHist::addBinWithError(Axis_t binCenter, Double_t n, Double_t elow, Double_t ehigh, Double_t binWidth,
473 Double_t xErrorFrac, Bool_t correctForBinWidth, Double_t scaleFactor)
474{
475 Double_t scale= 1;
476 if(binWidth > 0 && correctForBinWidth) {
477 scale= _nominalBinWidth/binWidth;
478 }
479 _entries+= n;
480 Int_t index= GetN();
481
482 Double_t dx(0.5*binWidth) ;
483 SetPoint(index,binCenter,n*scale*scaleFactor);
484 SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,elow*scale*scaleFactor,ehigh*scale*scaleFactor);
485 updateYAxisLimits(scale*(n-elow));
486 updateYAxisLimits(scale*(n+ehigh));
487}
488
489
490
491
492////////////////////////////////////////////////////////////////////////////////
493/// Add a bin to this histogram with the specified bin contents
494/// and error. The bin width is used to set the relative scale of
495/// bins with different widths.
496
497void RooHist::addBinWithXYError(Axis_t binCenter, Double_t n, Double_t exlow, Double_t exhigh, Double_t eylow, Double_t eyhigh,
498 Double_t scaleFactor)
499{
500 _entries+= n;
501 Int_t index= GetN();
502
503 SetPoint(index,binCenter,n*scaleFactor);
504 SetPointError(index,exlow,exhigh,eylow*scaleFactor,eyhigh*scaleFactor);
505 updateYAxisLimits(scaleFactor*(n-eylow));
506 updateYAxisLimits(scaleFactor*(n+eyhigh));
507}
508
509
510
511
512
513////////////////////////////////////////////////////////////////////////////////
514/// Add a bin to this histogram with the value (n1-n2)/(n1+n2)
515/// using an error bar calculated with Binomial statistics.
516
517void RooHist::addAsymmetryBin(Axis_t binCenter, Int_t n1, Int_t n2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
518{
519 Double_t scale= 1;
520 if(binWidth > 0) scale= _nominalBinWidth/binWidth;
521 Int_t index= GetN();
522
523 // calculate Binomial errors for this bin
524 Double_t ym,yp,dx(0.5*binWidth);
525 if(!RooHistError::instance().getBinomialIntervalAsym(n1,n2,ym,yp,_nSigma)) {
526 coutE(Plotting) << "RooHist::addAsymmetryBin: unable to calculate binomial error for bin with " << n1 << "," << n2 << " events" << endl;
527 return;
528 }
529
530 Double_t a= (Double_t)(n1-n2)/(n1+n2);
531 SetPoint(index,binCenter,a*scaleFactor);
532 SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
533 updateYAxisLimits(scale*yp);
534 updateYAxisLimits(scale*ym);
535}
536
537
538
539////////////////////////////////////////////////////////////////////////////////
540/// Add a bin to this histogram with the value (n1-n2)/(n1+n2)
541/// using an error bar calculated with Binomial statistics.
542
543void RooHist::addAsymmetryBinWithError(Axis_t binCenter, Double_t n1, Double_t n2, Double_t en1, Double_t en2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
544{
545 Double_t scale= 1;
546 if(binWidth > 0) scale= _nominalBinWidth/binWidth;
547 Int_t index= GetN();
548
549 // calculate Binomial errors for this bin
550 Double_t ym,yp,dx(0.5*binWidth);
551 Double_t a= (Double_t)(n1-n2)/(n1+n2);
552
553 Double_t error = 2*sqrt( pow(en1,2)*pow(n2,2) + pow(en2,2)*pow(n1,2) ) / pow(n1+n2,2) ;
554 ym=a-error ;
555 yp=a+error ;
556
557 SetPoint(index,binCenter,a*scaleFactor);
558 SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
559 updateYAxisLimits(scale*yp);
560 updateYAxisLimits(scale*ym);
561}
562
563
564
565////////////////////////////////////////////////////////////////////////////////
566/// Add a bin to this histogram with the value n1/(n1+n2)
567/// using an error bar calculated with Binomial statistics.
568
569void RooHist::addEfficiencyBin(Axis_t binCenter, Int_t n1, Int_t n2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
570{
571 Double_t scale= 1;
572 if(binWidth > 0) scale= _nominalBinWidth/binWidth;
573 Int_t index= GetN();
574
575 Double_t a= (Double_t)(n1)/(n1+n2);
576
577 // calculate Binomial errors for this bin
578 Double_t ym,yp,dx(0.5*binWidth);
579 if(!RooHistError::instance().getBinomialIntervalEff(n1,n2,ym,yp,_nSigma)) {
580 coutE(Plotting) << "RooHist::addEfficiencyBin: unable to calculate binomial error for bin with " << n1 << "," << n2 << " events" << endl;
581 return;
582 }
583
584 SetPoint(index,binCenter,a*scaleFactor);
585 SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
586 updateYAxisLimits(scale*yp);
587 updateYAxisLimits(scale*ym);
588}
589
590
591
592////////////////////////////////////////////////////////////////////////////////
593/// Add a bin to this histogram with the value n1/(n1+n2)
594/// using an error bar calculated with Binomial statistics.
595
596void RooHist::addEfficiencyBinWithError(Axis_t binCenter, Double_t n1, Double_t n2, Double_t en1, Double_t en2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
597{
598 Double_t scale= 1;
599 if(binWidth > 0) scale= _nominalBinWidth/binWidth;
600 Int_t index= GetN();
601
602 Double_t a= (Double_t)(n1)/(n1+n2);
603
604 Double_t error = sqrt( pow(en1,2)*pow(n2,2) + pow(en2,2)*pow(n1,2) ) / pow(n1+n2,2) ;
605
606 // calculate Binomial errors for this bin
607 Double_t ym,yp,dx(0.5*binWidth);
608 ym=a-error ;
609 yp=a+error ;
610
611
612 SetPoint(index,binCenter,a*scaleFactor);
613 SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
614 updateYAxisLimits(scale*yp);
615 updateYAxisLimits(scale*ym);
616}
617
618
619
620////////////////////////////////////////////////////////////////////////////////
621/// Destructor
622
624{
625}
626
627
628
629////////////////////////////////////////////////////////////////////////////////
630/// Return kTRUE if binning of this RooHist is identical to that of 'other'
631
633{
634 // First check if number of bins is the same
635 if (GetN() != other.GetN()) {
636 return kFALSE ;
637 }
638
639 // Next require that all bin centers are the same
640 Int_t i ;
641 for (i=0 ; i<GetN() ; i++) {
642 Double_t x1,x2,y1,y2 ;
643
644 GetPoint(i,x1,y1) ;
645 other.GetPoint(i,x2,y2) ;
646
647 if (fabs(x1-x2)>1e-10) {
648 return kFALSE ;
649 }
650
651 }
652
653 return kTRUE ;
654}
655
656
657
658////////////////////////////////////////////////////////////////////////////////
659/// Return kTRUE if contents of this RooHist is identical within given
660/// relative tolerance to that of 'other'
661
663{
664 // Make temporary TH1s output of RooHists to perform Kolmogorov test
666 TH1F h_self("h_self","h_self",GetN(),0,1) ;
667 TH1F h_other("h_other","h_other",GetN(),0,1) ;
669
670 for (Int_t i=0 ; i<GetN() ; i++) {
671 h_self.SetBinContent(i+1,GetY()[i]) ;
672 h_other.SetBinContent(i+1,other.GetY()[i]) ;
673 }
674
675 Double_t M = h_self.KolmogorovTest(&h_other,"M") ;
676 if (M>tol) {
677 Double_t kprob = h_self.KolmogorovTest(&h_other) ;
678 cout << "RooHist::isIdentical() tolerance exceeded M=" << M << " (tol=" << tol << "), corresponding prob = " << kprob << endl ;
679 return kFALSE ;
680 }
681
682 return kTRUE ;
683}
684
685
686
687////////////////////////////////////////////////////////////////////////////////
688/// Print info about this histogram to the specified output stream.
689///
690/// Standard: number of entries
691/// Shape: error CL and maximum value
692/// Verbose: print our bin contents and errors
693
694void RooHist::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const
695{
697 os << indent << "--- RooHist ---" << endl;
698 Int_t n= GetN();
699 os << indent << " Contains " << n << " bins" << endl;
700 if(verbose) {
701 os << indent << " Errors calculated at" << _nSigma << "-sigma CL" << endl;
702 os << indent << " Bin Contents:" << endl;
703 for(Int_t i= 0; i < n; i++) {
704 os << indent << setw(3) << i << ") x= " << fX[i];
705 if(fEXhigh[i] > 0 || fEXlow[i] > 0) {
706 os << " +" << fEXhigh[i] << " -" << fEXlow[i];
707 }
708 os << " , y = " << fY[i] << " +" << fEYhigh[i] << " -" << fEYlow[i] << endl;
709 }
710 }
711}
712
713
714
715////////////////////////////////////////////////////////////////////////////////
716/// Print name of RooHist
717
718void RooHist::printName(ostream& os) const
719{
720 os << GetName() ;
721}
722
723
724
725////////////////////////////////////////////////////////////////////////////////
726/// Print title of RooHist
727
728void RooHist::printTitle(ostream& os) const
729{
730 os << GetTitle() ;
731}
732
733
734
735////////////////////////////////////////////////////////////////////////////////
736/// Print class name of RooHist
737
738void RooHist::printClassName(ostream& os) const
739{
740 os << IsA()->GetName() ;
741}
742
743
744
745////////////////////////////////////////////////////////////////////////////////
746/// Create and return RooHist containing residuals w.r.t to given curve.
747/// If normalize is true, the residuals are normalized by the histogram
748/// errors creating a RooHist with pull values
749
750RooHist* RooHist::makeResidHist(const RooCurve& curve, bool normalize, bool useAverage) const
751{
752
753 // Copy all non-content properties from hist1
754 RooHist* hist = new RooHist(_nominalBinWidth) ;
755 if (normalize) {
756 hist->SetName(Form("pull_%s_%s",GetName(),curve.GetName())) ;
757 hist->SetTitle(Form("Pull of %s and %s",GetTitle(),curve.GetTitle())) ;
758 } else {
759 hist->SetName(Form("resid_%s_%s",GetName(),curve.GetName())) ;
760 hist->SetTitle(Form("Residual of %s and %s",GetTitle(),curve.GetTitle())) ;
761 }
762
763 // Determine range of curve
764 Double_t xstart,xstop,y ;
765 curve.GetPoint(0,xstart,y) ;
766 curve.GetPoint(curve.GetN()-1,xstop,y) ;
767
768 // Add histograms, calculate Poisson confidence interval on sum value
769 for(Int_t i=0 ; i<GetN() ; i++) {
770 Double_t x,point;
771 GetPoint(i,x,point) ;
772
773 // Only calculate pull for bins inside curve range
774 if (x<xstart || x>xstop) continue ;
775
776 Double_t yy ;
777 if (useAverage) {
778 Double_t exl = GetErrorXlow(i);
779 Double_t exh = GetErrorXhigh(i) ;
780 if (exl<=0 ) exl = GetErrorX(i);
781 if (exh<=0 ) exh = GetErrorX(i);
782 if (exl<=0 ) exl = 0.5*getNominalBinWidth();
783 if (exh<=0 ) exh = 0.5*getNominalBinWidth();
784 yy = point - curve.average(x-exl,x+exh) ;
785 } else {
786 yy = point - curve.interpolate(x) ;
787 }
788
789 Double_t dyl = GetErrorYlow(i) ;
790 Double_t dyh = GetErrorYhigh(i) ;
791 if (normalize) {
792 Double_t norm = (yy>0?dyl:dyh);
793 if (norm==0.) {
794 coutW(Plotting) << "RooHist::makeResisHist(" << GetName() << ") WARNING: point " << i << " has zero error, setting residual to zero" << endl ;
795 yy=0 ;
796 dyh=0 ;
797 dyl=0 ;
798 } else {
799 yy /= norm;
800 dyh /= norm;
801 dyl /= norm;
802 }
803 }
804 hist->addBinWithError(x,yy,dyl,dyh);
805 }
806 return hist ;
807}
#define f(i)
Definition: RSha256.hxx:104
#define e(i)
Definition: RSha256.hxx:103
static const double x2[5]
static const double x1[5]
#define coutW(a)
Definition: RooMsgService.h:32
#define coutE(a)
Definition: RooMsgService.h:33
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Axis_t
Definition: RtypesCore.h:74
double Double_t
Definition: RtypesCore.h:57
double Stat_t
Definition: RtypesCore.h:75
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
static void indent(ostringstream &buf, int indent_level)
char name[80]
Definition: TGX11.cxx:109
double pow(double, double)
double sqrt(double)
char * Form(const char *fmt,...)
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
A RooCurve is a one-dimensional graphical representation of a real-valued function.
Definition: RooCurve.h:32
Double_t average(Double_t lo, Double_t hi) const
Return average curve value in [xFirst,xLast] by integrating curve between points and dividing by xLas...
Definition: RooCurve.cxx:612
Double_t interpolate(Double_t x, Double_t tolerance=1e-10) const
Return linearly interpolated value of curve at xvalue.
Definition: RooCurve.cxx:698
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
static const RooHistError & instance()
Return a reference to a singleton object that is created the first time this method is called.
A RooHist is a graphical representation of binned data based on the TGraphAsymmErrors class.
Definition: RooHist.h:27
void addAsymmetryBinWithError(Axis_t binCenter, Double_t n1, Double_t n2, Double_t en1, Double_t en2, Double_t binWidth=0, Double_t xErrorFrac=1.0, Double_t scaleFactor=1.0)
Add a bin to this histogram with the value (n1-n2)/(n1+n2) using an error bar calculated with Binomia...
Definition: RooHist.cxx:543
RooHist * makeResidHist(const RooCurve &curve, bool normalize=false, bool useAverage=false) const
Create and return RooHist containing residuals w.r.t to given curve.
Definition: RooHist.cxx:750
Double_t _nominalBinWidth
Definition: RooHist.h:89
void addEfficiencyBin(Axis_t binCenter, Int_t n1, Int_t n2, Double_t binWidth=0, Double_t xErrorFrac=1.0, Double_t scaleFactor=1.0)
Add a bin to this histogram with the value n1/(n1+n2) using an error bar calculated with Binomial sta...
Definition: RooHist.cxx:569
Double_t getFitRangeNEvt() const
Return the number of events of the dataset associated with this RooHist.
Definition: RooHist.cxx:351
Double_t _rawEntries
Definition: RooHist.h:92
virtual void printTitle(std::ostream &os) const
Print title of RooHist.
Definition: RooHist.cxx:728
Int_t roundBin(Double_t y)
Return the nearest positive integer to the input value and print a warning if an adjustment is requir...
Definition: RooHist.cxx:402
void initialize()
Perform common initialization for all constructors.
Definition: RooHist.cxx:339
RooHist()
Default constructor.
Definition: RooHist.cxx:50
virtual void printName(std::ostream &os) const
Print name of RooHist.
Definition: RooHist.cxx:718
void addBinWithError(Axis_t binCenter, Double_t n, Double_t elow, Double_t ehigh, Double_t binWidth=0, Double_t xErrorFrac=1.0, Bool_t correctForBinWidth=kTRUE, Double_t scaleFactor=1.0)
Add a bin to this histogram with the specified bin contents and error.
Definition: RooHist.cxx:472
virtual void printMultiline(std::ostream &os, Int_t content, Bool_t verbose=kFALSE, TString indent="") const
Print info about this histogram to the specified output stream.
Definition: RooHist.cxx:694
void addEfficiencyBinWithError(Axis_t binCenter, Double_t n1, Double_t n2, Double_t en1, Double_t en2, Double_t binWidth=0, Double_t xErrorFrac=1.0, Double_t scaleFactor=1.0)
Add a bin to this histogram with the value n1/(n1+n2) using an error bar calculated with Binomial sta...
Definition: RooHist.cxx:596
void addBin(Axis_t binCenter, Double_t n, Double_t binWidth=0, Double_t xErrorFrac=1.0, Double_t scaleFactor=1.0)
Add a bin to this histogram with the specified integer bin contents and using an error bar calculated...
Definition: RooHist.cxx:422
Double_t getNominalBinWidth() const
Definition: RooHist.h:71
virtual void printClassName(std::ostream &os) const
Print class name of RooHist.
Definition: RooHist.cxx:738
Double_t _entries
Definition: RooHist.h:91
Double_t getFitRangeBinW() const
Return (average) bin width of this RooHist.
Definition: RooHist.cxx:391
void addAsymmetryBin(Axis_t binCenter, Int_t n1, Int_t n2, Double_t binWidth=0, Double_t xErrorFrac=1.0, Double_t scaleFactor=1.0)
Add a bin to this histogram with the value (n1-n2)/(n1+n2) using an error bar calculated with Binomia...
Definition: RooHist.cxx:517
Bool_t hasIdenticalBinning(const RooHist &other) const
Return kTRUE if binning of this RooHist is identical to that of 'other'.
Definition: RooHist.cxx:632
virtual ~RooHist()
Destructor.
Definition: RooHist.cxx:623
void addBinWithXYError(Axis_t binCenter, Double_t n, Double_t exlow, Double_t exhigh, Double_t eylow, Double_t eyhigh, Double_t scaleFactor=1.0)
Add a bin to this histogram with the specified bin contents and error.
Definition: RooHist.cxx:497
Double_t _nSigma
Definition: RooHist.h:90
Bool_t isIdentical(const RooHist &other, Double_t tol=1e-6) const
Return kTRUE if contents of this RooHist is identical within given relative tolerance to that of 'oth...
Definition: RooHist.cxx:662
void updateYAxisLimits(Double_t y)
Definition: RooPlotable.h:33
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Print detailed information.
Definition: RooPlotable.cxx:42
void setYAxisLabel(const char *label)
Definition: RooPlotable.h:32
const char * getYAxisLabel() const
Definition: RooPlotable.h:31
Lightweight RooAbsFunction implementation that applies a constant scale factor to another RooAbsFunc.
Definition: RooScaledFunc.h:21
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
Class to manage histogram axis.
Definition: TAxis.h:30
Double_t GetXmax() const
Definition: TAxis.h:134
Double_t GetXmin() const
Definition: TAxis.h:133
Int_t GetNbins() const
Definition: TAxis.h:121
const char * GetTitle() const
Returns title of object.
Definition: TAxis.h:129
TGraph with asymmetric error bars.
Double_t * fEXhigh
[fNpoints] array of X high errors
virtual void SetPointError(Double_t exl, Double_t exh, Double_t eyl, Double_t eyh)
Set ex and ey values for point pointed by the mouse.
Double_t GetErrorYhigh(Int_t i) const
Get high error on Y.
Double_t * fEYhigh
[fNpoints] array of Y high errors
Double_t GetErrorYlow(Int_t i) const
Get low error on Y.
Double_t GetErrorXlow(Int_t i) const
Get low error on X.
Double_t * fEYlow
[fNpoints] array of Y low errors
Double_t GetErrorXhigh(Int_t i) const
Get high error on X.
Double_t * fEXlow
[fNpoints] array of X low errors
Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorX(Int_t bin) const
This function is called by GraphFitChisquare.
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2269
Double_t * GetY() const
Definition: TGraph.h:131
virtual void SetName(const char *name="")
Set graph name.
Definition: TGraph.cxx:2308
Int_t GetN() const
Definition: TGraph.h:123
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2324
Double_t * fY
[fNpoints] array of Y points
Definition: TGraph.h:48
Double_t * fX
[fNpoints] array of X points
Definition: TGraph.h:47
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
Get x and y values for point number i.
Definition: TGraph.cxx:1593
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
The TH1 histogram class.
Definition: TH1.h:56
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8597
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1226
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
TAxis * GetYaxis()
Definition: TH1.h:317
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8678
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:8619
virtual Double_t KolmogorovTest(const TH1 *h2, Option_t *option="") const
Statistical test of compatibility in shape between this histogram and h2, using Kolmogorov test.
Definition: TH1.cxx:7684
TString fName
Definition: TNamed.h:32
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
TString & Append(const char *cs)
Definition: TString.h:559
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
@ InputArguments
Definition: RooGlobalFunc.h:68
auto * a
Definition: textangle.C:12
static long int sum(long int i)
Definition: Factory.cxx:2275