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