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
22Graphical 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;
224 Int_t n = hist1.GetN();
225 for(i=0 ; i<n ; i++) {
226 double x1;
227 double y1;
228 double x2;
229 double y2;
230 double dx1;
231 hist1.GetPoint(i,x1,y1) ;
232 dx1 = hist1.GetErrorX(i) ;
233 hist2.GetPoint(i,x2,y2) ;
234 addBin(x1,roundBin(wgt1*y1+wgt2*y2),2*dx1/xErrorFrac,xErrorFrac) ;
235 }
236
237 } else {
238 // Add histograms with SumW2 errors
239
240 // Add histograms, calculate combined sum-of-weights error
241 Int_t i;
242 Int_t n = hist1.GetN();
243 for(i=0 ; i<n ; i++) {
244 double x1;
245 double y1;
246 double x2;
247 double y2;
248 double dx1;
249 double dy1;
250 double dy2;
251 hist1.GetPoint(i,x1,y1) ;
252 dx1 = hist1.GetErrorX(i) ;
253 dy1 = hist1.GetErrorY(i) ;
254 dy2 = hist2.GetErrorY(i) ;
255 hist2.GetPoint(i,x2,y2) ;
256 double dy = sqrt(wgt1*wgt1*dy1*dy1+wgt2*wgt2*dy2*dy2) ;
257 addBinWithError(x1,wgt1*y1+wgt2*y2,dy,dy,2*dx1/xErrorFrac,xErrorFrac) ;
258 }
259 }
260
261}
262
263
264////////////////////////////////////////////////////////////////////////////////
265/// Create histogram from a pdf or function. Errors are computed based on the fit result provided.
266///
267/// This signature is intended for unfolding/deconvolution scenarios,
268/// where a pdf is constructed as "data minus background" and is thus
269/// intended to be displayed as "data" (or at least data-like).
270/// Usage of this signature is triggered by the draw style "P" in RooAbsReal::plotOn.
271///
272/// More details.
273/// \param[in] f The function to be plotted.
274/// \param[in] x The variable on the x-axis
275/// \param[in] xErrorFrac Size of the error in x as a fraction of the bin width
276/// \param[in] scaleFactor arbitrary scaling of the y-values
277/// \param[in] normVars variables over which to normalize
278/// \param[in] fr fit result
279RooHist::RooHist(const RooAbsReal &f, RooAbsRealLValue &x, double xErrorFrac, double scaleFactor,
280 const RooArgSet *normVars, const RooFitResult *fr)
281 : _nSigma(1), _rawEntries(-1)
282{
283 // grab the function's name and title
284 SetName(f.GetName());
285 std::string title{f.GetTitle()};
286 SetTitle(title.c_str());
287 // append " ( [<funit> ][/ <xunit> ])" to our y-axis label if necessary
288 if(0 != strlen(f.getUnit()) || 0 != strlen(x.getUnit())) {
289 title += " ( ";
290 if(0 != strlen(f.getUnit())) {
291 title += f.getUnit();
292 title += " ";
293 }
294 if(0 != strlen(x.getUnit())) {
295 title += "/ ";
296 title += x.getUnit();
297 title += " ";
298 }
299 title += ")";
300 }
301 setYAxisLabel(title.c_str());
302
303 RooProduct scaledFunc{"scaled_func", "scaled_func", {f, RooFit::RooConst(scaleFactor)}};
304 std::unique_ptr<RooAbsFunc> funcPtr{scaledFunc.bindVars(x, normVars, true)};
305
306 // calculate the points to add to our curve
307 int xbins = x.numBins();
308 RooArgSet nset;
309 if(normVars) nset.add(*normVars);
310 for(int i=0; i<xbins; ++i){
311 double xval = x.getBinning().binCenter(i);
312 double xwidth = x.getBinning().binWidth(i);
313 Axis_t xval_ax = xval;
314 double yval = (*funcPtr)(&xval);
315 double yerr = std::sqrt(yval);
316 if(fr) yerr = f.getPropagatedError(*fr,nset);
317 addBinWithError(xval_ax,yval,yerr,yerr,xwidth,xErrorFrac,false,scaleFactor) ;
318 _entries += yval;
319 }
320 _nominalBinWidth = 1.;
321}
322
323
324////////////////////////////////////////////////////////////////////////////////
325/// Perform common initialization for all constructors.
326
328{
330}
331
332
333////////////////////////////////////////////////////////////////////////////////
334/// Return the number of events of the dataset associated with this RooHist.
335/// This is the number of events in the RooHist itself, unless a different
336/// value was specified through setRawEntries()
337
339{
340 return (_rawEntries==-1 ? _entries : _rawEntries) ;
341}
342
343
344////////////////////////////////////////////////////////////////////////////////
345/// Calculate integral of histogram in given range
346
347double RooHist::getFitRangeNEvt(double xlo, double xhi) const
348{
349 double sum(0) ;
350 for (int i=0 ; i<GetN() ; i++) {
351 double x;
352 double y;
353
354 GetPoint(i,x,y) ;
355
356 if (x>=xlo && x<=xhi) {
357 // We have to use the original weights of the histogram, because the
358 // scaled points have nothing to do anymore with event weights in the
359 // case of non-uniform binning. For backwards compatibility with the
360 // RooHist version 1, we first need to check if the `_originalWeights`
361 // member is filled.
362 sum += _originalWeights.empty() ? y : _originalWeights[i];
363 }
364 }
365
366 if (_rawEntries!=-1) {
367 coutW(Plotting) << "RooHist::getFitRangeNEvt() WARNING: The number of normalisation events associated to histogram " << GetName() << " is not equal to number of events in this histogram."
368 << "\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"
369 << "\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:"
370 << "\n\t\t\t data->plotOn(frame01,CutRange(\"SB1\"));"
371 << "\n\t\t\t const double nData = data->sumEntries(\"\", \"SB1\"); //or the cut string such as sumEntries(\"x > 0.\");"
372 << "\n\t\t\t model.plotOn(frame01, RooFit::Normalization(nData, RooAbsReal::NumEvent), ProjectionRange(\"SB1\"));" << std::endl ;
374 }
375
376 return sum ;
377}
378
379
380////////////////////////////////////////////////////////////////////////////////
381/// Return the nearest positive integer to the input value
382/// and print a warning if an adjustment is required.
383
385{
386 if(y < 0) {
387 coutW(Plotting) << fName << "::roundBin: rounding negative bin contents to zero: " << y << std::endl;
388 return 0;
389 }
390 Int_t n= (Int_t)(y+0.5);
391 if(std::abs(y-n)>1e-6) {
392 coutW(Plotting) << fName << "::roundBin: rounding non-integer bin contents: " << y << std::endl;
393 }
394 return n;
395}
396
397
398void RooHist::addPoint(Axis_t binCenter, double y, double yscale, double exlow, double exhigh, double eylow, double eyhigh)
399{
400 const int index = GetN();
401 SetPoint(index, binCenter, y*yscale);
402
403 // If the scale is negative, the low and high errors must be swapped
404 if(std::abs(yscale) < 0) {
405 std::swap(eylow, eyhigh);
406 }
407
408 SetPointError(index, exlow, exhigh, std::abs(yscale) * eylow, std::abs(yscale) * eyhigh);
409
410 updateYAxisLimits(yscale * (y - eylow));
411 updateYAxisLimits(yscale * (y + eyhigh));
412
413 // We also track the original weights of the histogram, because if we only
414 // have info on the scaled points it's not possible anymore to compute the
415 // number of events in a subrange of the RooHist.
416 _originalWeights.resize(index + 1);
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 n, double binWidth, double xErrorFrac, double scaleFactor)
427{
428 if (n<0) {
429 coutW(Plotting) << "RooHist::addBin(" << GetName() << ") WARNING: negative entry set to zero when Poisson error bars are requested" << std::endl ;
430 }
431
432 double scale= 1;
433 if(binWidth > 0) {
434 scale= _nominalBinWidth/binWidth;
435 }
436 _entries+= n;
437
438 // calculate Poisson errors for this bin
439 double ym;
440 double yp;
441 double dx(0.5 * binWidth);
442
443 if (std::abs((double)((n-Int_t(n))>1e-5))) {
444 // need interpolation
445 double ym1(0);
446 double yp1(0);
447 double ym2(0);
448 double yp2(0);
449 Int_t n1 = Int_t(n) ;
450 Int_t n2 = n1+1 ;
451 if(!RooHistError::instance().getPoissonInterval(n1,ym1,yp1,_nSigma) ||
452 !RooHistError::instance().getPoissonInterval(n2,ym2,yp2,_nSigma)) {
453 coutE(Plotting) << "RooHist::addBin: unable to add bin with " << n << " events" << std::endl;
454 }
455 ym = ym1 + (n-n1)*(ym2-ym1) ;
456 yp = yp1 + (n-n1)*(yp2-yp1) ;
457 coutW(Plotting) << "RooHist::addBin(" << GetName()
458 << ") WARNING: non-integer bin entry " << n << " with Poisson errors, interpolating between Poisson errors of adjacent integer" << std::endl ;
459 } else {
460 // integer case
461 if(!RooHistError::instance().getPoissonInterval(Int_t(n),ym,yp,_nSigma)) {
462 coutE(Plotting) << "RooHist::addBin: unable to add bin with " << n << " events" << std::endl;
463 return;
464 }
465 }
466
467 addPoint(binCenter,n, scale*scaleFactor,dx*xErrorFrac,dx*xErrorFrac, n-ym, yp-n);
468}
469
470
471
472////////////////////////////////////////////////////////////////////////////////
473/// Add a bin to this histogram with the specified bin contents
474/// and error. The bin width is used to set the relative scale of
475/// bins with different widths.
476
477void RooHist::addBinWithError(Axis_t binCenter, double n, double elow, double ehigh, double binWidth,
478 double xErrorFrac, bool correctForBinWidth, double scaleFactor)
479{
480 double scale= 1;
481 if(binWidth > 0 && correctForBinWidth) {
482 scale= _nominalBinWidth/binWidth;
483 }
484 _entries+= n;
485
486 double dx(0.5*binWidth) ;
487 addPoint(binCenter,n, scale*scaleFactor,dx*xErrorFrac,dx*xErrorFrac, elow, ehigh);
488}
489
490
491
492
493////////////////////////////////////////////////////////////////////////////////
494/// Add a bin to this histogram with the specified bin contents
495/// and error. The bin width is used to set the relative scale of
496/// bins with different widths.
497
498void RooHist::addBinWithXYError(Axis_t binCenter, double n, double exlow, double exhigh, double eylow, double eyhigh,
499 double scaleFactor)
500{
501 _entries+= n;
502
503 addPoint(binCenter, n, scaleFactor,exlow,exhigh, eylow, eyhigh);
504}
505
506
507
508
509
510////////////////////////////////////////////////////////////////////////////////
511/// Add a bin to this histogram with the value (n1-n2)/(n1+n2)
512/// using an error bar calculated with Binomial statistics.
513
514void RooHist::addAsymmetryBin(Axis_t binCenter, Int_t n1, Int_t n2, double binWidth, double xErrorFrac, double scaleFactor)
515{
516 // calculate Binomial errors for this bin
517 double ym;
518 double yp;
519 double dx(0.5 * binWidth);
520 if(!RooHistError::instance().getBinomialIntervalAsym(n1,n2,ym,yp,_nSigma)) {
521 coutE(Plotting) << "RooHist::addAsymmetryBin: unable to calculate binomial error for bin with " << n1 << "," << n2 << " events" << std::endl;
522 return;
523 }
524
525 const Int_t denominator = n1 + n2;
526 double a = 0 == denominator ? 0. : (double)(n1 - n2) / (denominator);
527 addPoint(binCenter, a, scaleFactor,dx*xErrorFrac,dx*xErrorFrac, a-ym, yp-a);
528}
529
530
531
532////////////////////////////////////////////////////////////////////////////////
533/// Add a bin to this histogram with the value (n1-n2)/(n1+n2)
534/// using an error bar calculated with Binomial statistics.
535
536void RooHist::addAsymmetryBinWithError(Axis_t binCenter, double n1, double n2, double en1, double en2, double binWidth, double xErrorFrac, double scaleFactor)
537{
538 // calculate Binomial errors for this bin
539 double ym;
540 double yp;
541 double dx(0.5 * binWidth);
542 double a= (double)(n1-n2)/(n1+n2);
543
544 double error = 2*sqrt( pow(en1,2)*pow(n2,2) + pow(en2,2)*pow(n1,2) ) / pow(n1+n2,2) ;
545 ym=a-error ;
546 yp=a+error ;
547
548 addPoint(binCenter,a, scaleFactor, dx*xErrorFrac,dx*xErrorFrac, a-ym, yp-a);
549}
550
551
552
553////////////////////////////////////////////////////////////////////////////////
554/// Add a bin to this histogram with the value n1/(n1+n2)
555/// using an error bar calculated with Binomial statistics.
556
557void RooHist::addEfficiencyBin(Axis_t binCenter, Int_t n1, Int_t n2, double binWidth, double xErrorFrac, double scaleFactor)
558{
559 double a= (double)(n1)/(n1+n2);
560
561 // calculate Binomial errors for this bin
562 double ym;
563 double yp;
564 double dx(0.5 * binWidth);
565 if(!RooHistError::instance().getBinomialIntervalEff(n1,n2,ym,yp,_nSigma)) {
566 coutE(Plotting) << "RooHist::addEfficiencyBin: unable to calculate binomial error for bin with " << n1 << "," << n2 << " events" << std::endl;
567 return;
568 }
569
570 addPoint(binCenter,a, scaleFactor,dx*xErrorFrac,dx*xErrorFrac, a-ym, yp-a);
571}
572
573
574
575////////////////////////////////////////////////////////////////////////////////
576/// Add a bin to this histogram with the value n1/(n1+n2)
577/// using an error bar calculated with Binomial statistics.
578
579void RooHist::addEfficiencyBinWithError(Axis_t binCenter, double n1, double n2, double en1, double en2, double binWidth, double xErrorFrac, double scaleFactor)
580{
581 double a= (double)(n1)/(n1+n2);
582
583 double error = sqrt( pow(en1,2)*pow(n2,2) + pow(en2,2)*pow(n1,2) ) / pow(n1+n2,2) ;
584
585 // calculate Binomial errors for this bin
586 double ym;
587 double yp;
588 double dx(0.5 * binWidth);
589 ym=a-error ;
590 yp=a+error ;
591
592
593 addPoint(binCenter,a, scaleFactor,dx*xErrorFrac,dx*xErrorFrac, a-ym, yp-a);
594}
595
596
597////////////////////////////////////////////////////////////////////////////////
598/// Return true if binning of this RooHist is identical to that of 'other'
599
600bool RooHist::hasIdenticalBinning(const RooHist& other) const
601{
602 // First check if number of bins is the same
603 if (GetN() != other.GetN()) {
604 return false ;
605 }
606
607 // Next require that all bin centers are the same
608 Int_t i ;
609 for (i=0 ; i<GetN() ; i++) {
610 double x1;
611 double x2;
612 double y1;
613 double y2;
614
615 GetPoint(i,x1,y1) ;
616 other.GetPoint(i,x2,y2) ;
617
618 if (std::abs(x1-x2) > 1e-10 * _nominalBinWidth) {
619 return false ;
620 }
621
622 }
623
624 return true ;
625}
626
627
628
629////////////////////////////////////////////////////////////////////////////////
630/// Return true if contents of this RooHist is identical within given
631/// relative tolerance to that of 'other'
632
633bool RooHist::isIdentical(const RooHist& other, double tol, bool verbose) const
634{
635 // Make temporary TH1s output of RooHists to perform Kolmogorov test
636 TH1::AddDirectory(false) ;
637 TH1F h_self("h_self","h_self",GetN(),0,1) ;
638 TH1F h_other("h_other","h_other",GetN(),0,1) ;
639 TH1::AddDirectory(true) ;
640
641 for (Int_t i=0 ; i<GetN() ; i++) {
642 h_self.SetBinContent(i+1,GetY()[i]) ;
643 h_other.SetBinContent(i+1,other.GetY()[i]) ;
644 }
645
646 double M = h_self.KolmogorovTest(&h_other,"M") ;
647 if (M>tol) {
648 double kprob = h_self.KolmogorovTest(&h_other) ;
649 if(verbose) std::cout << "RooHist::isIdentical() tolerance exceeded M=" << M << " (tol=" << tol << "), corresponding prob = " << kprob << std::endl ;
650 return false ;
651 }
652
653 return true ;
654}
655
656
657
658////////////////////////////////////////////////////////////////////////////////
659/// Print info about this histogram to the specified output stream.
660///
661/// Standard: number of entries
662/// Shape: error CL and maximum value
663/// Verbose: print our bin contents and errors
664
665void RooHist::printMultiline(std::ostream& os, Int_t contents, bool verbose, TString indent) const
666{
667 RooPlotable::printMultiline(os,contents,verbose,indent);
668 os << indent << "--- RooHist ---" << std::endl;
669 Int_t n= GetN();
670 os << indent << " Contains " << n << " bins" << std::endl;
671 if(verbose) {
672 os << indent << " Errors calculated at" << _nSigma << "-sigma CL" << std::endl;
673 os << indent << " Bin Contents:" << std::endl;
674 for(Int_t i= 0; i < n; i++) {
675 os << indent << std::setw(3) << i << ") x= " << fX[i];
676 if(fEXhigh[i] > 0 || fEXlow[i] > 0) {
677 os << " +" << fEXhigh[i] << " -" << fEXlow[i];
678 }
679 os << " , y = " << fY[i] << " +" << fEYhigh[i] << " -" << fEYlow[i] << std::endl;
680 }
681 }
682}
683
684
685
686////////////////////////////////////////////////////////////////////////////////
687/// Print name of RooHist
688
689void RooHist::printName(std::ostream& os) const
690{
691 os << GetName() ;
692}
693
694
695
696////////////////////////////////////////////////////////////////////////////////
697/// Print title of RooHist
698
699void RooHist::printTitle(std::ostream& os) const
700{
701 os << GetTitle() ;
702}
703
704
705
706////////////////////////////////////////////////////////////////////////////////
707/// Print class name of RooHist
708
709void RooHist::printClassName(std::ostream& os) const
710{
711 os << ClassName() ;
712}
713
714
715std::unique_ptr<RooHist> RooHist::createEmptyResidHist(const RooCurve& curve, bool normalize) const
716{
717 // Copy all non-content properties from hist1
718 auto hist = std::make_unique<RooHist>(_nominalBinWidth) ;
719 const std::string name = GetName() + std::string("_") + curve.GetName();
720 const std::string title = GetTitle() + std::string(" and ") + curve.GetTitle();
721 hist->SetName(((normalize ? "pull_" : "resid_") + name).c_str()) ;
722 hist->SetTitle(((normalize ? "Pull of " : "Residual of ") + title).c_str()) ;
723
724 return hist;
725}
726
727
728void RooHist::fillResidHist(RooHist & residHist, const RooCurve& curve,bool normalize, bool useAverage) const
729{
730 // Determine range of curve
731 double xstart;
732 double xstop;
733 double y;
734 curve.GetPoint(0,xstart,y) ;
735 curve.GetPoint(curve.GetN()-1,xstop,y) ;
736
737 // Add histograms, calculate Poisson confidence interval on sum value
738 for(Int_t i=0 ; i<GetN() ; i++) {
739 double x;
740 double point;
741 GetPoint(i,x,point) ;
742
743 // Only calculate pull for bins inside curve range
744 if (x<xstart || x>xstop) continue ;
745
746 double yy ;
747 if (useAverage) {
748 double exl = GetErrorXlow(i);
749 double exh = GetErrorXhigh(i) ;
750 if (exl<=0 ) exl = GetErrorX(i);
751 if (exh<=0 ) exh = GetErrorX(i);
752 if (exl<=0 ) exl = 0.5*getNominalBinWidth();
753 if (exh<=0 ) exh = 0.5*getNominalBinWidth();
754 yy = point - curve.average(x-exl,x+exh) ;
755 } else {
756 yy = point - curve.interpolate(x) ;
757 }
758
759 double dyl = GetErrorYlow(i) ;
760 double dyh = GetErrorYhigh(i) ;
761 if (normalize) {
762 double norm = (yy>0?dyl:dyh);
763 if (norm==0.) {
764 coutW(Plotting) << "RooHist::makeResisHist(" << GetName() << ") WARNING: point " << i << " has zero error, setting residual to zero" << std::endl;
765 yy=0 ;
766 dyh=0 ;
767 dyl=0 ;
768 } else {
769 yy /= norm;
770 dyh /= norm;
771 dyl /= norm;
772 }
773 }
774 residHist.addBinWithError(x,yy,dyl,dyh);
775 }
776}
777
778
779////////////////////////////////////////////////////////////////////////////////
780/// Create and return RooHist containing residuals w.r.t to given curve.
781/// If normalize is true, the residuals are normalized by the histogram
782/// errors creating a RooHist with pull values
783
784RooHist* RooHist::makeResidHist(const RooCurve& curve, bool normalize, bool useAverage) const
785{
786 RooHist* hist = createEmptyResidHist(curve, normalize).release();
787 fillResidHist(*hist, curve, normalize, useAverage);
788 return hist ;
789}
#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:382
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:2489
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:24
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:655
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:580
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.
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:784
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:709
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:557
void fillResidHist(RooHist &residHist, const RooCurve &curve, bool normalize=false, bool useAverage=false) const
Definition RooHist.cxx:728
void initialize()
Perform common initialization for all constructors.
Definition RooHist.cxx:327
void addBinWithError(Axis_t binCenter, double n, double elow, double ehigh, double binWidth=0, double xErrorFrac=1.0, bool correctForBinWidth=true, double scaleFactor=1.0)
Add a bin to this histogram with the specified bin contents and error.
Definition RooHist.cxx:477
void addEfficiencyBinWithError(Axis_t binCenter, double n1, double n2, double en1, double en2, double binWidth=0, double xErrorFrac=1.0, double scaleFactor=1.0)
Add a bin to this histogram with the value n1/(n1+n2) using an error bar calculated with Binomial sta...
Definition RooHist.cxx:579
RooHist()
Definition RooHist.h:31
void printName(std::ostream &os) const override
Print name of RooHist.
Definition RooHist.cxx:689
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:665
std::unique_ptr< RooHist > createEmptyResidHist(const RooCurve &curve, bool normalize=false) const
Definition RooHist.cxx:715
void printTitle(std::ostream &os) const override
Print title of RooHist.
Definition RooHist.cxx:699
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:536
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:384
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:633
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:514
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:426
void addPoint(Axis_t binCenter, double y, double yscale, double exlow, double exhigh, double eylow, double eyhigh)
Definition RooHist.cxx:398
bool hasIdenticalBinning(const RooHist &other) const
Return true if binning of this RooHist is identical to that of 'other'.
Definition RooHist.cxx:600
double getFitRangeNEvt() const override
Return the number of events of the dataset associated with this RooHist.
Definition RooHist.cxx:338
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:498
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
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:2342
Double_t * GetY() const
Definition TGraph.h:140
Int_t GetN() const
Definition TGraph.h:132
Double_t * fY
[fNpoints] array of Y points
Definition TGraph.h:48
void SetName(const char *name="") override
Set graph name.
Definition TGraph.cxx:2381
Double_t * fX
[fNpoints] array of X points
Definition TGraph.h:47
void SetTitle(const char *title="") override
Change (i.e.
Definition TGraph.cxx:2397
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:1533
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:622
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition TH1.cxx:9162
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9084
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition TH1.cxx:1294
TAxis * GetXaxis()
Definition TH1.h:324
virtual Int_t GetNbinsX() const
Definition TH1.h:297
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:9243
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5082
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition TH1.cxx:9184
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:8199
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:225
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