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 "RooFit.h"
40
41#include "RooCurve.h"
42#include "RooHist.h"
43#include "RooAbsReal.h"
44#include "RooArgSet.h"
45#include "RooRealVar.h"
46#include "RooRealIntegral.h"
47#include "RooRealBinding.h"
48#include "RooScaledFunc.h"
49#include "RooMsgService.h"
50
51#include "Riostream.h"
52#include "TClass.h"
53#include "TMath.h"
54#include "TAxis.h"
55#include "TMatrixD.h"
56#include "TVectorD.h"
57#include <iomanip>
58#include <deque>
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] scale1 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 Int_t step;
330 Double_t dx= (xhi-xlo)/(minPoints-1.);
331 Double_t *yval= new Double_t[minPoints];
332
333 // Get list of initial x values. If function provides sampling hint use that,
334 // otherwise use default binning of frame
335 list<Double_t>* xval = samplingHint ;
336 if (!xval) {
337 xval = new list<Double_t> ;
338 for(step= 0; step < minPoints; step++) {
339 xval->push_back(xlo + step*dx) ;
340 }
341 }
342
343
344 Double_t ymax(-1e30), ymin(1e30) ;
345
346 step=0 ;
347 for(list<Double_t>::iterator iter = xval->begin() ; iter!=xval->end() ; ++iter,++step) {
348 Double_t xx = *iter ;
349 if (step==minPoints-1) xx-=1e-15 ;
350
351 yval[step]= func(&xx);
352 if (_showProgress) {
353 ccoutP(Plotting) << "." ;
354 cout.flush() ;
355 }
356
358 if (numee>=0) {
359 coutW(Plotting) << "At observable [x]=" << xx << " " ;
361 }
362 if (doEEVal) {
363 yval[step]=eeVal ;
364 }
365 }
367
368
369 if (yval[step]>ymax) ymax=yval[step] ;
370 if (yval[step]<ymin) ymin=yval[step] ;
371 }
372 Double_t yrangeEst=(ymax-ymin) ;
373
374 // store points of the coarse scan and calculate any refinements necessary
375 Double_t minDx= resolution*(xhi-xlo);
376 Double_t x1,x2= xlo;
377
378 if (wmode==Extended) {
379 // Add two points to make curve jump from 0 to yval at the left end of the plotting range.
380 // This ensures that filled polygons are drawn properly.
381 addPoint(xlo-dx*1.001, 0);
382 addPoint(xlo-dx,yval[0]) ;
383 } else if (wmode==Straight) {
384 addPoint(xlo,0) ;
385 }
386
387 addPoint(xlo,yval[0]);
388
389 list<Double_t>::iterator iter2 = xval->begin() ;
390 x1 = *iter2 ;
391 step=1 ;
392 while(true) {
393 x1= x2;
394 ++iter2 ;
395 if (iter2==xval->end()) {
396 break ;
397 }
398 x2= *iter2 ;
399 if (prec<0) {
400 // If precision is <0, no attempt at recursive interpolation is made
401 addPoint(x2,yval[step]) ;
402 } else {
403 addRange(func,x1,x2,yval[step-1],yval[step],prec*yrangeEst,minDx,numee,doEEVal,eeVal);
404 }
405 step++ ;
406 }
407 addPoint(xhi,yval[minPoints-1]) ;
408
409 if (wmode==Extended) {
410 addPoint(xhi+dx,yval[minPoints-1]) ;
411 addPoint(xhi+dx*1.001, 0);
412 } else if (wmode==Straight) {
413 addPoint(xhi,0) ;
414 }
415
416 // cleanup
417 delete [] yval;
418 if (xval != samplingHint) {
419 delete xval ;
420 }
421
422}
423
424
425////////////////////////////////////////////////////////////////////////////////
426/// Fill the range (x1,x2) with points calculated using func(&x). No point will
427/// be added at x1, and a point will always be added at x2. The density of points
428/// will be calculated so that the maximum deviation from a straight line
429/// approximation is prec*(ymax-ymin) down to the specified minimum horizontal spacing.
430
432 Double_t y1, Double_t y2, Double_t minDy, Double_t minDx,
433 Int_t numee, Bool_t doEEVal, Double_t eeVal)
434{
435 // Explicitly skip empty ranges to eliminate point duplication
436 if (fabs(x2-x1)<1e-20) {
437 return ;
438 }
439
440 // calculate our value at the midpoint of this range
441 Double_t xmid= 0.5*(x1+x2);
442 Double_t ymid= func(&xmid);
443 if (_showProgress) {
444 ccoutP(Plotting) << "." ;
445 cout.flush() ;
446 }
447
449 if (numee>=0) {
450 coutW(Plotting) << "At observable [x]=" << xmid << " " ;
452 }
453 if (doEEVal) {
454 ymid=eeVal ;
455 }
456 }
458
459 // test if the midpoint is sufficiently close to a straight line across this interval
460 Double_t dy= ymid - 0.5*(y1+y2);
461 if((xmid - x1 >= minDx) && fabs(dy)>0 && fabs(dy) >= minDy) {
462 // fill in each subrange
463 addRange(func,x1,xmid,y1,ymid,minDy,minDx,numee,doEEVal,eeVal);
464 addRange(func,xmid,x2,ymid,y2,minDy,minDx,numee,doEEVal,eeVal);
465 }
466 else {
467 // add the endpoint
468 addPoint(x2,y2);
469 }
470}
471
472
473////////////////////////////////////////////////////////////////////////////////
474/// Add a point with the specified coordinates. Update our y-axis limits.
475
477{
478// cout << "RooCurve("<< GetName() << ") adding point at (" << x << "," << y << ")" << endl ;
479 Int_t next= GetN();
480 SetPoint(next, x, y);
482}
483
484
485////////////////////////////////////////////////////////////////////////////////
486/// Return the number of events associated with the plotable object,
487/// it is always 1 for curves
488
490 return 1;
491}
492
493
494////////////////////////////////////////////////////////////////////////////////
495/// Return the number of events associated with the plotable object,
496/// in the given range. It is always 1 for curves
497
499{
500 return 1 ;
501}
502
503
504////////////////////////////////////////////////////////////////////////////////
505/// Get the bin width associated with this plotable object.
506/// It is alwats zero for curves
507
509 return 0 ;
510}
511
512
513
514////////////////////////////////////////////////////////////////////////////////
515
516void RooCurve::printName(ostream& os) const
517//
518{
519 // Print the name of this curve
520 os << GetName() ;
521}
522
523
524////////////////////////////////////////////////////////////////////////////////
525/// Print the title of this curve
526
527void RooCurve::printTitle(ostream& os) const
528{
529 os << GetTitle() ;
530}
531
532
533////////////////////////////////////////////////////////////////////////////////
534/// Print the class name of this curve
535
536void RooCurve::printClassName(ostream& os) const
537{
538 os << IsA()->GetName() ;
539}
540
541
542
543////////////////////////////////////////////////////////////////////////////////
544/// Print the details of this curve
545
546void RooCurve::printMultiline(ostream& os, Int_t /*contents*/, Bool_t /*verbose*/, TString indent) const
547{
548 os << indent << "--- RooCurve ---" << endl ;
549 Int_t n= GetN();
550 os << indent << " Contains " << n << " points" << endl;
551 os << indent << " Graph points:" << endl;
552 for(Int_t i= 0; i < n; i++) {
553 os << indent << setw(3) << i << ") x = " << fX[i] << " , y = " << fY[i] << endl;
554 }
555}
556
557
558
559////////////////////////////////////////////////////////////////////////////////
560/// Calculate the chi^2/NDOF of this curve with respect to the histogram
561/// 'hist' accounting nFitParam floating parameters in case the curve
562/// was the result of a fit
563
564Double_t RooCurve::chiSquare(const RooHist& hist, Int_t nFitParam) const
565{
566 Int_t i,np = hist.GetN() ;
567 Double_t x,y,eyl,eyh,exl,exh ;
568
569 // Find starting and ending bin of histogram based on range of RooCurve
570 Double_t xstart,xstop ;
571
572 GetPoint(0,xstart,y) ;
573 GetPoint(GetN()-1,xstop,y) ;
574
575 Int_t nbin(0) ;
576
577 Double_t chisq(0) ;
578 for (i=0 ; i<np ; i++) {
579
580 // Retrieve histogram contents
581 ((RooHist&)hist).GetPoint(i,x,y) ;
582
583 // Check if point is in range of curve
584 if (x<xstart || x>xstop) continue ;
585
586 eyl = hist.GetEYlow()[i] ;
587 eyh = hist.GetEYhigh()[i] ;
588 exl = hist.GetEXlow()[i] ;
589 exh = hist.GetEXhigh()[i] ;
590
591 // Integrate function over this bin
592 Double_t avg = average(x-exl,x+exh) ;
593
594 // Add pull^2 to chisq
595 if (y!=0) {
596 Double_t pull = (y>avg) ? ((y-avg)/eyl) : ((y-avg)/eyh) ;
597 chisq += pull*pull ;
598 nbin++ ;
599 }
600 }
601
602 // Return chisq/nDOF
603 return chisq / (nbin-nFitParam) ;
604}
605
606
607
608////////////////////////////////////////////////////////////////////////////////
609/// Return average curve value in [xFirst,xLast] by integrating curve between points
610/// and dividing by xLast-xFirst
611
613{
614 if (xFirst>=xLast) {
615 coutE(InputArguments) << "RooCurve::average(" << GetName()
616 << ") invalid range (" << xFirst << "," << xLast << ")" << endl ;
617 return 0 ;
618 }
619
620 // Find Y values and begin and end points
621 Double_t yFirst = interpolate(xFirst,1e-10) ;
622 Double_t yLast = interpolate(xLast,1e-10) ;
623
624 // Find first and last mid points
625 Int_t ifirst = findPoint(xFirst,1e10) ;
626 Int_t ilast = findPoint(xLast,1e10) ;
627 Double_t xFirstPt,yFirstPt,xLastPt,yLastPt ;
628 const_cast<RooCurve&>(*this).GetPoint(ifirst,xFirstPt,yFirstPt) ;
629 const_cast<RooCurve&>(*this).GetPoint(ilast,xLastPt,yLastPt) ;
630
631 Double_t tolerance=1e-3*(xLast-xFirst) ;
632
633 // Handle trivial scenario -- no midway points, point only at or outside given range
634 if (ilast-ifirst==1 &&(xFirstPt-xFirst)<-1*tolerance && (xLastPt-xLast)>tolerance) {
635 return 0.5*(yFirst+yLast) ;
636 }
637
638 // If first point closest to xFirst is at xFirst or before xFirst take the next point
639 // as the first midway point
640 if ((xFirstPt-xFirst)<-1*tolerance) {
641 ifirst++ ;
642 const_cast<RooCurve&>(*this).GetPoint(ifirst,xFirstPt,yFirstPt) ;
643 }
644
645 // If last point closest to yLast is at yLast or beyond yLast the the previous point
646 // as the last midway point
647 if ((xLastPt-xLast)>tolerance) {
648 ilast-- ;
649 const_cast<RooCurve&>(*this).GetPoint(ilast,xLastPt,yLastPt) ;
650 }
651
652 Double_t sum(0),x1,y1,x2,y2 ;
653
654 // Trapezoid integration from lower edge to first midpoint
655 sum += (xFirstPt-xFirst)*(yFirst+yFirstPt)/2 ;
656
657 // Trapezoid integration between midpoints
658 Int_t i ;
659 for (i=ifirst ; i<ilast ; i++) {
660 const_cast<RooCurve&>(*this).GetPoint(i,x1,y1) ;
661 const_cast<RooCurve&>(*this).GetPoint(i+1,x2,y2) ;
662 sum += (x2-x1)*(y1+y2)/2 ;
663 }
664
665 // Trapezoid integration from last midpoint to upper edge
666 sum += (xLast-xLastPt)*(yLastPt+yLast)/2 ;
667 return sum/(xLast-xFirst) ;
668}
669
670
671
672////////////////////////////////////////////////////////////////////////////////
673/// Find the nearest point to xvalue. Return -1 if distance
674/// exceeds tolerance
675
677{
678 Double_t delta(std::numeric_limits<double>::max()),x,y ;
679 Int_t i,n = GetN() ;
680 Int_t ibest(-1) ;
681 for (i=0 ; i<n ; i++) {
682 ((RooCurve&)*this).GetPoint(i,x,y) ;
683 if (fabs(xvalue-x)<delta) {
684 delta = fabs(xvalue-x) ;
685 ibest = i ;
686 }
687 }
688
689 return (delta<tolerance)?ibest:-1 ;
690}
691
692
693////////////////////////////////////////////////////////////////////////////////
694/// Return linearly interpolated value of curve at xvalue. If distance
695/// to nearest point is less than tolerance, return nearest point value
696/// instead
697
699{
700 // Find best point
701 int n = GetN() ;
702 int ibest = findPoint(xvalue,1e10) ;
703
704 // Get position of best point
705 Double_t xbest, ybest ;
706 const_cast<RooCurve*>(this)->GetPoint(ibest,xbest,ybest) ;
707
708 // Handle trivial case of being dead on
709 if (fabs(xbest-xvalue)<tolerance) {
710 return ybest ;
711 }
712
713 // Get nearest point on other side w.r.t. xvalue
714 Double_t xother,yother, retVal(0) ;
715 if (xbest<xvalue) {
716 if (ibest==n-1) {
717 // Value beyond end requested -- return value of last point
718 return ybest ;
719 }
720 const_cast<RooCurve*>(this)->GetPoint(ibest+1,xother,yother) ;
721 if (xother==xbest) return ybest ;
722 retVal = ybest + (yother-ybest)*(xvalue-xbest)/(xother-xbest) ;
723
724 } else {
725 if (ibest==0) {
726 // Value before 1st point requested -- return value of 1st point
727 return ybest ;
728 }
729 const_cast<RooCurve*>(this)->GetPoint(ibest-1,xother,yother) ;
730 if (xother==xbest) return ybest ;
731 retVal = yother + (ybest-yother)*(xvalue-xother)/(xbest-xother) ;
732 }
733
734 return retVal ;
735}
736
737
738
739
740////////////////////////////////////////////////////////////////////////////////
741/// Construct filled RooCurve represented error band that captures alpha% of the variations
742/// of the curves passed through argument variations, where the percentage alpha corresponds to
743/// the central interval fraction of a significance Z
744
745RooCurve* RooCurve::makeErrorBand(const vector<RooCurve*>& variations, Double_t Z) const
746{
747 RooCurve* band = new RooCurve ;
748 band->SetName(Form("%s_errorband",GetName())) ;
749 band->SetLineWidth(1) ;
750 band->SetFillColor(kCyan) ;
751 band->SetLineColor(kCyan) ;
752
753 vector<double> bandLo(GetN()) ;
754 vector<double> bandHi(GetN()) ;
755 for (int i=0 ; i<GetN() ; i++) {
756 calcBandInterval(variations,i,Z,bandLo[i],bandHi[i],kFALSE) ;
757 }
758
759 for (int i=0 ; i<GetN() ; i++) {
760 band->addPoint(GetX()[i],bandLo[i]) ;
761 }
762 for (int i=GetN()-1 ; i>=0 ; i--) {
763 band->addPoint(GetX()[i],bandHi[i]) ;
764 }
765 // if the axis of the old graph is alphanumeric, copy the labels to the new one as well
766 if(this->GetXaxis() && this->GetXaxis()->IsAlphanumeric()){
767 band->GetXaxis()->Set(this->GetXaxis()->GetNbins(),this->GetXaxis()->GetXmin(),this->GetXaxis()->GetXmax());
768 for(int i=0; i<this->GetXaxis()->GetNbins(); ++i){
769 band->GetXaxis()->SetBinLabel(i+1,this->GetXaxis()->GetBinLabel(i+1));
770 }
771 }
772
773 return band ;
774}
775
776
777
778
779////////////////////////////////////////////////////////////////////////////////
780/// Construct filled RooCurve represented error band represent the error added in quadrature defined by the curves arguments
781/// plusVar and minusVar corresponding to one-sigma variations of each parameter. The resulting error band, combined used the correlation matrix C
782/// is multiplied with the significance parameter Z to construct the equivalent of a Z sigma error band (in Gaussian approximation)
783
784RooCurve* RooCurve::makeErrorBand(const vector<RooCurve*>& plusVar, const vector<RooCurve*>& minusVar, const TMatrixD& C, Double_t Z) const
785{
786
787 RooCurve* band = new RooCurve ;
788 band->SetName(Form("%s_errorband",GetName())) ;
789 band->SetLineWidth(1) ;
790 band->SetFillColor(kCyan) ;
791 band->SetLineColor(kCyan) ;
792
793 vector<double> bandLo(GetN()) ;
794 vector<double> bandHi(GetN()) ;
795 for (int i=0 ; i<GetN() ; i++) {
796 calcBandInterval(plusVar,minusVar,i,C,Z,bandLo[i],bandHi[i]) ;
797 }
798
799 for (int i=0 ; i<GetN() ; i++) {
800 band->addPoint(GetX()[i],bandLo[i]) ;
801 }
802 for (int i=GetN()-1 ; i>=0 ; i--) {
803 band->addPoint(GetX()[i],bandHi[i]) ;
804 }
805
806 // if the axis of the old graph is alphanumeric, copy the labels to the new one as well
807 if(this->GetXaxis() && this->GetXaxis()->IsAlphanumeric()){
808 band->GetXaxis()->Set(this->GetXaxis()->GetNbins(),this->GetXaxis()->GetXmin(),this->GetXaxis()->GetXmax());
809 for(int i=0; i<this->GetXaxis()->GetNbins(); ++i){
810 band->GetXaxis()->SetBinLabel(i+1,this->GetXaxis()->GetBinLabel(i+1));
811 }
812 }
813
814 return band ;
815}
816
817
818
819
820
821////////////////////////////////////////////////////////////////////////////////
822/// Retrieve variation points from curves
823
824void 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
825{
826 vector<double> y_plus(plusVar.size()), y_minus(minusVar.size()) ;
827 Int_t j(0) ;
828 for (vector<RooCurve*>::const_iterator iter=plusVar.begin() ; iter!=plusVar.end() ; ++iter) {
829 y_plus[j++] = (*iter)->interpolate(GetX()[i]) ;
830 }
831 j=0 ;
832 for (vector<RooCurve*>::const_iterator iter=minusVar.begin() ; iter!=minusVar.end() ; ++iter) {
833 y_minus[j++] = (*iter)->interpolate(GetX()[i]) ;
834 }
835 Double_t y_cen = GetY()[i] ;
836 Int_t n = j ;
837
838 // Make vector of variations
839 TVectorD F(plusVar.size()) ;
840 for (j=0 ; j<n ; j++) {
841 F[j] = (y_plus[j]-y_minus[j])/2 ;
842 }
843
844 // Calculate error in linear approximation from variations and correlation coefficient
845 Double_t sum = F*(C*F) ;
846
847 lo= y_cen + sqrt(sum) ;
848 hi= y_cen - sqrt(sum) ;
849}
850
851
852
853////////////////////////////////////////////////////////////////////////////////
854
855void RooCurve::calcBandInterval(const vector<RooCurve*>& variations,Int_t i,Double_t Z, Double_t& lo, Double_t& hi, Bool_t approxGauss) const
856{
857 vector<double> y(variations.size()) ;
858 Int_t j(0) ;
859 for (vector<RooCurve*>::const_iterator iter=variations.begin() ; iter!=variations.end() ; ++iter) {
860 y[j++] = (*iter)->interpolate(GetX()[i]) ;
861}
862
863 if (!approxGauss) {
864 // Construct central 68% interval from variations collected at each point
865 Double_t pvalue = TMath::Erfc(Z/sqrt(2.)) ;
866 Int_t delta = Int_t( y.size()*(pvalue)/2 + 0.5) ;
867 sort(y.begin(),y.end()) ;
868 lo = y[delta] ;
869 hi = y[y.size()-delta] ;
870 } else {
871 // Estimate R.M.S of variations at each point and use that as Gaussian sigma
872 Double_t sum_y(0), sum_ysq(0) ;
873 for (unsigned int k=0 ; k<y.size() ; k++) {
874 sum_y += y[k] ;
875 sum_ysq += y[k]*y[k] ;
876 }
877 sum_y /= y.size() ;
878 sum_ysq /= y.size() ;
879
880 Double_t rms = sqrt(sum_ysq - (sum_y*sum_y)) ;
881 lo = GetY()[i] - Z*rms ;
882 hi = GetY()[i] + Z*rms ;
883 }
884}
885
886
887
888
889////////////////////////////////////////////////////////////////////////////////
890/// Return true if curve is identical to other curve allowing for given
891/// absolute tolerance on each point compared point.
892
894{
895 // Determine X range and Y range
896 Int_t n= min(GetN(),other.GetN());
897 Double_t xmin(1e30), xmax(-1e30), ymin(1e30), ymax(-1e30) ;
898 for(Int_t i= 0; i < n; i++) {
899 if (fX[i]<xmin) xmin=fX[i] ;
900 if (fX[i]>xmax) xmax=fX[i] ;
901 if (fY[i]<ymin) ymin=fY[i] ;
902 if (fY[i]>ymax) ymax=fY[i] ;
903 }
904 Double_t Yrange=ymax-ymin ;
905
906 Bool_t ret(kTRUE) ;
907 for(Int_t i= 2; i < n-2; i++) {
908 Double_t yTest = interpolate(other.fX[i],1e-10) ;
909 Double_t rdy = fabs(yTest-other.fY[i])/Yrange ;
910 if (rdy>tol) {
911
912// cout << "xref = " << other.fX[i] << " yref = " << other.fY[i] << " xtest = " << fX[i] << " ytest = " << fY[i]
913// << " ytestInt[other.fX] = " << interpolate(other.fX[i],1e-10) << endl ;
914
915 cout << "RooCurve::isIdentical[" << i << "] Y tolerance exceeded (" << rdy << ">" << tol
916 << "), X=" << other.fX[i] << "(" << fX[i] << ")" << " Ytest=" << yTest << " Yref=" << other.fY[i] << " range = " << Yrange << endl ;
917 ret=kFALSE ;
918 }
919 }
920
921 return ret ;
922}
923
924
#define f(i)
Definition: RSha256.hxx:104
#define e(i)
Definition: RSha256.hxx:103
static const double x2[5]
static const double x1[5]
#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:43
unsigned int UInt_t
Definition: RtypesCore.h:44
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
@ kCyan
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
static void indent(ostringstream &buf, int indent_level)
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float type_of_call hi(const int &, const int &)
float ymax
Definition: THbookFile.cxx:93
double sqrt(double)
char * Form(const char *fmt,...)
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
Bool_t isValid() const
Definition: RooAbsFunc.h:33
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:60
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:28
A RooCurve is a one-dimensional graphical representation of a real-valued function.
Definition: RooCurve.h:32
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:476
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:564
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:612
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:431
Bool_t _showProgress
Definition: RooCurve.h:91
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
virtual ~RooCurve()
Destructor.
Definition: RooCurve.cxx:248
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:745
virtual void printName(std::ostream &os) const
Print name of object.
Definition: RooCurve.cxx:516
@ Extended
Definition: RooCurve.h:35
@ Straight
Definition: RooCurve.h:35
RooCurve()
Default constructor.
Definition: RooCurve.cxx:68
Bool_t isIdentical(const RooCurve &other, Double_t tol=1e-6) const
Return true if curve is identical to other curve allowing for given absolute tolerance on each point ...
Definition: RooCurve.cxx:893
Double_t getFitRangeNEvt() const
Return the number of events associated with the plotable object, it is always 1 for curves.
Definition: RooCurve.cxx:489
virtual void printClassName(std::ostream &os) const
Print the class name of this curve.
Definition: RooCurve.cxx:536
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Print the details of this curve.
Definition: RooCurve.cxx:546
Double_t interpolate(Double_t x, Double_t tolerance=1e-10) const
Return linearly interpolated value of curve at xvalue.
Definition: RooCurve.cxx:698
virtual void printTitle(std::ostream &os) const
Print the title of this curve.
Definition: RooCurve.cxx:527
Int_t findPoint(Double_t value, Double_t tolerance=1e-10) const
Find the nearest point to xvalue.
Definition: RooCurve.cxx:676
Double_t getFitRangeBinW() const
Get the bin width associated with this plotable object.
Definition: RooCurve.cxx:508
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:855
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:49
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:21
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:820
Double_t GetXmax() const
Definition: TAxis.h:134
const char * GetBinLabel(Int_t bin) const
Return label for bin.
Definition: TAxis.cxx:437
virtual void Set(Int_t nbins, Double_t xmin, Double_t xmax)
Initialize axis with fix bins.
Definition: TAxis.cxx:728
Bool_t IsAlphanumeric()
Definition: TAxis.h:84
Double_t GetXmin() const
Definition: TAxis.h:133
Int_t GetNbins() const
Definition: TAxis.h:121
Double_t * GetEXhigh() const
Double_t * GetEYlow() const
Double_t * GetEXlow() const
Double_t * GetEYhigh() const
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:2269
Double_t * GetY() const
Definition: TGraph.h:131
virtual void SetName(const char *name="")
Set graph name.
Definition: TGraph.cxx:2308
Int_t GetN() const
Definition: TGraph.h:123
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2324
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:2389
Double_t * GetX() const
Definition: TGraph.h:130
TAxis * GetXaxis() const
Get x axis of the graph.
Definition: TGraph.cxx:1626
Double_t * fX
[fNpoints] array of X points
Definition: TGraph.h:47
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:1593
TString fName
Definition: TNamed.h:32
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
TString & Append(const char *cs)
Definition: TString.h:559
return c1
Definition: legend1.C:41
Double_t y[n]
Definition: legend1.C:17
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< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
@ InputArguments
Definition: RooGlobalFunc.h:68
Double_t Erfc(Double_t x)
Compute the complementary error function erfc(x).
Definition: TMath.cxx:194
static long int sum(long int i)
Definition: Factory.cxx:2275