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