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