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 "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///
63/// RooXYChi2Var constructor with function and X-Y values dataset
64///
65/// If the function is a pdf, it must be extendable. in this case, hhe value of
66/// the function that defines the chi^2 in this form is takes as the p.d.f.
67/// times the expected number of events
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 values that correspond to the Y and its error
78///
79
80RooXYChi2Var::RooXYChi2Var(const char *name, const char *title, RooAbsReal &func, RooDataSet &xydata, bool integrate)
81 : RooXYChi2Var{name, title, func, xydata, nullptr, integrate, makeRooAbsTestStatisticCfg()}
82{
83}
84
85////////////////////////////////////////////////////////////////////////////////
86///
87/// RooXYChi2Var constructor with function and X-Y values dataset.
88///
89/// If the function is a pdf, it must be extendable. in this case, hhe value of
90/// the function that defines the chi^2 in this form is takes as the p.d.f.
91/// times the expected number of events
92///
93/// An X-Y dataset is a weighted dataset with one or more observables X where given yvar is interpreted
94/// 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.
95///
96/// To store errors associated with the x and y values in a RooDataSet, call RooRealVar::setAttribute("StoreError")
97/// on each X-type observable for which the error should be stored and add datapoints to the dataset as follows
98///
99/// RooDataSet::add(xset,yval,yerr) where xset is the RooArgSet of x observables (with or without errors) and yval and yerr
100/// are the double values that correspond to the Y and its error
101///
102
103RooXYChi2Var::RooXYChi2Var(const char *name, const char *title, RooAbsReal &func, RooDataSet &xydata, RooRealVar &yvar,
104 bool integrate)
105 : RooXYChi2Var{name, title, func, xydata, &yvar, integrate, makeRooAbsTestStatisticCfg()}
106{
107}
108
109
110/// \cond ROOFIT_INTERNAL
111// For internal use in RooAbsReal::createChi2().
112RooXYChi2Var::RooXYChi2Var(const char *name, const char *title, RooAbsReal &func, RooAbsData &data, RooRealVar *yvar,
113 bool integrate, RooAbsTestStatistic::Configuration const &cfg)
114 : RooAbsOptTestStatistic(name, title, func, data, RooArgSet(), cfg),
115 _integrate(integrate),
116 _intConfig(*defaultIntegratorConfig())
117{
118 bool isPdf = dynamic_cast<RooAbsPdf const *>(&func) != nullptr;
119
120 if (isPdf) {
121 auto &extPdf = static_cast<RooAbsPdf const &>(func);
122 if (!extPdf.canBeExtended()) {
123 throw std::runtime_error(
124 Form("RooXYChi2Var::RooXYChi2Var(%s) ERROR: Input p.d.f. must be extendible", GetName()));
125 }
126 }
127
128 _extended = isPdf;
129 _yvar = yvar ? static_cast<RooRealVar *>(_dataClone->get()->find(yvar->GetName())) : nullptr;
130
131 initialize();
132}
133/// \endcond ROOFIT_INTERNAL
134
135
136////////////////////////////////////////////////////////////////////////////////
137/// Copy constructor
138
141 _extended(other._extended),
142 _integrate(other._integrate),
143 _yvar(other._yvar ? static_cast<RooRealVar *>(_dataClone->get()->find(other._yvar->GetName())) : nullptr),
145{
146
147 initialize() ;
148
149}
150
151
152////////////////////////////////////////////////////////////////////////////////
153/// Common constructor initialization
154
156{
157 // If this tests statistic is not a "Slave" in the RooAbsOptTestStatistic
158 // framework, it doesn't do any actual computation and no initialization is
159 // needed. It would not even work, because _funcObsSet would be a nullptr.
160 if(operMode() != Slave) return;
161
162 for(RooAbsArg * arg : *_funcObsSet) {
163 if (auto* var = dynamic_cast<RooRealVar*>(arg)) {
164 _rrvArgs.add(*var) ;
165 }
166 }
167 if (_yvar) {
168 _rrvArgs.add(*_yvar) ;
169 }
170
171 // Define alternate numeric integrator configuration for bin integration
172 // We expect bin contents to very only very slowly so a non-adaptive
173 // Gauss-Kronrod integrator is expected to perform well (if RooFitMore is available)
176#ifdef R__HAS_MATHMORE
177 _intConfig.method1D().setLabel("RooGaussKronrodIntegrator1D") ;
178#endif
179 _intConfig.methodND().setLabel("RooAdaptiveIntegratorND") ;
180
182
183}
184
185
186
187////////////////////////////////////////////////////////////////////////////////
188/// Initialize bin content integrator
189
191{
192 if (!_funcInt) {
193 _funcInt = std::unique_ptr<RooAbsReal>{_funcClone->createIntegral(_rrvArgs,_rrvArgs,_intConfig,"bin")};
194 for(auto * x : static_range_cast<RooRealVar*>(_rrvArgs)) {
195 _binList.push_back(&x->getBinning("bin",false,true)) ;
196 }
197 }
198
199}
200
201
203
204
205////////////////////////////////////////////////////////////////////////////////
206/// Calculate contribution to internal error due to error on 'x' coordinates
207/// at point i
208
209double RooXYChi2Var::xErrorContribution(double ydata) const
210{
211 double ret(0) ;
212
213 for(auto * var : static_range_cast<RooRealVar*>(_rrvArgs)) {
214
215 if (var->hasAsymError()) {
216
217 // Get value at central X
218 double cxval = var->getVal() ;
219 double xerrLo = -var->getAsymErrorLo() ;
220 double xerrHi = var->getAsymErrorHi() ;
221 double xerr = (xerrLo+xerrHi)/2 ;
222
223 // Get value at X-eps
224 var->setVal(cxval - xerr/100) ;
225 double fxmin = fy() ;
226
227 // Get value at X+eps
228 var->setVal(cxval + xerr/100) ;
229 double fxmax = fy() ;
230
231 // Calculate slope
232 double slope = (fxmax-fxmin)/(2*xerr/100.) ;
233
234// cout << "xerrHi = " << xerrHi << " xerrLo = " << xerrLo << " slope = " << slope << endl ;
235
236 // Asymmetric X error, decide which one to use
237 if ((ydata>cxval && fxmax>fxmin) || (ydata<=cxval && fxmax<=fxmin)) {
238 // Use right X error
239 ret += pow(xerrHi*slope,2) ;
240 } else {
241 // Use left X error
242 ret += pow(xerrLo*slope,2) ;
243 }
244
245 } else if (var->hasError()) {
246
247 // Get value at central X
248 double cxval = var->getVal() ;
249 double xerr = var->getError() ;
250
251 // Get value at X-eps
252 var->setVal(cxval - xerr/100) ;
253 double fxmin = fy() ;
254
255 // Get value at X+eps
256 var->setVal(cxval + xerr/100) ;
257 double fxmax = fy() ;
258
259 // Calculate slope
260 double slope = (fxmax-fxmin)/(2*xerr/100.) ;
261
262// cout << var << " " ;
263// var->Print() ;
264
265// cout << var->GetName() << " xerr = " << xerr << " slope = " << slope << endl ;
266
267 // Symmetric X error
268 ret += pow(xerr*slope,2) ;
269 }
270 }
271 return ret ;
272}
273
274
275
276
277////////////////////////////////////////////////////////////////////////////////
278/// Return function value requested bu present configuration
279///
280/// If integration is required, the function value integrated
281/// over the bin volume divided by the bin volume is returned,
282/// otherwise the value at the bin center is returned.
283/// The bin volume is defined by the error on the 'X' coordinates
284///
285/// If an extended p.d.f. is used as function, its value is
286/// also multiplied by the expected number of events here
287
288double RooXYChi2Var::fy() const
289{
290 // Get function value
291 double yfunc ;
292 if (!_integrate) {
293 yfunc = _funcClone->getVal(_dataClone->get()) ;
294 } else {
295 double volume(1) ;
296 auto rrvIter = _rrvArgs.begin();
297 for (auto iter = _binList.begin() ; iter != _binList.end() ; ++iter) {
298 auto* x = static_cast<RooRealVar*>(*rrvIter);
299 double xmin = x->getVal() + x->getErrorLo() ;
300 double xmax = x->getVal() + x->getErrorHi() ;
301 (*iter)->setRange(xmin,xmax) ;
302 x->setShapeDirty() ;
303 volume *= (xmax - xmin) ;
304 ++rrvIter;
305 }
306 double ret = _funcInt->getVal() ;
307 return ret / volume ;
308 }
309 if (_extended) {
310 RooAbsPdf* pdf = static_cast<RooAbsPdf*>(_funcClone) ;
311 // Multiply with expected number of events
312 yfunc *= pdf->expectedEvents(_dataClone->get()) ;
313 }
314 return yfunc ;
315}
316
317
318
319////////////////////////////////////////////////////////////////////////////////
320/// Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize
321
322double RooXYChi2Var::evaluatePartition(std::size_t firstEvent, std::size_t lastEvent, std::size_t stepSize) const
323{
324 double result(0);
325 double carry(0);
326
327 // Loop over bins of dataset
328 RooDataSet* xydata = static_cast<RooDataSet*>(_dataClone) ;
329
330 for (auto i=firstEvent ; i<lastEvent ; i+=stepSize) {
331
332 // get the data values for this event
333 xydata->get(i);
334
335 // Get function value
336 double yfunc = fy() ;
337
338 // Get data value and error
339 double ydata ;
340 double eylo;
341 double eyhi;
342 if (_yvar) {
343 ydata = _yvar->getVal() ;
344 eylo = -1*_yvar->getErrorLo() ;
345 eyhi = _yvar->getErrorHi() ;
346 } else {
347 ydata = xydata->weight() ;
348 xydata->weightError(eylo,eyhi) ;
349 }
350
351 // Calculate external error
352 double eExt = yfunc-ydata ;
353
354 // Pick upper or lower error bar depending on sign of external error
355 double eInt = (eExt>0) ? eyhi : eylo ;
356
357 // Add contributions due to error in x coordinates
358 double eIntX2 = _integrate ? 0 : xErrorContribution(ydata) ;
359
360// cout << "fy = " << yfunc << " eExt = " << eExt << " eInt = " << eInt << " eIntX2 = " << eIntX2 << endl ;
361
362 // Return 0 if eInt=0, special handling in MINUIT will follow
363 if (eInt==0.) {
364 coutE(Eval) << "RooXYChi2Var::RooXYChi2Var(" << GetName() << ") INFINITY ERROR: data point " << i
365 << " has zero error, but function is not zero (f=" << yfunc << ")" << endl ;
366 return 0 ;
367 }
368
369 // Add chi2 term
370 double term = eExt*eExt/(eInt*eInt+ eIntX2);
371 double y = term - carry;
372 double t = result + y;
373 carry = (t - result) - y;
374 result = t;
375 }
376
377 _evalCarry = carry;
378 return result ;
379}
380
381
382
384{
385 // Inform base class that observable yvar cannot be optimized away from the dataset
386 if (_yvar) return RooArgSet(*_yvar) ;
387 return RooArgSet() ;
388}
#define e(i)
Definition RSha256.hxx:103
RooAbsData * _dataClone
Pointer to internal clone if input data.
#define coutE(a)
bool _extended
Definition RooNLLVar.h:43
RooRealVar * _yvar
Y variable if so designated.
RooNumIntConfig _intConfig
Numeric integrator configuration for integration of function over bin.
bool _integrate
Is integration over the bin volume requested.
#define ClassImp(name)
Definition Rtypes.h:377
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:77
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
TIterator begin()
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:55
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:57
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