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