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