ROOT   Reference Guide
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 *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
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
45using std::endl;
46
47namespace {
48 RooAbsTestStatistic::Configuration makeRooAbsTestStatisticCfg() {
50 cfg.verbose = false;
51 return cfg;
52 }
53}
54
55
56////////////////////////////////////////////////////////////////////////////////
57///
58/// RooXYChi2Var constructor with function and X-Y values dataset
59///
60/// If the function is a pdf, it must be extendable. in this case, hhe value of
61/// the function that defines the chi^2 in this form is takes as the p.d.f.
62/// times the expected number of events
63///
64/// An X-Y dataset is a weighted dataset with one or more observables X where the weight is interpreted
65/// as the Y value and the weight error is interpreted as the Y value error. The weight must have an
66/// non-zero error defined at each point for the chi^2 calculation to be meaningful.
67///
68/// To store errors associated with the x and y values in a RooDataSet, call RooRealVar::setAttribute("StoreError")
69/// on each X-type observable for which the error should be stored and add datapoints to the dataset as follows
70///
71/// RooDataSet::add(xset,yval,yerr) where xset is the RooArgSet of x observables (with or without errors) and yval and yerr
72/// are the double values that correspond to the Y and its error
73///
74
75RooXYChi2Var::RooXYChi2Var(const char *name, const char *title, RooAbsReal &func, RooDataSet &xydata, bool integrate)
76 : RooXYChi2Var{name, title, func, xydata, nullptr, integrate, makeRooAbsTestStatisticCfg()}
77{
78}
79
80////////////////////////////////////////////////////////////////////////////////
81///
82/// RooXYChi2Var constructor with function and X-Y values dataset.
83///
84/// If the function is a pdf, it must be extendable. in this case, hhe value of
85/// the function that defines the chi^2 in this form is takes as the p.d.f.
86/// times the expected number of events
87///
88/// An X-Y dataset is a weighted dataset with one or more observables X where given yvar is interpreted
89/// 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.
90///
91/// To store errors associated with the x and y values in a RooDataSet, call RooRealVar::setAttribute("StoreError")
92/// on each X-type observable for which the error should be stored and add datapoints to the dataset as follows
93///
94/// RooDataSet::add(xset,yval,yerr) where xset is the RooArgSet of x observables (with or without errors) and yval and yerr
95/// are the double values that correspond to the Y and its error
96///
97
98RooXYChi2Var::RooXYChi2Var(const char *name, const char *title, RooAbsReal &func, RooDataSet &xydata, RooRealVar &yvar,
99 bool integrate)
100 : RooXYChi2Var{name, title, func, xydata, &yvar, integrate, makeRooAbsTestStatisticCfg()}
101{
102}
103
104
105/// \cond ROOFIT_INTERNAL
106// For internal use in RooAbsReal::createChi2().
107RooXYChi2Var::RooXYChi2Var(const char *name, const char *title, RooAbsReal &func, RooAbsData &data, RooRealVar *yvar,
108 bool integrate, RooAbsTestStatistic::Configuration const &cfg)
109 : RooAbsOptTestStatistic(name, title, func, data, RooArgSet(), cfg),
110 _integrate(integrate),
111 _intConfig(*defaultIntegratorConfig())
112{
113 bool isPdf = dynamic_cast<RooAbsPdf const *>(&func) != nullptr;
114
115 if (isPdf) {
116 auto &extPdf = static_cast<RooAbsPdf const &>(func);
117 if (!extPdf.canBeExtended()) {
118 throw std::runtime_error(
119 Form("RooXYChi2Var::RooXYChi2Var(%s) ERROR: Input p.d.f. must be extendible", GetName()));
120 }
121 }
122
123 _extended = isPdf;
124 _yvar = yvar ? static_cast<RooRealVar *>(_dataClone->get()->find(yvar->GetName())) : nullptr;
125
126 initialize();
127}
128/// \endcond ROOFIT_INTERNAL
129
130
131////////////////////////////////////////////////////////////////////////////////
132/// Copy constructor
133
136 _extended(other._extended),
137 _integrate(other._integrate),
138 _yvar(other._yvar ? static_cast<RooRealVar *>(_dataClone->get()->find(other._yvar->GetName())) : nullptr),
139 _intConfig(other._intConfig)
140{
141
142 initialize() ;
143
144}
145
146
147////////////////////////////////////////////////////////////////////////////////
148/// Common constructor initialization
149
151{
152 // If this tests statistic is not a "Slave" in the RooAbsOptTestStatistic
153 // framework, it doesn't do any actual computation and no initialization is
154 // needed. It would not even work, because _funcObsSet would be a nullptr.
155 if(operMode() != Slave) return;
156
157 for(RooAbsArg * arg : *_funcObsSet) {
158 if (auto* var = dynamic_cast<RooRealVar*>(arg)) {
160 }
161 }
162 if (_yvar) {
164 }
165
166 // Define alternate numeric integrator configuration for bin integration
167 // We expect bin contents to very only very slowly so a non-adaptive
168 // Gauss-Kronrod integrator is expected to perform well (if RooFitMore is available)
171#ifdef R__HAS_MATHMORE
172 _intConfig.method1D().setLabel("RooGaussKronrodIntegrator1D") ;
173#endif
175
177
178}
179
180
181
182////////////////////////////////////////////////////////////////////////////////
183/// Initialize bin content integrator
184
186{
187 if (!_funcInt) {
188 _funcInt = std::unique_ptr<RooAbsReal>{_funcClone->createIntegral(_rrvArgs,_rrvArgs,_intConfig,"bin")};
189 for(auto * x : static_range_cast<RooRealVar*>(_rrvArgs)) {
190 _binList.push_back(&x->getBinning("bin",false,true)) ;
191 }
192 }
193
194}
195
196
198
199
200////////////////////////////////////////////////////////////////////////////////
201/// Calculate contribution to internal error due to error on 'x' coordinates
202/// at point i
203
204double RooXYChi2Var::xErrorContribution(double ydata) const
205{
206 double ret(0) ;
207
208 for(auto * var : static_range_cast<RooRealVar*>(_rrvArgs)) {
209
210 if (var->hasAsymError()) {
211
212 // Get value at central X
213 double cxval = var->getVal() ;
214 double xerrLo = -var->getAsymErrorLo() ;
215 double xerrHi = var->getAsymErrorHi() ;
216 double xerr = (xerrLo+xerrHi)/2 ;
217
218 // Get value at X-eps
219 var->setVal(cxval - xerr/100) ;
220 double fxmin = fy() ;
221
222 // Get value at X+eps
223 var->setVal(cxval + xerr/100) ;
224 double fxmax = fy() ;
225
226 // Calculate slope
227 double slope = (fxmax-fxmin)/(2*xerr/100.) ;
228
229// cout << "xerrHi = " << xerrHi << " xerrLo = " << xerrLo << " slope = " << slope << std::endl ;
230
231 // Asymmetric X error, decide which one to use
232 if ((ydata>cxval && fxmax>fxmin) || (ydata<=cxval && fxmax<=fxmin)) {
233 // Use right X error
234 ret += pow(xerrHi*slope,2) ;
235 } else {
236 // Use left X error
237 ret += pow(xerrLo*slope,2) ;
238 }
239
240 } else if (var->hasError()) {
241
242 // Get value at central X
243 double cxval = var->getVal() ;
244 double xerr = var->getError() ;
245
246 // Get value at X-eps
247 var->setVal(cxval - xerr/100) ;
248 double fxmin = fy() ;
249
250 // Get value at X+eps
251 var->setVal(cxval + xerr/100) ;
252 double fxmax = fy() ;
253
254 // Calculate slope
255 double slope = (fxmax-fxmin)/(2*xerr/100.) ;
256
257// cout << var << " " ;
258// var->Print() ;
259
260// cout << var->GetName() << " xerr = " << xerr << " slope = " << slope << std::endl ;
261
262 // Symmetric X error
263 ret += pow(xerr*slope,2) ;
264 }
265 }
266 return ret ;
267}
268
269
270
271
272////////////////////////////////////////////////////////////////////////////////
273/// Return function value requested bu present configuration
274///
275/// If integration is required, the function value integrated
276/// over the bin volume divided by the bin volume is returned,
277/// otherwise the value at the bin center is returned.
278/// The bin volume is defined by the error on the 'X' coordinates
279///
280/// If an extended p.d.f. is used as function, its value is
281/// also multiplied by the expected number of events here
282
283double RooXYChi2Var::fy() const
284{
285 // Get function value
286 double yfunc ;
287 if (!_integrate) {
288 yfunc = _funcClone->getVal(_dataClone->get()) ;
289 } else {
290 double volume(1) ;
291 auto rrvIter = _rrvArgs.begin();
292 for (auto iter = _binList.begin() ; iter != _binList.end() ; ++iter) {
293 auto* x = static_cast<RooRealVar*>(*rrvIter);
294 double xmin = x->getVal() + x->getErrorLo() ;
295 double xmax = x->getVal() + x->getErrorHi() ;
296 (*iter)->setRange(xmin,xmax) ;
297 x->setShapeDirty() ;
298 volume *= (xmax - xmin) ;
299 ++rrvIter;
300 }
301 double ret = _funcInt->getVal() ;
302 return ret / volume ;
303 }
304 if (_extended) {
305 RooAbsPdf* pdf = static_cast<RooAbsPdf*>(_funcClone) ;
306 // Multiply with expected number of events
307 yfunc *= pdf->expectedEvents(_dataClone->get()) ;
308 }
309 return yfunc ;
310}
311
312
313
314////////////////////////////////////////////////////////////////////////////////
315/// Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize
316
317double RooXYChi2Var::evaluatePartition(std::size_t firstEvent, std::size_t lastEvent, std::size_t stepSize) const
318{
319 double result(0);
320 double carry(0);
321
322 // Loop over bins of dataset
323 RooDataSet* xydata = static_cast<RooDataSet*>(_dataClone) ;
324
325 for (auto i=firstEvent ; i<lastEvent ; i+=stepSize) {
326
327 // get the data values for this event
328 xydata->get(i);
329
330 // Get function value
331 double yfunc = fy() ;
332
333 // Get data value and error
334 double ydata ;
335 double eylo;
336 double eyhi;
337 if (_yvar) {
338 ydata = _yvar->getVal() ;
339 eylo = -1*_yvar->getErrorLo() ;
340 eyhi = _yvar->getErrorHi() ;
341 } else {
342 ydata = xydata->weight() ;
343 xydata->weightError(eylo,eyhi) ;
344 }
345
346 // Calculate external error
347 double eExt = yfunc-ydata ;
348
349 // Pick upper or lower error bar depending on sign of external error
350 double eInt = (eExt>0) ? eyhi : eylo ;
351
352 // Add contributions due to error in x coordinates
353 double eIntX2 = _integrate ? 0 : xErrorContribution(ydata) ;
354
355// cout << "fy = " << yfunc << " eExt = " << eExt << " eInt = " << eInt << " eIntX2 = " << eIntX2 << std::endl ;
356
357 // Return 0 if eInt=0, special handling in MINUIT will follow
358 if (eInt==0.) {
359 coutE(Eval) << "RooXYChi2Var::RooXYChi2Var(" << GetName() << ") INFINITY ERROR: data point " << i
360 << " has zero error, but function is not zero (f=" << yfunc << ")" << std::endl ;
361 return 0 ;
362 }
363
365 double term = eExt*eExt/(eInt*eInt+ eIntX2);
366 double y = term - carry;
367 double t = result + y;
368 carry = (t - result) - y;
369 result = t;
370 }
371
372 _evalCarry = carry;
373 return result ;
374}
375
376
377
379{
380 // Inform base class that observable yvar cannot be optimized away from the dataset
381 if (_yvar) return RooArgSet(*_yvar) ;
382 return RooArgSet() ;
383}
#define e(i)
Definition RSha256.hxx:103
#define coutE(a)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
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
float xmax
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:79
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
const_iterator begin() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
Abstract base class for test statistics objects that evaluate a function or PDF at each point of a gi...
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.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
double _evalCarry
! carry of Kahan sum in evaluatePartition
GOFOpMode operMode() const
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
bool setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
Container class to hold unbinned data.
Definition RooDataSet.h:33
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.
double weight() const override
Return event weight of current event.
RooCategory & methodND()
void setEpsRel(double newEpsRel)
Set relative convergence criteria (convergence if std::abs(Err)/abs(Int)<newEpsRel)
RooCategory & method1D()
void setEpsAbs(double newEpsAbs)
Set absolute convergence criteria (convergence if std::abs(Err)<newEpsAbs)
Variable that can be changed from the outside.
Definition RooRealVar.h:37
double getErrorLo() const
Definition RooRealVar.h:67
double getErrorHi() const
Definition RooRealVar.h:68
RooXYChi2Var implements a simple chi^2 calculation from an unbinned dataset with values x,...
double fy() const
Return function value requested bu present configuration.
std::unique_ptr< RooAbsReal > _funcInt
! Function integral
std::list< RooAbsBinning * > _binList
! Bin ranges
bool _extended
Is the input function and extended p.d.f.
void initIntegrator()
Initialize bin content integrator.
RooRealVar * _yvar
Y variable if so designated.
RooXYChi2Var(const char *name, const char *title, RooAbsReal &func, RooDataSet &data, bool integrate=false)
RooXYChi2Var constructor with function and X-Y values dataset.
RooArgSet requiredExtraObservables() const override
double 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.
bool _integrate
Is integration over the bin volume requested.
RooArgSet _rrvArgs
Set of real-valued observables.
double xErrorContribution(double ydata) const
Calculate contribution to internal error due to error on 'x' coordinates at point i.
void initialize()
Common constructor initialization.
RooNumIntConfig _intConfig
Numeric integrator configuration for integration of function over bin.
~RooXYChi2Var() override
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17