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