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