Logo ROOT   6.14/05
Reference Guide
RooXYChi2Var.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 //
19 // Class RooXYChi2Var implements a simple chi^2 calculation from a unbinned
20 // dataset with values x,y with errors on y (and optionally on x) and a function.
21 // The function can be either a RooAbsReal, or an extended RooAbsPdf where
22 // the function value is calculated as the probability density times the
23 // expected number of events
24 // The chi^2 is calculated as
25 //
26 // / (Data[y]-) - func \+2
27 // Sum[point] | ------------------ |
28 // \ Data[ErrY] /
29 //
30 //
31 
32 #include "RooFit.h"
33 
34 #include "RooXYChi2Var.h"
35 #include "RooDataSet.h"
36 #include "RooAbsReal.h"
37 
38 #include "Riostream.h"
39 
40 #include "RooRealVar.h"
41 
43 #include "RooAbsDataStore.h"
44 #include "RooRealBinding.h"
45 #include "RooNumIntFactory.h"
46 
47 using namespace std;
48 
50 ;
51 
52 
53 ////////////////////////////////////////////////////////////////////////////////
54 /// coverity[UNINIT_CTOR]
55 
57 {
58  _funcInt = 0 ;
59  _rrvIter = _rrvArgs.createIterator() ;
60 }
61 
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 ///
65 /// RooXYChi2Var constructor with function and X-Y values dataset
66 ///
67 /// An X-Y dataset is a weighted dataset with one or more observables X where the weight is interpreted
68 /// as the Y value and the weight error is interpreted as the Y value error. The weight must have an
69 /// non-zero error defined at each point for the chi^2 calculation to be meaningful.
70 ///
71 /// To store errors associated with the x and y values in a RooDataSet, call RooRealVar::setAttribute("StoreError")
72 /// on each X-type observable for which the error should be stored and add datapoints to the dataset as follows
73 ///
74 /// RooDataSet::add(xset,yval,yerr) where xset is the RooArgSet of x observables (with or without errors) and yval and yerr
75 /// are the Double_t values that correspond to the Y and its error
76 ///
77 
78 RooXYChi2Var::RooXYChi2Var(const char *name, const char* title, RooAbsReal& func, RooDataSet& xydata, Bool_t integrate) :
79  RooAbsOptTestStatistic(name,title,func,xydata,RooArgSet(),0,0,1,RooFit::Interleave,0,0),
80  _extended(kFALSE),
81  _integrate(integrate),
82  _intConfig(*defaultIntegratorConfig()),
83  _funcInt(0)
84 {
85  _extended = kFALSE ;
86  _yvar = 0 ;
87 
88  initialize() ;
89 }
90 
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 ///
94 /// RooXYChi2Var constructor with function and X-Y values dataset
95 ///
96 /// An X-Y dataset is a weighted dataset with one or more observables X where given yvar is interpreted
97 /// as the Y value. The Y variable must have a non-zero error defined at each point for the chi^2 calculation to be meaningful.
98 ///
99 /// To store errors associated with the x and y values in a RooDataSet, call RooRealVar::setAttribute("StoreError")
100 /// on each X-type observable for which the error should be stored and add datapoints to the dataset as follows
101 ///
102 /// RooDataSet::add(xset,yval,yerr) where xset is the RooArgSet of x observables (with or without errors) and yval and yerr
103 /// are the Double_t values that correspond to the Y and its error
104 ///
105 
106 RooXYChi2Var::RooXYChi2Var(const char *name, const char* title, RooAbsReal& func, RooDataSet& xydata, RooRealVar& yvar, Bool_t integrate) :
107  RooAbsOptTestStatistic(name,title,func,xydata,RooArgSet(),0,0,1,RooFit::Interleave,0,0),
108  _extended(kFALSE),
109  _integrate(integrate),
111  _funcInt(0)
112 {
113  _extended = kFALSE ;
114  _yvar = (RooRealVar*) _dataClone->get()->find(yvar.GetName()) ;
115 
116  initialize() ;
117 }
118 
119 
120 ////////////////////////////////////////////////////////////////////////////////
121 ///
122 /// RooXYChi2Var constructor with an extended p.d.f. and X-Y values dataset
123 /// The value of the function that defines the chi^2 in this form is takes as
124 /// the p.d.f. times the expected number of events
125 ///
126 /// An X-Y dataset is a weighted dataset with one or more observables X where the weight is interpreted
127 /// as the Y value and the weight error is interpreted as the Y value error. The weight must have an
128 /// non-zero error defined at each point for the chi^2 calculation to be meaningful.
129 ///
130 /// To store errors associated with the x and y values in a RooDataSet, call RooRealVar::setAttribute("StoreError")
131 /// on each X-type observable for which the error should be stored and add datapoints to the dataset as follows
132 ///
133 /// RooDataSet::add(xset,yval,yerr) where xset is the RooArgSet of x observables (with or without errors) and yval and yerr
134 /// are the Double_t values that correspond to the Y and its error
135 ///
136 
137 RooXYChi2Var::RooXYChi2Var(const char *name, const char* title, RooAbsPdf& extPdf, RooDataSet& xydata, Bool_t integrate) :
138  RooAbsOptTestStatistic(name,title,extPdf,xydata,RooArgSet(),0,0,1,RooFit::Interleave,0,0),
139  _extended(kTRUE),
140  _integrate(integrate),
142  _funcInt(0)
143 {
144  if (!extPdf.canBeExtended()) {
145  throw(string(Form("RooXYChi2Var::ctor(%s) ERROR: Input p.d.f. must be an extendible",GetName()))) ;
146  }
147  _yvar = 0 ;
148  initialize() ;
149 }
150 
151 
152 
153 
154 ////////////////////////////////////////////////////////////////////////////////
155 ///
156 /// RooXYChi2Var constructor with an extended p.d.f. and X-Y values dataset
157 /// The value of the function that defines the chi^2 in this form is takes as
158 /// the p.d.f. times the expected number of events
159 ///
160 /// An X-Y dataset is a weighted dataset with one or more observables X where the weight is interpreted
161 /// as the Y value and the weight error is interpreted as the Y value error. The weight must have an
162 /// non-zero error defined at each point for the chi^2 calculation to be meaningful.
163 ///
164 /// To store errors associated with the x and y values in a RooDataSet, call RooRealVar::setAttribute("StoreError")
165 /// on each X-type observable for which the error should be stored and add datapoints to the dataset as follows
166 ///
167 /// RooDataSet::add(xset,yval,yerr) where xset is the RooArgSet of x observables (with or without errors) and yval and yerr
168 /// are the Double_t values that correspond to the Y and its error
169 ///
170 
171 RooXYChi2Var::RooXYChi2Var(const char *name, const char* title, RooAbsPdf& extPdf, RooDataSet& xydata, RooRealVar& yvar, Bool_t integrate) :
172  RooAbsOptTestStatistic(name,title,extPdf,xydata,RooArgSet(),0,0,1,RooFit::Interleave,0,0),
173  _extended(kTRUE),
174  _integrate(integrate),
176  _funcInt(0)
177 {
178  if (!extPdf.canBeExtended()) {
179  throw(string(Form("RooXYChi2Var::ctor(%s) ERROR: Input p.d.f. must be an extendible",GetName()))) ;
180  }
181  _yvar = (RooRealVar*) _dataClone->get()->find(yvar.GetName()) ;
182  initialize() ;
183 }
184 
185 
186 
187 
188 ////////////////////////////////////////////////////////////////////////////////
189 /// Copy constructor
190 
191 RooXYChi2Var::RooXYChi2Var(const RooXYChi2Var& other, const char* name) :
192  RooAbsOptTestStatistic(other,name),
193  _extended(other._extended),
194  _integrate(other._integrate),
195  _intConfig(other._intConfig),
196  _funcInt(0)
197 {
198  _yvar = other._yvar ? (RooRealVar*) _dataClone->get()->find(other._yvar->GetName()) : 0 ;
199  initialize() ;
200 
201 }
202 
203 
204 
205 
206 ////////////////////////////////////////////////////////////////////////////////
207 /// Common constructor initialization
208 
210 {
212  RooAbsArg* arg ;
213  while((arg=(RooAbsArg*)iter->Next())) {
214  RooRealVar* var = dynamic_cast<RooRealVar*>(arg) ;
215  if (var) {
216  _rrvArgs.add(*var) ;
217  }
218  }
219  if (_yvar) {
220  _rrvArgs.add(*_yvar) ;
221  }
222  delete iter ;
224 
225  // Define alternate numeric integrator configuration for bin integration
226  // We expect bin contents to very only very slowly so a non-adaptive
227  // Gauss-Kronrod integrator is expected to perform well.
228  _intConfig.setEpsRel(1e-7) ;
229  _intConfig.setEpsAbs(1e-7) ;
230  _intConfig.method1D().setLabel("RooGaussKronrodIntegrator1D") ;
231  _intConfig.methodND().setLabel("RooAdaptiveIntegratorND") ;
232 
233  initIntegrator() ;
234 
235 }
236 
237 
238 
239 ////////////////////////////////////////////////////////////////////////////////
240 /// Initialize bin content integrator
241 
243 {
244  if (!_funcInt) {
246  _rrvIter->Reset() ;
247  RooRealVar* x ;
248  while((x=(RooRealVar*)_rrvIter->Next())) {
249  _binList.push_back(&x->getBinning("bin",kFALSE,kTRUE)) ;
250  }
251  }
252 
253 }
254 
255 
256 
257 ////////////////////////////////////////////////////////////////////////////////
258 /// Destructor
259 
261 {
262  delete _rrvIter ;
263  if (_funcInt) delete _funcInt ;
264 }
265 
266 
267 
268 
269 ////////////////////////////////////////////////////////////////////////////////
270 /// Calculate contribution to internal error due to error on 'x' coordinates
271 /// at point i
272 
274 {
275  RooRealVar* var ;
276  Double_t ret(0) ;
277 
278  _rrvIter->Reset() ;
279  while((var=(RooRealVar*)_rrvIter->Next())) {
280 
281  if (var->hasAsymError()) {
282 
283  // Get value at central X
284  Double_t cxval = var->getVal() ;
285  Double_t xerrLo = -var->getAsymErrorLo() ;
286  Double_t xerrHi = var->getAsymErrorHi() ;
287  Double_t xerr = (xerrLo+xerrHi)/2 ;
288 
289  // Get value at X-eps
290  var->setVal(cxval - xerr/100) ;
291  Double_t fxmin = fy() ;
292 
293  // Get value at X+eps
294  var->setVal(cxval + xerr/100) ;
295  Double_t fxmax = fy() ;
296 
297  // Calculate slope
298  Double_t slope = (fxmax-fxmin)/(2*xerr/100.) ;
299 
300 // cout << "xerrHi = " << xerrHi << " xerrLo = " << xerrLo << " slope = " << slope << endl ;
301 
302  // Asymmetric X error, decide which one to use
303  if ((ydata>cxval && fxmax>fxmin) || (ydata<=cxval && fxmax<=fxmin)) {
304  // Use right X error
305  ret += pow(xerrHi*slope,2) ;
306  } else {
307  // Use left X error
308  ret += pow(xerrLo*slope,2) ;
309  }
310 
311  } else if (var->hasError()) {
312 
313  // Get value at central X
314  Double_t cxval = var->getVal() ;
315  Double_t xerr = var->getError() ;
316 
317  // Get value at X-eps
318  var->setVal(cxval - xerr/100) ;
319  Double_t fxmin = fy() ;
320 
321  // Get value at X+eps
322  var->setVal(cxval + xerr/100) ;
323  Double_t fxmax = fy() ;
324 
325  // Calculate slope
326  Double_t slope = (fxmax-fxmin)/(2*xerr/100.) ;
327 
328 // cout << var << " " ;
329 // var->Print() ;
330 
331 // cout << var->GetName() << " xerr = " << xerr << " slope = " << slope << endl ;
332 
333  // Symmetric X error
334  ret += pow(xerr*slope,2) ;
335  }
336  }
337  return ret ;
338 }
339 
340 
341 
342 
343 ////////////////////////////////////////////////////////////////////////////////
344 /// Return function value requested bu present configuration
345 ///
346 /// If integration is required, the function value integrated
347 /// over the bin volume divided by the bin volume is returned,
348 /// otherwise the value at the bin center is returned.
349 /// The bin volume is defined by the error on the 'X' coordinates
350 ///
351 /// If an extended p.d.f. is used as function, its value is
352 /// also multiplied by the expected number of events here
353 
355 {
356  // Get function value
357  Double_t yfunc ;
358  if (!_integrate) {
359  yfunc = _funcClone->getVal(_dataClone->get()) ;
360  } else {
361  Double_t volume(1) ;
362  _rrvIter->Reset() ;
363  for (list<RooAbsBinning*>::const_iterator iter = _binList.begin() ; iter != _binList.end() ; ++iter) {
364  RooRealVar* x = (RooRealVar*) _rrvIter->Next() ;
365  Double_t xmin = x->getVal() + x->getErrorLo() ;
366  Double_t xmax = x->getVal() + x->getErrorHi() ;
367  (*iter)->setRange(xmin,xmax) ;
368  x->setShapeDirty() ;
369  volume *= (xmax - xmin) ;
370  }
371  Double_t ret = _funcInt->getVal() ;
372  return ret / volume ;
373  }
374  if (_extended) {
375  RooAbsPdf* pdf = (RooAbsPdf*) _funcClone ;
376  // Multiply with expected number of events
377  yfunc *= pdf->expectedEvents(_dataClone->get()) ;
378  }
379  return yfunc ;
380 }
381 
382 
383 
384 ////////////////////////////////////////////////////////////////////////////////
385 /// Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize
386 
387 Double_t RooXYChi2Var::evaluatePartition(Int_t firstEvent, Int_t lastEvent, Int_t stepSize) const
388 {
389  Double_t result(0), carry(0);
390 
391  // Loop over bins of dataset
392  RooDataSet* xydata = (RooDataSet*) _dataClone ;
393 
394  _dataClone->store()->recalculateCache( _projDeps, firstEvent, lastEvent, stepSize,kFALSE ) ;
395 
396  for (Int_t i=firstEvent ; i<lastEvent ; i+=stepSize) {
397 
398  // get the data values for this event
399  xydata->get(i);
400 
401  if (!xydata->valid()) {
402  continue ;
403  }
404 
405 // cout << "xydata = " << endl ;
406 // xydata->get()->Print("v") ;
407  //xydata->store()->dump() ;
408 
409  // Get function value
410  Double_t yfunc = fy() ;
411 
412  // Get data value and error
413  Double_t ydata ;
414  Double_t eylo,eyhi ;
415  if (_yvar) {
416  ydata = _yvar->getVal() ;
417  eylo = -1*_yvar->getErrorLo() ;
418  eyhi = _yvar->getErrorHi() ;
419  } else {
420  ydata = xydata->weight() ;
421  xydata->weightError(eylo,eyhi) ;
422  }
423 
424  // Calculate external error
425  Double_t eExt = yfunc-ydata ;
426 
427  // Pick upper or lower error bar depending on sign of external error
428  Double_t eInt = (eExt>0) ? eyhi : eylo ;
429 
430  // Add contributions due to error in x coordinates
431  Double_t eIntX2 = _integrate ? 0 : xErrorContribution(ydata) ;
432 
433 // cout << "fy = " << yfunc << " eExt = " << eExt << " eInt = " << eInt << " eIntX2 = " << eIntX2 << endl ;
434 
435  // Return 0 if eInt=0, special handling in MINUIT will follow
436  if (eInt==0.) {
437  coutE(Eval) << "RooXYChi2Var::RooXYChi2Var(" << GetName() << ") INFINITY ERROR: data point " << i
438  << " has zero error, but function is not zero (f=" << yfunc << ")" << endl ;
439  return 0 ;
440  }
441 
442  // Add chi2 term
443  Double_t term = eExt*eExt/(eInt*eInt+ eIntX2);
444  Double_t y = term - carry;
445  Double_t t = result + y;
446  carry = (t - result) - y;
447  result = t;
448  }
449 
450  _evalCarry = carry;
451  return result ;
452 }
453 
454 
455 
457 {
458  // Inform base class that observable yvar cannot be optimized away from the dataset
459  if (_yvar) return RooArgSet(*_yvar) ;
460  return RooArgSet() ;
461 }
462 
463 
464 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * _rrvIter
Definition: RooXYChi2Var.h:75
TIterator * createIterator(Bool_t dir=kIterForward) const
#define coutE(a)
Definition: RooMsgService.h:34
float xmin
Definition: THbookFile.cxx:93
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:86
virtual void Reset()=0
static RooNumIntConfig * defaultIntegratorConfig()
Returns the default numeric integration configuration for all RooAbsReals.
RooArgSet requiredExtraObservables() const
virtual const RooArgSet * get() const
Definition: RooAbsData.h:79
Bool_t _extended
Definition: RooXYChi2Var.h:70
RooCategory & methodND()
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
const RooAbsBinning & getBinning(const char *name=0, Bool_t verbose=kTRUE, Bool_t createOnTheFly=kFALSE) const
Return binning definition with name.
Definition: RooRealVar.cxx:266
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual Double_t evaluatePartition(Int_t firstEvent, Int_t lastEvent, Int_t stepSize) const
Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize.
STL namespace.
void initialize()
Iterator over set of real-valued observables.
Iterator abstract base class.
Definition: TIterator.h:30
virtual Bool_t setLabel(const char *label, Bool_t printError=kTRUE)
Set value by specifying the name of the desired state If printError is set, a message will be printed...
RooXYChi2Var()
coverity[UNINIT_CTOR]
Double_t getErrorLo() const
Definition: RooRealVar.h:62
RooCategory & method1D()
Double_t x[n]
Definition: legend1.C:17
Bool_t hasAsymError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:59
Double_t _evalCarry
avoids loss of precision
RooRealVar * _yvar
Definition: RooXYChi2Var.h:73
double pow(double, double)
std::list< RooAbsBinning * > _binList
Function integral.
Definition: RooXYChi2Var.h:85
friend class RooArgSet
Definition: RooAbsArg.h:471
Bool_t hasError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:54
virtual void weightError(Double_t &lo, Double_t &hi, ErrorType etype=SumW2) const
Return asymmetric error on weight. (Dummy implementation returning zero)
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events from this p.d.f for use in extended likelihood calculations.
Definition: RooAbsPdf.cxx:2953
virtual Double_t weight() const
Return event weight of current event.
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:204
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:501
RooNumIntConfig _intConfig
Definition: RooXYChi2Var.h:83
RooAbsReal * _funcInt
Definition: RooXYChi2Var.h:84
Double_t getAsymErrorHi() const
Definition: RooRealVar.h:58
Bool_t _integrate
Definition: RooXYChi2Var.h:71
virtual ~RooXYChi2Var()
Destructor.
void initIntegrator()
Initialize bin content integrator.
char * Form(const char *fmt,...)
const UInt_t eInt[256]
float xmax
Definition: THbookFile.cxx:93
RooAbsDataStore * store()
Definition: RooAbsData.h:55
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:216
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event &#39;index&#39;.
const Bool_t kFALSE
Definition: RtypesCore.h:88
#define ClassImp(name)
Definition: Rtypes.h:359
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Bool_t valid() const
Definition: RooAbsData.h:85
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
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t xErrorContribution(Double_t ydata) const
Calculate contribution to internal error due to error on &#39;x&#39; coordinates at point i...
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
void setShapeDirty() const
Definition: RooAbsArg.h:442
virtual TObject * Next()=0
Double_t fy() const
Return function value requested bu present configuration.
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
virtual void recalculateCache(const RooArgSet *, Int_t, Int_t, Int_t, Bool_t)
void setEpsAbs(Double_t newEpsAbs)
Set absolute convergence criteria (convergence if abs(Err)<newEpsAbs)
Double_t getError() const
Definition: RooRealVar.h:53
RooArgSet _rrvArgs
Definition: RooXYChi2Var.h:74
const Bool_t kTRUE
Definition: RtypesCore.h:87
Double_t getAsymErrorLo() const
Definition: RooRealVar.h:57
Double_t getErrorHi() const
Definition: RooRealVar.h:63
char name[80]
Definition: TGX11.cxx:109
void setEpsRel(Double_t newEpsRel)
Set relative convergence criteria (convergence if abs(Err)/abs(Int)<newEpsRel)
RooAbsOptTestStatistic is the abstract base class for test statistics objects that evaluate a functio...