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