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 <iomanip>
58 #include <deque>
59 
60 using namespace std ;
61 
63 
64 
65 ////////////////////////////////////////////////////////////////////////////////
66 /// Default constructor.
67 
68 RooCurve::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(),
87  RooPlotable(),
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++) {
153  updateYAxisLimits(fY[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 
167 RooCurve::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++) {
181  updateYAxisLimits(fY[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 
201 RooCurve::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
230  Double_t last(-RooNumber::infinity()) ;
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 
307 void RooCurve::addPoints(const RooAbsFunc &func, Double_t xlo, Double_t xhi,
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 
357  if (RooAbsReal::numEvalErrors()>0) {
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.00000001,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,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 
448  if (RooAbsReal::numEvalErrors()>0) {
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 
516 void 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 
527 void RooCurve::printTitle(ostream& os) const
528 {
529  os << GetTitle() ;
530 }
531 
532 
533 ////////////////////////////////////////////////////////////////////////////////
534 /// Print the class name of this curve
535 
536 void RooCurve::printClassName(ostream& os) const
537 {
538  os << IsA()->GetName() ;
539 }
540 
541 
542 
543 ////////////////////////////////////////////////////////////////////////////////
544 /// Print the details of this curve
545 
546 void 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 
564 Double_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 
676 Int_t RooCurve::findPoint(Double_t xvalue, Double_t tolerance) const
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 
745 RooCurve* 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 
784 RooCurve* 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 
824 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
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 
855 void 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 
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:527
n
const Int_t n
Definition: legend1.C:16
RooScaledFunc.h
RooCurve::getFitRangeBinW
Double_t getFitRangeBinW() const
Get the bin width associated with this plotable object.
Definition: RooCurve.cxx:508
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:77
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
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:3825
e
#define e(i)
Definition: RSha256.hxx:121
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:612
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:271
TVectorD.h
TGraph::SetTitle
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2324
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:122
RooFit.h
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:68
RooCurve::~RooCurve
virtual ~RooCurve()
Destructor.
Definition: RooCurve.cxx:248
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:676
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:54
xmax
float xmax
Definition: THbookFile.cxx:95
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:516
Int_t
int Int_t
Definition: RtypesCore.h:45
TNamed::fName
TString fName
Definition: TNamed.h:38
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
Definition: RooAbsReal.h:61
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:564
TString
Definition: TString.h:136
TMatrixT
Definition: TMatrixDfwd.h:22
TAxis::IsAlphanumeric
Bool_t IsAlphanumeric()
Definition: TAxis.h:84
RooCurve::Extended
@ Extended
Definition: RooCurve.h:35
RooAbsReal::clearEvalErrorLog
static void clearEvalErrorLog()
Clear the stack of evaluation error messages.
Definition: RooAbsReal.cxx:3805
bool
RooAbsFunc
Definition: RooAbsFunc.h:23
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:1626
kCyan
@ kCyan
Definition: Rtypes.h:66
TGraph::GetX
Double_t * GetX() const
Definition: TGraph.h:130
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:745
SetLineWidth
f1 SetLineWidth(4)
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:546
TGraph::SetName
virtual void SetName(const char *name="")
Set graph name.
Definition: TGraph.cxx:2308
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:2389
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
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:68
RooFit::Plotting
@ Plotting
Definition: RooGlobalFunc.h:67
RooCurve::initialize
void initialize()
Perform initialization that is common to all curves.
Definition: RooCurve.cxx:257
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:489
RooCurve
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:131
sqrt
double sqrt(double)
RooRealVar.h
RooHist
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:2269
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:1593
ymin
float ymin
Definition: THbookFile.cxx:95
TAxis::SetBinLabel
virtual void SetBinLabel(Int_t bin, const char *label)
Set label for bin.
Definition: TAxis.cxx:823
sum
static long int sum(long int i)
Definition: Factory.cxx:2272
TMath::Erfc
Double_t Erfc(Double_t x)
Compute the complementary error function erfc(x).
Definition: TMath.cxx:194
TVectorT
Definition: TMatrixTBase.h:78
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:33
RooCurve::printClassName
virtual void printClassName(std::ostream &os) const
Print the class name of this curve.
Definition: RooCurve.cxx:536
Double_t
double Double_t
Definition: RtypesCore.h:59
TGraph
Definition: TGraph.h:41
RooPlotable::updateYAxisLimits
void updateYAxisLimits(Double_t y)
Definition: RooPlotable.h:33
ccoutW
#define ccoutW(a)
Definition: RooMsgService.h:40
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:698
TAxis.h
TGraph::GetN
Int_t GetN() const
Definition: TGraph.h:123
name
char name[80]
Definition: TGX11.cxx:110
RooScaledFunc
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
Definition: RooPlotable.h:26
RooRealIntegral.h
SetLineColor
lv SetLineColor(kBlue)
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
TGraphAsymmErrors::GetEYlow
Double_t * GetEYlow() const
Definition: TGraphAsymmErrors.h:79
RooAbsReal::numEvalErrors
static Int_t numEvalErrors()
Return the number of logged evaluation errors since the last clearing.
Definition: RooAbsReal.cxx:3867
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:476
TMatrixD.h
Riostream.h
RooAbsRealLValue
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:431
TGraphAsymmErrors::GetEXhigh
Double_t * GetEXhigh() const
Definition: TGraphAsymmErrors.h:78
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:307
RooCurve::_showProgress
Bool_t _showProgress
Definition: RooCurve.h:91
TMath.h
RooArgSet
Definition: RooArgSet.h:28
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:855
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:893
c1
return c1
Definition: legend1.C:41
TGraphAsymmErrors::GetEYhigh
Double_t * GetEYhigh() const
Definition: TGraphAsymmErrors.h:80