ROOT   Reference Guide
Searching...
No Matches
RooCurve.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 *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
15 *****************************************************************************/
16
17/**
18\file RooCurve.cxx
19\class RooCurve
20\ingroup Roofitcore
21
22One-dimensional graphical representation of a real-valued function.
23A curve is approximated by straight line segments with end points chosen to give
24a "good" approximation to the true curve. The goodness of the approximation is
25controlled by a precision and a resolution parameter.
26
27A RooCurve derives from TGraph, so it can either be drawn as a line (default) or
28as points:
29
30RooPlot *p = y.plotOn(x.frame());
31p->getAttMarker("curve_y")->SetMarkerStyle(20);
32p->setDrawOptions("curve_y","PL");
33p->Draw();
34
35
36To retrieve a RooCurve from a RooPlot, use RooPlot::getCurve().
37**/
38
39#include "RooCurve.h"
40#include "RooHist.h"
41#include "RooAbsReal.h"
42#include "RooArgSet.h"
43#include "RooRealVar.h"
44#include "RooRealBinding.h"
45#include "RooMsgService.h"
46#include "RooProduct.h"
47#include "RooConstVar.h"
48
49#include "Riostream.h"
50#include "TMath.h"
51#include "TAxis.h"
52#include "TMatrixD.h"
53#include "TVectorD.h"
54#include "Math/Util.h"
55#include <iomanip>
56#include <deque>
57#include <algorithm>
58
59using std::endl, std::ostream, std::list, std::vector, std::cout, std::min;
60
62
63namespace {
64
65// Helpers to manage points
66struct Point {
67 double x;
68 double y;
69};
70
71inline Point getPoint(TGraph const &gr, int i)
72{
73 Point p;
74 gr.GetPoint(i, p.x, p.y);
75 return p;
76}
77
78} // namespace
79
80
81////////////////////////////////////////////////////////////////////////////////
82/// Default constructor.
83
85{
86 initialize();
87}
88
89
90////////////////////////////////////////////////////////////////////////////////
91/// Create a 1-dim curve of the value of the specified real-valued expression
92/// as a function of x. Use the optional precision parameter to control
93/// how precisely the smooth curve is rasterized. Use the optional argument set
94/// to specify how the expression should be normalized. Use the optional scale
95/// factor to rescale the expression after normalization.
96/// If shiftToZero is set, the entire curve is shifted down to make the lowest
97/// point of the curve go through zero.
98RooCurve::RooCurve(const RooAbsReal &f, RooAbsRealLValue &x, double xlo, double xhi, Int_t xbins, double scaleFactor,
99 const RooArgSet *normVars, double prec, double resolution, bool shiftToZero, WingMode wmode,
100 Int_t nEvalError, Int_t doEEVal, double eeVal, bool showProg)
101 : _showProgress(showProg)
102{
103 // grab the function's name and title
104 TString name(f.GetName());
105 SetName(name.Data());
106 TString title(f.GetTitle());
107 SetTitle(title.Data());
108 // append " ( [<funit> ][/ <xunit> ])" to our y-axis label if necessary
109 if(0 != strlen(f.getUnit()) || 0 != strlen(x.getUnit())) {
110 title.Append(" ( ");
111 if(0 != strlen(f.getUnit())) {
112 title.Append(f.getUnit());
113 title.Append(" ");
114 }
115 if(0 != strlen(x.getUnit())) {
116 title.Append("/ ");
117 title.Append(x.getUnit());
118 title.Append(" ");
119 }
120 title.Append(")");
121 }
122 setYAxisLabel(title.Data());
123
124 RooProduct scaledFunc{"scaled_func", "scaled_func", {f, RooFit::RooConst(scaleFactor)}};
125 std::unique_ptr<RooAbsFunc> funcPtr{scaledFunc.bindVars(x, normVars, true)};
126
127 // calculate the points to add to our curve
128 if(xbins > 0){
129 // regular mode - use the sampling hint to decide where to evaluate the pdf
130 std::unique_ptr<std::list<double>> hint{f.plotSamplingHint(x,xlo,xhi)};
132 if (_showProgress) {
133 ccoutP(Plotting) << endl ;
134 }
135 } else {
136 // if number of bins is set to <= 0, skip any interpolation and just evaluate the pdf at the bin centers
137 // this is useful when plotting a pdf like a histogram
138 int nBinsX = x.numBins();
139 for(int i=0; i<nBinsX; ++i){
140 double xval = x.getBinning().binCenter(i);
142 }
143 }
144 initialize();
145
146 if (shiftToZero) shiftCurveToZero() ;
147
149 for (int i=0 ; i<GetN() ; i++) {
151 }
152 this->Sort();
153}
154
155
156
157////////////////////////////////////////////////////////////////////////////////
158/// Create a 1-dim curve of the value of the specified real-valued
159/// expression as a function of x. Use the optional precision
160/// parameter to control how precisely the smooth curve is
161/// rasterized. If shiftToZero is set, the entire curve is shifted
162/// down to make the lowest point in of the curve go through zero.
163
164RooCurve::RooCurve(const char *name, const char *title, const RooAbsFunc &func,
165 double xlo, double xhi, UInt_t minPoints, double prec, double resolution,
166 bool shiftToZero, WingMode wmode, Int_t nEvalError, Int_t doEEVal, double eeVal)
167{
168 SetName(name);
169 SetTitle(title);
171 initialize();
172 if (shiftToZero) shiftCurveToZero() ;
173
175 for (int i=0 ; i<GetN() ; i++) {
177 }
178 this->Sort();
179}
180
181
182
183////////////////////////////////////////////////////////////////////////////////
184/// Constructor of a curve as sum of two other curves.
185/// \f[
186/// C_\mathrm{sum} = \mathrm{scale1}*c1 + \mathrm{scale2}*c2
187/// \f]
188///
189/// \param[in] name Name of the curve (to retrieve it from a plot)
190/// \param[in] title Title (for plotting).
191/// \param[in] c1 First curve.
192/// \param[in] c2 Second curve.
193/// \param[in] scale1 Scale y values for c1 by this factor.
194/// \param[in] scale2 Scale y values for c2 by this factor.
195
196RooCurve::RooCurve(const char* name, const char* title, const RooCurve& c1, const RooCurve& c2, double scale1, double scale2)
197{
198 initialize() ;
199 SetName(name) ;
200 SetTitle(title) ;
201
202 // Make deque of points in X
203 std::deque<double> pointList ;
204
205 // Add X points of C1
206 Int_t n1 = c1.GetN();
207 for (int i1=0 ; i1<n1 ; i1++) {
208 pointList.push_back(c1.GetPointX(i1));
209 }
210
211 // Add X points of C2
212 Int_t n2 = c2.GetN();
213 for (int i2=0 ; i2<n2 ; i2++) {
214 pointList.push_back(c2.GetPointX(i2));
215 }
216
217 // Sort X points
218 std::sort(pointList.begin(),pointList.end()) ;
219
220 // Loop over X points
221 double last(-RooNumber::infinity()) ;
222 for (auto point : pointList) {
223
224 if ((point-last)>1e-10) {
225 // Add OR of points to new curve, skipping duplicate points within tolerance
227 }
228 last = point ;
229 }
230
231 this->Sort();
232}
233
234
235RooCurve::~RooCurve() = default;
236
237
238////////////////////////////////////////////////////////////////////////////////
239/// Perform initialization that is common to all curves
240
242{
243 // set default line width in pixels
244 SetLineWidth(3);
245 // set default line color
247}
248
249
250
251////////////////////////////////////////////////////////////////////////////////
252/// Find lowest point in curve and move all points in curve so that
253/// lowest point will go exactly through zero
254
256{
257 double minVal = std::numeric_limits<double>::infinity();
258 double maxVal = -std::numeric_limits<double>::infinity();
259
260 // First iteration, find current lowest point
261 for (int i = 1; i < GetN() - 1; i++) {
262 double y = GetPointY(i);
263 minVal = std::min(y, minVal);
264 maxVal = std::max(y, maxVal);
265 }
266
267 // Second iteration, lower all points by minVal
268 for (int i = 1; i < GetN() - 1; i++) {
269 Point point = getPoint(*this, i);
270 SetPoint(i, point.x, point.y - minVal);
271 }
272
273 setYAxisLimits(0, maxVal - minVal);
274}
275
276
277
278////////////////////////////////////////////////////////////////////////////////
279/// Add points calculated with the specified function, over the range (xlo,xhi).
280/// Add at least minPoints equally spaced points, and add sufficient points so that
281/// the maximum deviation from the final straight-line segments is prec*(ymax-ymin),
282/// down to a minimum horizontal spacing of resolution*(xhi-xlo).
283
284void RooCurve::addPoints(const RooAbsFunc &func, double xlo, double xhi,
285 Int_t minPoints, double prec, double resolution, WingMode wmode,
286 Int_t numee, bool doEEVal, double eeVal, list<double>* samplingHint)
287{
288 // check the inputs
289 if(!func.isValid()) {
290 coutE(InputArguments) << fName << "::addPoints: input function is not valid" << endl;
291 return;
292 }
293 if(minPoints <= 0 || xhi <= xlo) {
295 return;
296 }
297
298 // Perform a coarse scan of the function to estimate its y range.
299 // Save the results so we do not have to re-evaluate at the scan points.
300
301 // Adjust minimum number of points to external sampling hint if used
302 if (samplingHint) {
303 minPoints = samplingHint->size() ;
304 }
305
306 double dx= (xhi-xlo)/(minPoints-1.);
307 const double epsilon = (xhi - xlo) * relativeXEpsilon();
308 std::vector<double> yval(minPoints);
309
310 // Get list of initial x values. If function provides sampling hint use that,
311 // otherwise use default binning of frame
312 std::vector<double> xval;
313 if (!samplingHint) {
314 for(int step= 0; step < minPoints; step++) {
315 xval.push_back(xlo + step*dx) ;
316 }
317 } else {
318 std::copy(samplingHint->begin(), samplingHint->end(), std::back_inserter(xval));
319 }
320
321 for (unsigned int step=0; step < xval.size(); ++step) {
322 double xx = xval[step];
323 if (step == static_cast<unsigned int>(minPoints-1))
324 xx -= 1e-9 * dx;
325
326 yval[step]= func(&xx);
327 if (_showProgress) {
328 ccoutP(Plotting) << "." ;
329 cout.flush() ;
330 }
331
333 if (numee>=0) {
334 coutW(Plotting) << "At observable [x]=" << xx << " " ;
335 RooAbsReal::printEvalErrors(ccoutW(Plotting),numee) ;
336 }
337 if (doEEVal) {
338 yval[step]=eeVal ;
339 }
340 }
342 }
343
344 const double ymax = *std::max_element(yval.begin(), yval.end());
345 const double ymin = *std::min_element(yval.begin(), yval.end());
346 double yrangeEst=(ymax-ymin) ;
347
348 // store points of the coarse scan and calculate any refinements necessary
349 double minDx= resolution*(xhi-xlo);
350 double x1;
351 double x2 = xlo;
352
353 if (wmode==Extended) {
354 // Add two points to make curve jump from 0 to yval at the left end of the plotting range.
355 // This ensures that filled polygons are drawn properly. The first point needs to be to the
356 // left of the second, so it's shifted by 1/1000 more than the second.
359 } else if (wmode==Straight) {
361 }
362
364
365 auto iter2 = xval.begin() ;
366 x1 = *iter2 ;
367 int step=1 ;
368 while(true) {
369 x1= x2;
370 ++iter2 ;
371 if (iter2==xval.end()) {
372 break ;
373 }
374 x2= *iter2 ;
375 if (prec<0) {
376 // If precision is <0, no attempt at recursive interpolation is made
378 } else {
380 }
381 step++ ;
382 }
384
385 if (wmode==Extended) {
386 // Add two points to close polygon. The order matters. Since they are sorted in x later, the second
387 // point is shifted by 1/1000 more than the second-to-last point.
390 } else if (wmode==Straight) {
392 }
393}
394
395
396////////////////////////////////////////////////////////////////////////////////
397/// Fill the range (x1,x2) with points calculated using func(&x). No point will
398/// be added at x1, and a point will always be added at x2. The density of points
399/// will be calculated so that the maximum deviation from a straight line
400/// approximation is prec*(ymax-ymin) down to the specified minimum horizontal spacing.
401
402void RooCurve::addRange(const RooAbsFunc& func, double x1, double x2,
403 double y1, double y2, double minDy, double minDx,
404 int numee, bool doEEVal, double eeVal, double epsilon)
405{
406 // Explicitly skip empty ranges to eliminate point duplication
407 if (std::abs(x2-x1) <= epsilon) {
408 return ;
409 }
410
411 // calculate our value at the midpoint of this range
412 double xmid= 0.5*(x1+x2);
413 double ymid= func(&xmid);
414 if (_showProgress) {
415 ccoutP(Plotting) << "." ;
416 cout.flush() ;
417 }
418
420 if (numee>=0) {
421 coutW(Plotting) << "At observable [x]=" << xmid << " " ;
422 RooAbsReal::printEvalErrors(ccoutW(Plotting),numee) ;
423 }
424 if (doEEVal) {
425 ymid=eeVal ;
426 }
427 }
429
430 // test if the midpoint is sufficiently close to a straight line across this interval
431 double dy= ymid - 0.5*(y1+y2);
432 if((xmid - x1 >= minDx) && std::abs(dy)>0 && std::abs(dy) >= minDy) {
433 // fill in each subrange
436 }
437 else {
440 }
441}
442
443
444////////////////////////////////////////////////////////////////////////////////
445/// Add a point with the specified coordinates. Update our y-axis limits.
446
448{
449// cout << "RooCurve("<< GetName() << ") adding point at (" << x << "," << y << ")" << endl ;
450 Int_t next= GetN();
451 SetPoint(next, x, y);
453}
454
455
456////////////////////////////////////////////////////////////////////////////////
457/// Return the number of events associated with the plotable object,
458/// it is always 1 for curves
459
461 return 1;
462}
463
464
465////////////////////////////////////////////////////////////////////////////////
466/// Return the number of events associated with the plotable object,
467/// in the given range. It is always 1 for curves
468
469double RooCurve::getFitRangeNEvt(double, double) const
470{
471 return 1 ;
472}
473
474
475////////////////////////////////////////////////////////////////////////////////
476/// Get the bin width associated with this plotable object.
477/// It is alwats zero for curves
478
480 return 0 ;
481}
482
483
484
485////////////////////////////////////////////////////////////////////////////////
486
487void RooCurve::printName(ostream& os) const
488//
489{
490 // Print the name of this curve
491 os << GetName() ;
492}
493
494
495////////////////////////////////////////////////////////////////////////////////
496/// Print the title of this curve
497
498void RooCurve::printTitle(ostream& os) const
499{
500 os << GetTitle() ;
501}
502
503
504////////////////////////////////////////////////////////////////////////////////
505/// Print the class name of this curve
506
507void RooCurve::printClassName(ostream& os) const
508{
509 os << ClassName() ;
510}
511
512
513
514////////////////////////////////////////////////////////////////////////////////
515/// Print the details of this curve
516
517void RooCurve::printMultiline(ostream& os, Int_t /*contents*/, bool /*verbose*/, TString indent) const
518{
519 os << indent << "--- RooCurve ---" << endl ;
520 Int_t n= GetN();
521 os << indent << " Contains " << n << " points" << endl;
522 os << indent << " Graph points:" << endl;
523 for(Int_t i= 0; i < n; i++) {
524 os << indent << std::setw(3) << i << ") x = " << fX[i] << " , y = " << fY[i] << endl;
525 }
526}
527
528
529
530////////////////////////////////////////////////////////////////////////////////
531/// Calculate the chi^2/NDOF of this curve with respect to the histogram
532/// 'hist' accounting nFitParam floating parameters in case the curve
533/// was the result of a fit
534
535double RooCurve::chiSquare(const RooHist& hist, Int_t nFitParam) const
536{
537 Int_t np = hist.GetN();
538
539 // Find starting and ending bin of histogram based on range of RooCurve
540 double xstart = GetPointX(0);
541 double xstop = GetPointX(GetN()-1);
542
543 Int_t nbin(0) ;
544
546 for (int i=0 ; i<np ; i++) {
547
548 // Retrieve histogram contents
549 Point point = getPoint(hist, i);
550
551 // Check if point is in range of curve
552 if (point.x<xstart || point.x>xstop) continue ;
553
554 double eyl = hist.GetEYlow()[i] ;
555 double eyh = hist.GetEYhigh()[i] ;
556 double exl = hist.GetEXlow()[i] ;
557 double exh = hist.GetEXhigh()[i] ;
558
559 // Integrate function over this bin
560 double avg = average(point.x-exl,point.x+exh) ;
561
562 // Add pull^2 to chisq
563 if (point.y!=0) {
564 double pull = (point.y>avg) ? ((point.y-avg)/eyl) : ((point.y-avg)/eyh) ;
565 chisq += pull*pull ;
566 nbin++ ;
567 }
568 }
569
570 // Return chisq/nDOF
571 return chisq.Sum() / (nbin-nFitParam) ;
572}
573
574
575
576////////////////////////////////////////////////////////////////////////////////
577/// Return average curve value in [xFirst,xLast] by integrating curve between points
578/// and dividing by xLast-xFirst
579
580double RooCurve::average(double xFirst, double xLast) const
581{
582 if (xFirst>=xLast) {
583 coutE(InputArguments) << "RooCurve::average(" << GetName()
584 << ") invalid range (" << xFirst << "," << xLast << ")" << endl ;
585 return 0 ;
586 }
587
588 // Find Y values and begin and end points
589 double yFirst = interpolate(xFirst,1e-10) ;
590 double yLast = interpolate(xLast,1e-10) ;
591
592 // Find first and last mid points
593 Int_t ifirst = findPoint(xFirst, std::numeric_limits<double>::infinity());
594 Int_t ilast = findPoint(xLast, std::numeric_limits<double>::infinity());
595
596 // Make sure the midpoints are actually in the interval
597 while (GetPointX(ifirst) < xFirst) {
598 ++ifirst;
599 }
600 while (GetPointX(ilast) > xLast) {
601 --ilast;
602 }
603
604 // Handle trivial scenario -- no midway points, point only at or outside given range
605 if (ilast < ifirst) {
606 return 0.5*(yFirst+yLast) ;
607 }
608
609 Point firstPt = getPoint(*this, ifirst);
610 Point lastPt = getPoint(*this, ilast);
611
612 // Trapezoid integration from lower edge to first midpoint
613 double sum = 0.5 * (firstPt.x-xFirst)*(yFirst+firstPt.y);
614
615 // Trapezoid integration between midpoints
616 for (int i=ifirst ; i<ilast ; i++) {
617 Point p1 = getPoint(*this, i) ;
618 Point p2 = getPoint(*this, i+1) ;
619 sum += 0.5 * (p2.x-p1.x)*(p1.y+p2.y);
620 }
621
622 // Trapezoid integration from last midpoint to upper edge
623 sum += 0.5 * (xLast-lastPt.x)*(lastPt.y+yLast);
624 return sum/(xLast-xFirst) ;
625}
626
627
628
629////////////////////////////////////////////////////////////////////////////////
630/// Find the nearest point to xvalue. Return -1 if distance
631/// exceeds tolerance
632
633Int_t RooCurve::findPoint(double xvalue, double tolerance) const
634{
635 double delta(std::numeric_limits<double>::max());
636 Int_t n = GetN();
637 Int_t ibest(-1) ;
638 for (int i=0 ; i<n ; i++) {
639 double x = GetPointX(i);
640 if (std::abs(xvalue-x)<delta) {
641 delta = std::abs(xvalue-x) ;
642 ibest = i ;
643 }
644 }
645
646 return (delta<tolerance)?ibest:-1 ;
647}
648
649
650////////////////////////////////////////////////////////////////////////////////
651/// Return linearly interpolated value of curve at xvalue. If distance
652/// to nearest point is less than tolerance, return nearest point value
654
655double RooCurve::interpolate(double xvalue, double tolerance) const
656{
657 // Find best point
658 int n = GetN() ;
659 int ibest = findPoint(xvalue,1e10) ;
660
661 // Get position of best point
662 Point pbest = getPoint(*this, ibest);
663
664 // Handle trivial case of being dead on
665 if (std::abs(pbest.x-xvalue)<tolerance) {
666 return pbest.y;
667 }
668
669 // Get nearest point on other side w.r.t. xvalue
670 double retVal(0);
671 if (pbest.x<xvalue) {
672 if (ibest==n-1) {
673 // Value beyond end requested -- return value of last point
674 return pbest.y ;
675 }
676 Point pother = getPoint(*this, ibest+1);
677 if (pother.x==pbest.x) return pbest.y ;
678 retVal = pbest.y + (pother.y-pbest.y)*(xvalue-pbest.x)/(pother.x-pbest.x) ;
679
680 } else {
681 if (ibest==0) {
682 // Value before 1st point requested -- return value of 1st point
683 return pbest.y ;
684 }
685 Point pother = getPoint(*this, ibest-1);
686 if (pother.x==pbest.x) return pbest.y ;
687 retVal = pother.y + (pbest.y-pother.y)*(xvalue-pother.x)/(pbest.x-pother.x) ;
688 }
689
690 return retVal ;
691}
692
693
694
695
696////////////////////////////////////////////////////////////////////////////////
697/// Construct filled RooCurve represented error band that captures alpha% of the variations
698/// of the curves passed through argument variations, where the percentage alpha corresponds to
699/// the central interval fraction of a significance Z
700
701RooCurve* RooCurve::makeErrorBand(const vector<RooCurve*>& variations, double Z) const
702{
703 RooCurve* band = new RooCurve ;
704 band->SetName((std::string(GetName()) + "_errorband").c_str());
705 band->SetLineWidth(1) ;
706 band->SetFillColor(kCyan) ;
707 band->SetLineColor(kCyan) ;
708
709 vector<double> bandLo(GetN()) ;
710 vector<double> bandHi(GetN()) ;
711 for (int i=0 ; i<GetN() ; i++) {
712 calcBandInterval(variations,i,Z,bandLo[i],bandHi[i],false) ;
713 }
714
715 for (int i=0 ; i<GetN() ; i++) {
717 }
718 for (int i=GetN()-1 ; i>=0 ; i--) {
720 }
721 // if the axis of the old graph is alphanumeric, copy the labels to the new one as well
722 if(this->GetXaxis() && this->GetXaxis()->IsAlphanumeric()){
723 band->GetXaxis()->Set(this->GetXaxis()->GetNbins(),this->GetXaxis()->GetXmin(),this->GetXaxis()->GetXmax());
724 for(int i=0; i<this->GetXaxis()->GetNbins(); ++i){
725 band->GetXaxis()->SetBinLabel(i+1,this->GetXaxis()->GetBinLabel(i+1));
726 }
727 }
728
729 return band ;
730}
731
732
733
734
735////////////////////////////////////////////////////////////////////////////////
736/// Construct filled RooCurve represented error band represent the error added in quadrature defined by the curves arguments
737/// plusVar and minusVar corresponding to one-sigma variations of each parameter. The resulting error band, combined used the correlation matrix C
738/// is multiplied with the significance parameter Z to construct the equivalent of a Z sigma error band (in Gaussian approximation)
739
740RooCurve* RooCurve::makeErrorBand(const vector<RooCurve*>& plusVar, const vector<RooCurve*>& minusVar, const TMatrixD& C, double Z) const
741{
742
743 RooCurve* band = new RooCurve ;
744 band->SetName((std::string(GetName()) + "_errorband").c_str());
745 band->SetLineWidth(1) ;
746 band->SetFillColor(kCyan) ;
747 band->SetLineColor(kCyan) ;
748
749 vector<double> bandLo(GetN()) ;
750 vector<double> bandHi(GetN()) ;
751 for (int i=0 ; i<GetN() ; i++) {
752 calcBandInterval(plusVar,minusVar,i,C,Z,bandLo[i],bandHi[i]) ;
753 }
754
755 for (int i=0 ; i<GetN() ; i++) {
757 }
758 for (int i=GetN()-1 ; i>=0 ; i--) {
760 }
761
762 // if the axis of the old graph is alphanumeric, copy the labels to the new one as well
763 if(this->GetXaxis() && this->GetXaxis()->IsAlphanumeric()){
764 band->GetXaxis()->Set(this->GetXaxis()->GetNbins(),this->GetXaxis()->GetXmin(),this->GetXaxis()->GetXmax());
765 for(int i=0; i<this->GetXaxis()->GetNbins(); ++i){
766 band->GetXaxis()->SetBinLabel(i+1,this->GetXaxis()->GetBinLabel(i+1));
767 }
768 }
769
770 return band ;
771}
772
773
774
775
776
777////////////////////////////////////////////////////////////////////////////////
778/// Retrieve variation points from curves
779
780void RooCurve::calcBandInterval(const vector<RooCurve*>& plusVar, const vector<RooCurve*>& minusVar,Int_t i, const TMatrixD& C, double /*Z*/, double& lo, double& hi) const
781{
782 vector<double> y_plus(plusVar.size());
783 vector<double> y_minus(minusVar.size());
784 Int_t j(0) ;
785 for (vector<RooCurve*>::const_iterator iter=plusVar.begin() ; iter!=plusVar.end() ; ++iter) {
786 y_plus[j++] = (*iter)->interpolate(GetX()[i]) ;
787 }
788 j=0 ;
789 for (vector<RooCurve*>::const_iterator iter=minusVar.begin() ; iter!=minusVar.end() ; ++iter) {
790 y_minus[j++] = (*iter)->interpolate(GetX()[i]) ;
791 }
792 double y_cen = GetY()[i] ;
793 Int_t n = j ;
794
795 // Make vector of variations
796 TVectorD F(plusVar.size()) ;
797 for (j=0 ; j<n ; j++) {
798 F[j] = (y_plus[j]-y_minus[j])/2 ;
799 }
800
801 // Calculate error in linear approximation from variations and correlation coefficient
802 double sum = F*(C*F) ;
803
804 lo= y_cen + sqrt(sum) ;
805 hi= y_cen - sqrt(sum) ;
806}
807
808
809
810////////////////////////////////////////////////////////////////////////////////
811
812void RooCurve::calcBandInterval(const vector<RooCurve*>& variations,Int_t i,double Z, double& lo, double& hi, bool approxGauss) const
813{
814 vector<double> y(variations.size()) ;
815 Int_t j(0) ;
816 for (vector<RooCurve*>::const_iterator iter=variations.begin() ; iter!=variations.end() ; ++iter) {
817 y[j++] = (*iter)->interpolate(GetX()[i]) ;
818}
819
820 if (!approxGauss) {
821 // Construct central 68% interval from variations collected at each point
822 double pvalue = TMath::Erfc(Z/sqrt(2.)) ;
823 Int_t delta = Int_t( y.size()*(pvalue)/2 + 0.5) ;
824 sort(y.begin(),y.end()) ;
825 lo = y[delta] ;
826 hi = y[y.size()-delta] ;
827 } else {
828 // Estimate R.M.S of variations at each point and use that as Gaussian sigma
829 double sum_y(0);
830 double sum_ysq(0);
831 for (unsigned int k=0 ; k<y.size() ; k++) {
832 sum_y += y[k] ;
833 sum_ysq += y[k]*y[k] ;
834 }
835 sum_y /= y.size() ;
836 sum_ysq /= y.size() ;
837
838 double rms = sqrt(sum_ysq - (sum_y*sum_y)) ;
839 lo = GetY()[i] - Z*rms ;
840 hi = GetY()[i] + Z*rms ;
841 }
842}
843
844
845
846
847////////////////////////////////////////////////////////////////////////////////
848/// Return true if curve is identical to other curve allowing for given
849/// absolute tolerance on each point compared point.
850
851bool RooCurve::isIdentical(const RooCurve& other, double tol, bool verbose) const
852{
853 // Determine X range and Y range
854 Int_t n= min(GetN(),other.GetN());
855 double xmin(1e30);
856 double xmax(-1e30);
857 double ymin(1e30);
858 double ymax(-1e30);
859 for(Int_t i= 0; i < n; i++) {
860 if (fX[i]<xmin) xmin=fX[i] ;
861 if (fX[i]>xmax) xmax=fX[i] ;
862 if (fY[i]<ymin) ymin=fY[i] ;
863 if (fY[i]>ymax) ymax=fY[i] ;
864 }
865 const double Yrange=ymax-ymin ;
866
867 bool ret(true) ;
868 for(Int_t i= 2; i < n-2; i++) {
869 double yTest = interpolate(other.fX[i],1e-10) ;
870 double rdy = std::abs(yTest-other.fY[i])/Yrange ;
871 if (rdy>tol) {
872 ret = false;
873 if(!verbose) continue;
874 cout << "RooCurve::isIdentical[" << std::setw(3) << i << "] Y tolerance exceeded (" << std::setprecision(5) << std::setw(10) << rdy << ">" << tol << "),";
875 cout << " x,y=(" << std::right << std::setw(10) << fX[i] << "," << std::setw(10) << fY[i] << ")\tref: y="
876 << std::setw(10) << other.interpolate(fX[i], 1.E-15) << ". [Nearest point from ref: ";
877 auto j = other.findPoint(fX[i], 1.E10);
878 std::cout << "j=" << j << "\tx,y=(" << std::setw(10) << other.fX[j] << "," << std::setw(10) << other.fY[j] << ") ]" << "\trange=" << Yrange << std::endl;
879 }
880 }
881
882 return ret ;
883}
884
885
886
887////////////////////////////////////////////////////////////////////////////////
888/// Returns sampling hints for a histogram with given boundaries. This helper
889/// function is meant to be used by binned RooAbsReals to produce sampling
890/// hints that are working well with RooFits plotting.
891
892std::list<double> *
893RooCurve::plotSamplingHintForBinBoundaries(std::span<const double> boundaries, double xlo, double xhi)
894{
895 auto hint = new std::list<double>;
896
897 // Make sure the difference between two points around a bin boundary is
898 // larger than the relative epsilon for which the RooCurve considers two
899 // points as the same. Otherwise, the points right of the bin boundary would
900 // be skipped.
901 const double delta = (xhi - xlo) * RooCurve::relativeXEpsilon();
902
903 // Sample points right next to the plot limits
904 hint->push_back(xlo + delta);
905 hint->push_back(xhi - delta);
906
907 // Sample points very close to the left and right of the bin boundaries that
908 // are strictly in between the plot limits.
909 for (const double x : boundaries) {
910 if (x - xlo > delta && xhi - x > delta) {
911 hint->push_back(x - delta);
912 hint->push_back(x + delta);
913 }
914 }
915
916 hint->sort();
917
918 return hint;
919}
#define f(i)
Definition RSha256.hxx:104
#define e(i)
Definition RSha256.hxx:103
#define coutW(a)
#define ccoutP(a)
#define coutE(a)
#define ccoutW(a)
int Int_t
Definition RtypesCore.h:45
#define ClassImp(name)
Definition Rtypes.h:377
@ kCyan
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
static void indent(ostringstream &buf, int indent_level)
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t SetLineWidth
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t SetLineColor
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
float xmin
#define hi
float ymin
float xmax
float ymax
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition Util.h:122
T Sum() const
Definition Util.h:240
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
bool isValid() const
Definition RooAbsFunc.h:37
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).
static Int_t numEvalErrors()
Return the number of logged evaluation errors since the last clearing.
static void printEvalErrors(std::ostream &os=std::cout, Int_t maxPerNode=10000000)
Print all outstanding logged evaluation error on the given ostream.
static void clearEvalErrorLog()
Clear the stack of evaluation error messages.
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
void addPoints(const RooAbsFunc &func, double xlo, double xhi, Int_t minPoints, double prec, double resolution, WingMode wmode, Int_t numee=0, bool doEEVal=false, double eeVal=0.0, std::list< double > *samplingHint=nullptr)
Add points calculated with the specified function, over the range (xlo,xhi).
Definition RooCurve.cxx:284
void printTitle(std::ostream &os) const override
Print the title of this curve.
Definition RooCurve.cxx:498
void initialize()
Perform initialization that is common to all curves.
Definition RooCurve.cxx:241
~RooCurve() override
double getFitRangeBinW() const override
Get the bin width associated with this plotable object.
Definition RooCurve.cxx:479
static constexpr double relativeXEpsilon()
The distance between two points x1 and x2 relative to the full plot range below which two points are ...
Definition RooCurve.h:104
void addRange(const RooAbsFunc &func, double x1, double x2, double y1, double y2, double minDy, double minDx, int numee, bool doEEVal, double eeVal, double epsilon)
Fill the range (x1,x2) with points calculated using func(&x).
Definition RooCurve.cxx:402
void printName(std::ostream &os) const override
Print name of object.
Definition RooCurve.cxx:487
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Print the details of this curve.
Definition RooCurve.cxx:517
double interpolate(double x, double tolerance=1e-10) const
Return linearly interpolated value of curve at xvalue.
Definition RooCurve.cxx:655
static std::list< double > * plotSamplingHintForBinBoundaries(std::span< const double > boundaries, double xlo, double xhi)
Returns sampling hints for a histogram with given boundaries.
Definition RooCurve.cxx:893
void printClassName(std::ostream &os) const override
Print the class name of this curve.
Definition RooCurve.cxx:507
@ Extended
Definition RooCurve.h:39
@ Straight
Definition RooCurve.h:39
RooCurve()
Default constructor.
Definition RooCurve.cxx:84
bool _showProgress
! Show progress indication when adding points
Definition RooCurve.h:97
RooCurve * makeErrorBand(const std::vector< RooCurve * > &variations, double Z=1) const
Construct filled RooCurve represented error band that captures alpha% of the variations of the curves...
Definition RooCurve.cxx:701
double chiSquare(const RooHist &hist, int nFitParam) const
Calculate the chi^2/NDOF of this curve with respect to the histogram 'hist' accounting nFitParam floa...
Definition RooCurve.cxx:535
void calcBandInterval(const std::vector< RooCurve * > &variations, Int_t i, double Z, double &lo, double &hi, bool approxGauss) const
Definition RooCurve.cxx:812
double getFitRangeNEvt() const override
Return the number of events associated with the plotable object, it is always 1 for curves.
Definition RooCurve.cxx:460
void shiftCurveToZero()
Find lowest point in curve and move all points in curve so that lowest point will go exactly through ...
Definition RooCurve.cxx:255
Add a point with the specified coordinates. Update our y-axis limits.
Definition RooCurve.cxx:447
bool isIdentical(const RooCurve &other, double tol=1e-6, bool verbose=true) const
Return true if curve is identical to other curve allowing for given absolute tolerance on each point ...
Definition RooCurve.cxx:851
Int_t findPoint(double value, double tolerance=1e-10) const
Find the nearest point to xvalue.
Definition RooCurve.cxx:633
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
Graphical representation of binned data based on the TGraphAsymmErrors class.
Definition RooHist.h:29
static constexpr double infinity()
Return internal infinity representation.
Definition RooNumber.h:25
void updateYAxisLimits(double y)
Definition RooPlotable.h:30
void setYAxisLimits(double ymin, double ymax)
Definition RooPlotable.h:34
void setYAxisLabel(const char *label)
Definition RooPlotable.h:29
Represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:37
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
virtual void SetBinLabel(Int_t bin, const char *label)
Set label for bin.
Definition TAxis.cxx:886
Bool_t IsAlphanumeric() const
Definition TAxis.h:88
Double_t GetXmax() const
Definition TAxis.h:140
const char * GetBinLabel(Int_t bin) const
Return label for bin.
Definition TAxis.cxx:440
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition TAxis.cxx:794
Double_t GetXmin() const
Definition TAxis.h:139
Int_t GetNbins() const
Definition TAxis.h:125
Double_t * GetEXlow() const override
Double_t * GetEYhigh() const override
Double_t * GetEXhigh() const override
Double_t * GetEYlow() const override
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
virtual Double_t GetPointX(Int_t i) const
Get x value for point i.
Definition TGraph.cxx:1544
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
virtual void Sort(Bool_t(*greater)(const TGraph *, Int_t, Int_t)=&TGraph::CompareX, Bool_t ascending=kTRUE, Int_t low=0, Int_t high=-1111)
Sorts the points of this TGraph using in-place quicksort (see e.g.
Definition TGraph.cxx:2486
Double_t * GetX() const
Definition TGraph.h:139
void SetName(const char *name="") override
Set graph name.
Definition TGraph.cxx:2381
TAxis * GetXaxis() const
Get x axis of the graph.
Definition TGraph.cxx:1566
Double_t * fX
[fNpoints] array of X points
Definition TGraph.h:47
virtual Double_t GetPointY(Int_t i) const
Get y value for point i.
Definition TGraph.cxx:1555
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
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
const char * Data() const
Definition TString.h:376
TString & Append(const char *cs)
Definition TString.h:572
Definition event_demo.C:69
RooConstVar & RooConst(double val)
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TGraphErrors * gr
Definition legend1.C:25
return c2
Definition legend2.C:14
#define F(x, y, z)
Double_t Erfc(Double_t x)
Computes the complementary error function erfc(x).
Definition TMath.cxx:199
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345