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