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