Logo ROOT  
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
42//#include "RooGaussKronrodIntegrator1D.h"
43#include "RooAbsDataStore.h"
44#include "RooRealBinding.h"
45#include "RooNumIntFactory.h"
46
47using namespace std;
48
50;
51
52
53////////////////////////////////////////////////////////////////////////////////
54/// coverity[UNINIT_CTOR]
55
57{
58 _funcInt = 0 ;
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
78RooXYChi2Var::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{
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
106RooXYChi2Var::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),
110 _intConfig(*defaultIntegratorConfig()),
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
137RooXYChi2Var::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),
141 _intConfig(*defaultIntegratorConfig()),
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
171RooXYChi2Var::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),
175 _intConfig(*defaultIntegratorConfig()),
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
191RooXYChi2Var::RooXYChi2Var(const RooXYChi2Var& other, const char* 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 (if RooFitMore is available)
230#ifdef R__HAS_MATHMORE
231 _intConfig.method1D().setLabel("RooGaussKronrodIntegrator1D") ;
232#endif
233 _intConfig.methodND().setLabel("RooAdaptiveIntegratorND") ;
234
236
237}
238
239
240
241////////////////////////////////////////////////////////////////////////////////
242/// Initialize bin content integrator
243
245{
246 if (!_funcInt) {
248 _rrvIter->Reset() ;
249 RooRealVar* x ;
250 while((x=(RooRealVar*)_rrvIter->Next())) {
251 _binList.push_back(&x->getBinning("bin",kFALSE,kTRUE)) ;
252 }
253 }
254
255}
256
257
258
259////////////////////////////////////////////////////////////////////////////////
260/// Destructor
261
263{
264 delete _rrvIter ;
265 if (_funcInt) delete _funcInt ;
266}
267
268
269
270
271////////////////////////////////////////////////////////////////////////////////
272/// Calculate contribution to internal error due to error on 'x' coordinates
273/// at point i
274
276{
277 RooRealVar* var ;
278 Double_t ret(0) ;
279
280 _rrvIter->Reset() ;
281 while((var=(RooRealVar*)_rrvIter->Next())) {
282
283 if (var->hasAsymError()) {
284
285 // Get value at central X
286 Double_t cxval = var->getVal() ;
287 Double_t xerrLo = -var->getAsymErrorLo() ;
288 Double_t xerrHi = var->getAsymErrorHi() ;
289 Double_t xerr = (xerrLo+xerrHi)/2 ;
290
291 // Get value at X-eps
292 var->setVal(cxval - xerr/100) ;
293 Double_t fxmin = fy() ;
294
295 // Get value at X+eps
296 var->setVal(cxval + xerr/100) ;
297 Double_t fxmax = fy() ;
298
299 // Calculate slope
300 Double_t slope = (fxmax-fxmin)/(2*xerr/100.) ;
301
302// cout << "xerrHi = " << xerrHi << " xerrLo = " << xerrLo << " slope = " << slope << endl ;
303
304 // Asymmetric X error, decide which one to use
305 if ((ydata>cxval && fxmax>fxmin) || (ydata<=cxval && fxmax<=fxmin)) {
306 // Use right X error
307 ret += pow(xerrHi*slope,2) ;
308 } else {
309 // Use left X error
310 ret += pow(xerrLo*slope,2) ;
311 }
312
313 } else if (var->hasError()) {
314
315 // Get value at central X
316 Double_t cxval = var->getVal() ;
317 Double_t xerr = var->getError() ;
318
319 // Get value at X-eps
320 var->setVal(cxval - xerr/100) ;
321 Double_t fxmin = fy() ;
322
323 // Get value at X+eps
324 var->setVal(cxval + xerr/100) ;
325 Double_t fxmax = fy() ;
326
327 // Calculate slope
328 Double_t slope = (fxmax-fxmin)/(2*xerr/100.) ;
329
330// cout << var << " " ;
331// var->Print() ;
332
333// cout << var->GetName() << " xerr = " << xerr << " slope = " << slope << endl ;
334
335 // Symmetric X error
336 ret += pow(xerr*slope,2) ;
337 }
338 }
339 return ret ;
340}
341
342
343
344
345////////////////////////////////////////////////////////////////////////////////
346/// Return function value requested bu present configuration
347///
348/// If integration is required, the function value integrated
349/// over the bin volume divided by the bin volume is returned,
350/// otherwise the value at the bin center is returned.
351/// The bin volume is defined by the error on the 'X' coordinates
352///
353/// If an extended p.d.f. is used as function, its value is
354/// also multiplied by the expected number of events here
355
357{
358 // Get function value
359 Double_t yfunc ;
360 if (!_integrate) {
361 yfunc = _funcClone->getVal(_dataClone->get()) ;
362 } else {
363 Double_t volume(1) ;
364 _rrvIter->Reset() ;
365 for (list<RooAbsBinning*>::const_iterator iter = _binList.begin() ; iter != _binList.end() ; ++iter) {
367 Double_t xmin = x->getVal() + x->getErrorLo() ;
368 Double_t xmax = x->getVal() + x->getErrorHi() ;
369 (*iter)->setRange(xmin,xmax) ;
370 x->setShapeDirty() ;
371 volume *= (xmax - xmin) ;
372 }
373 Double_t ret = _funcInt->getVal() ;
374 return ret / volume ;
375 }
376 if (_extended) {
377 RooAbsPdf* pdf = (RooAbsPdf*) _funcClone ;
378 // Multiply with expected number of events
379 yfunc *= pdf->expectedEvents(_dataClone->get()) ;
380 }
381 return yfunc ;
382}
383
384
385
386////////////////////////////////////////////////////////////////////////////////
387/// Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize
388
389Double_t RooXYChi2Var::evaluatePartition(std::size_t firstEvent, std::size_t lastEvent, std::size_t stepSize) const
390{
391 Double_t result(0), carry(0);
392
393 // Loop over bins of dataset
394 RooDataSet* xydata = (RooDataSet*) _dataClone ;
395
396 _dataClone->store()->recalculateCache( _projDeps, firstEvent, lastEvent, stepSize,kFALSE ) ;
397
398 for (auto i=firstEvent ; i<lastEvent ; i+=stepSize) {
399
400 // get the data values for this event
401 xydata->get(i);
402
403 if (!xydata->valid()) {
404 continue ;
405 }
406
407// cout << "xydata = " << endl ;
408// xydata->get()->Print("v") ;
409 //xydata->store()->dump() ;
410
411 // Get function value
412 Double_t yfunc = fy() ;
413
414 // Get data value and error
415 Double_t ydata ;
416 Double_t eylo,eyhi ;
417 if (_yvar) {
418 ydata = _yvar->getVal() ;
419 eylo = -1*_yvar->getErrorLo() ;
420 eyhi = _yvar->getErrorHi() ;
421 } else {
422 ydata = xydata->weight() ;
423 xydata->weightError(eylo,eyhi) ;
424 }
425
426 // Calculate external error
427 Double_t eExt = yfunc-ydata ;
428
429 // Pick upper or lower error bar depending on sign of external error
430 Double_t eInt = (eExt>0) ? eyhi : eylo ;
431
432 // Add contributions due to error in x coordinates
433 Double_t eIntX2 = _integrate ? 0 : xErrorContribution(ydata) ;
434
435// cout << "fy = " << yfunc << " eExt = " << eExt << " eInt = " << eInt << " eIntX2 = " << eIntX2 << endl ;
436
437 // Return 0 if eInt=0, special handling in MINUIT will follow
438 if (eInt==0.) {
439 coutE(Eval) << "RooXYChi2Var::RooXYChi2Var(" << GetName() << ") INFINITY ERROR: data point " << i
440 << " has zero error, but function is not zero (f=" << yfunc << ")" << endl ;
441 return 0 ;
442 }
443
444 // Add chi2 term
445 Double_t term = eExt*eExt/(eInt*eInt+ eIntX2);
446 Double_t y = term - carry;
447 Double_t t = result + y;
448 carry = (t - result) - y;
449 result = t;
450 }
451
452 _evalCarry = carry;
453 return result ;
454}
455
456
457
459{
460 // Inform base class that observable yvar cannot be optimized away from the dataset
461 if (_yvar) return RooArgSet(*_yvar) ;
462 return RooArgSet() ;
463}
#define e(i)
Definition: RSha256.hxx:103
#define coutE(a)
Definition: RooMsgService.h:33
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
double pow(double, double)
char * Form(const char *fmt,...)
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
friend class RooArgSet
Definition: RooAbsArg.h:572
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual void recalculateCache(const RooArgSet *, Int_t, Int_t, Int_t, Bool_t)
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
RooAbsDataStore * store()
Definition: RooAbsData.h:63
virtual Bool_t valid() const
Definition: RooAbsData.h:93
RooAbsOptTestStatistic is the abstract base class for test statistics objects that evaluate a functio...
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
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:3283
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
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:560
Double_t _evalCarry
avoids loss of precision
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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:88
virtual Bool_t setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
virtual void weightError(Double_t &lo, Double_t &hi, ErrorType etype=SumW2) const override
Return asymmetric error on weight. (Dummy implementation returning zero)
Definition: RooDataSet.cxx:998
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
virtual Double_t weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:973
void setEpsAbs(Double_t newEpsAbs)
Set absolute convergence criteria (convergence if abs(Err)<newEpsAbs)
RooCategory & methodND()
void setEpsRel(Double_t newEpsRel)
Set relative convergence criteria (convergence if abs(Err)/abs(Int)<newEpsRel)
RooCategory & method1D()
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
Double_t getAsymErrorLo() const
Definition: RooRealVar.h:66
Double_t getErrorHi() const
Definition: RooRealVar.h:72
Bool_t hasAsymError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:68
Double_t getErrorLo() const
Definition: RooRealVar.h:71
Double_t getAsymErrorHi() const
Definition: RooRealVar.h:67
Bool_t hasError(Bool_t allowZero=kTRUE) const
Definition: RooRealVar.h:63
Double_t getError() const
Definition: RooRealVar.h:62
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:225
RooAbsReal * _funcInt
Definition: RooXYChi2Var.h:84
RooArgSet requiredExtraObservables() const
std::list< RooAbsBinning * > _binList
Function integral.
Definition: RooXYChi2Var.h:85
Bool_t _extended
Definition: RooXYChi2Var.h:70
void initIntegrator()
Initialize bin content integrator.
RooRealVar * _yvar
Definition: RooXYChi2Var.h:73
RooXYChi2Var()
coverity[UNINIT_CTOR]
virtual ~RooXYChi2Var()
Destructor.
TIterator * _rrvIter
Definition: RooXYChi2Var.h:75
RooArgSet _rrvArgs
Definition: RooXYChi2Var.h:74
virtual Double_t evaluatePartition(std::size_t firstEvent, std::size_t lastEvent, std::size_t stepSize) const
Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize.
Double_t fy() const
Return function value requested bu present configuration.
void initialize()
Iterator over set of real-valued observables.
Bool_t _integrate
Definition: RooXYChi2Var.h:71
RooNumIntConfig _intConfig
Definition: RooXYChi2Var.h:83
Double_t xErrorContribution(Double_t ydata) const
Calculate contribution to internal error due to error on 'x' coordinates at point i.
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const UInt_t eInt[256]
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
@ Interleave
Definition: RooGlobalFunc.h:70