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