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