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