Logo ROOT  
Reference Guide
RooChi2Var.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 RooChi2Var
19 \ingroup Roofitcore
20 \brief RooChi2Var implements a simple \f$ \chi^2 \f$ calculation from a binned dataset and a PDF.
21 *
22 * It calculates:
23 *
24 \f{align*}{
25 \chi^2 &= \sum_{\mathrm{bins}} \left( \frac{N_\mathrm{PDF,bin} - N_\mathrm{Data,bin}}{\Delta_\mathrm{bin}} \right)^2 \\
26 N_\mathrm{PDF,bin} &=
27 \begin{cases}
28 \mathrm{pdf}(\text{bin centre}) \cdot V_\mathrm{bin} \cdot N_\mathrm{Data,tot} &\text{normal PDF}\\
29 \mathrm{pdf}(\text{bin centre}) \cdot V_\mathrm{bin} \cdot N_\mathrm{Data,expected} &\text{extended PDF}
30 \end{cases} \\
31 \Delta_\mathrm{bin} &=
32 \begin{cases}
33 \sqrt{N_\mathrm{PDF,bin}} &\text{if } \mathtt{DataError == RooAbsData::Expected}\\
34 \mathtt{data{\rightarrow}weightError()} &\text{otherwise} \\
35 \end{cases}
36 \f}
37 * If the dataset doesn't have user-defined errors, errors are assumed to be \f$ \sqrt{N} \f$.
38 * In extended PDF mode, N_tot (total number of data events) is substituted with N_expected, the
39 * expected number of events that the PDF predicts.
40 *
41 * \note If the dataset has errors stored, empty bins will prevent the calculation of \f$ \chi^2 \f$, because those have
42 * zero error. This leads to messages like:
43 * ```
44 * [#0] ERROR:Eval -- RooChi2Var::RooChi2Var(chi2_GenPdf_data_hist) INFINITY ERROR: bin 2 has zero error
45 * ```
46 *
47 * \note In this case, one can use the expected errors of the PDF instead of the data errors:
48 * ```{.cpp}
49 * RooChi2Var chi2(..., ..., RooFit::DataError(RooAbsData::Expected), ...);
50 * ```
51 */
52
53#include "RooFit.h"
54
55#include "RooChi2Var.h"
56#include "RooDataHist.h"
57#include "RooAbsPdf.h"
58#include "RooCmdConfig.h"
59#include "RooMsgService.h"
60
61#include "Riostream.h"
62#include "TClass.h"
63
64#include "RooRealVar.h"
65#include "RooAbsDataStore.h"
66
67
68using namespace std;
69
70namespace {
71 template<class ...Args>
72 RooAbsTestStatistic::Configuration makeRooAbsTestStatisticCfgForFunc(Args const& ... args) {
74 cfg.rangeName = RooCmdConfig::decodeStringOnTheFly("RooChi2Var::RooChi2Var","RangeWithName",0,"",args...);
75 cfg.nCPU = RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","NumCPU",0,1,args...);
77 cfg.verbose = static_cast<bool>(RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","Verbose",0,1,args...));
78 cfg.cloneInputData = false;
79 cfg.integrateOverBinsPrecision = RooCmdConfig::decodeDoubleOnTheFly("RooChi2Var::RooChi2Var", "IntegrateBins", 0, -1., {args...});
80 return cfg;
81 }
82
83 template<class ...Args>
84 RooAbsTestStatistic::Configuration makeRooAbsTestStatisticCfgForPdf(Args const& ... args) {
85 auto cfg = makeRooAbsTestStatisticCfgForFunc(args...);
86 cfg.addCoefRangeName = RooCmdConfig::decodeStringOnTheFly("RooChi2Var::RooChi2Var","AddCoefRange",0,"",args...);
87 cfg.splitCutRange = static_cast<bool>(RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","SplitRange",0,0,args...));
88 return cfg;
89 }
90}
91
93;
94
96
97
98////////////////////////////////////////////////////////////////////////////////
99/// RooChi2Var constructor. Optional arguments are:
100/// \param[in] name Name of the PDF
101/// \param[in] title Title for plotting etc.
102/// \param[in] func Function
103/// \param[in] hdata Data histogram
104/// \param[in] argX Optional arguments according to table below.
105/// <table>
106/// <tr><th> Argument <th> Effect
107/// <tr><td>
108/// DataError() <td> Choose between Poisson errors and Sum-of-weights errors
109/// <tr><td>
110/// NumCPU() <td> Activate parallel processing feature
111/// <tr><td>
112/// Range() <td> Fit only selected region
113/// <tr><td>
114/// Verbose() <td> Verbose output of GOF framework
115/// <tr><td>
116/// IntegrateBins() <td> Integrate PDF within each bin. This sets the desired precision. Only useful for binned fits.
117
118RooChi2Var::RooChi2Var(const char *name, const char* title, RooAbsReal& func, RooDataHist& hdata,
119 const RooCmdArg& arg1,const RooCmdArg& arg2,const RooCmdArg& arg3,
120 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,
121 const RooCmdArg& arg7,const RooCmdArg& arg8,const RooCmdArg& arg9) :
122 RooAbsOptTestStatistic(name,title,func,hdata,_emptySet,
123 makeRooAbsTestStatisticCfgForFunc(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9))
124{
125 RooCmdConfig pc("RooChi2Var::RooChi2Var") ;
126 pc.defineInt("etype","DataError",0,(Int_t)RooDataHist::Auto) ;
127 pc.defineInt("extended","Extended",0,kFALSE) ;
128 pc.allowUndefined() ;
129
130 pc.process(arg1) ; pc.process(arg2) ; pc.process(arg3) ;
131 pc.process(arg4) ; pc.process(arg5) ; pc.process(arg6) ;
132 pc.process(arg7) ; pc.process(arg8) ; pc.process(arg9) ;
133
134 if (func.IsA()->InheritsFrom(RooAbsPdf::Class())) {
135 _funcMode = pc.getInt("extended") ? ExtendedPdf : Pdf ;
136 } else {
138 }
139 _etype = (RooDataHist::ErrorType) pc.getInt("etype") ;
140
143 }
144
145}
146
147
148////////////////////////////////////////////////////////////////////////////////
149/// RooChi2Var constructor. Optional arguments taken
150///
151/// \param[in] name Name of the PDF
152/// \param[in] title Title for plotting etc.
153/// \param[in] pdf PDF to fit
154/// \param[in] hdata Data histogram
155/// \param[in] argX Optional arguments according to table below.
156/// <table>
157/// <tr><th> Argument <th> Effect
158/// <tr><td>
159/// Extended() <td> Include extended term in calculation
160/// <tr><td>
161/// DataError() <td> Choose between Poisson errors and Sum-of-weights errors
162/// <tr><td>
163/// NumCPU() <td> Activate parallel processing feature
164/// <tr><td>
165/// Range() <td> Fit only selected region
166/// <tr><td>
167/// SumCoefRange() <td> Set the range in which to interpret the coefficients of RooAddPdf components
168/// <tr><td>
169/// SplitRange() <td> Fit range is split by index category of simultaneous PDF
170/// <tr><td>
171/// ConditionalObservables() <td> Define projected observables
172/// <tr><td>
173/// Verbose() <td> Verbose output of GOF framework
174/// <tr><td>
175/// IntegrateBins() <td> Integrate PDF within each bin. This sets the desired precision.
176
177RooChi2Var::RooChi2Var(const char *name, const char* title, RooAbsPdf& pdf, RooDataHist& hdata,
178 const RooCmdArg& arg1,const RooCmdArg& arg2,const RooCmdArg& arg3,
179 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,
180 const RooCmdArg& arg7,const RooCmdArg& arg8,const RooCmdArg& arg9) :
181 RooAbsOptTestStatistic(name,title,pdf,hdata,
182 *static_cast<const RooArgSet*>(RooCmdConfig::decodeObjOnTheFly("RooChi2Var::RooChi2Var","ProjectedObservables",0,&_emptySet,
183 arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9)),
184 makeRooAbsTestStatisticCfgForPdf(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9))
185{
186 RooCmdConfig pc("RooChi2Var::RooChi2Var") ;
187 pc.defineInt("extended","Extended",0,kFALSE) ;
188 pc.defineInt("etype","DataError",0,(Int_t)RooDataHist::Auto) ;
189 pc.allowUndefined() ;
190
191 pc.process(arg1) ; pc.process(arg2) ; pc.process(arg3) ;
192 pc.process(arg4) ; pc.process(arg5) ; pc.process(arg6) ;
193 pc.process(arg7) ; pc.process(arg8) ; pc.process(arg9) ;
194
195 _funcMode = pc.getInt("extended") ? ExtendedPdf : Pdf ;
196 _etype = (RooDataHist::ErrorType) pc.getInt("etype") ;
199 }
200}
201
202
203////////////////////////////////////////////////////////////////////////////////
204/// Constructor of a chi2 for given p.d.f. with respect given binned
205/// dataset. If cutRange is specified the calculation of the chi2 is
206/// restricted to that named range. If addCoefRange is specified, the
207/// interpretation of fractions for all component RooAddPdfs that do
208/// not have a frozen range interpretation is set to chosen range
209/// name. If nCPU is greater than one the chi^2 calculation is
210/// parallelized over the specified number of processors. If
211/// interleave is true the partitioning of event over processors
212/// follows a (i % n == i_set) strategy rather than a bulk
213/// partitioning strategy which may result in unequal load balancing
214/// in binned datasets with many (adjacent) zero bins. If
215/// splitCutRange is true the cutRange is used to construct an
216/// individual cutRange for each RooSimultaneous index category state
217/// name cutRange_{indexStateName}.
218
219RooChi2Var::RooChi2Var(const char *name, const char *title, RooAbsPdf& pdf, RooDataHist& hdata,
220 RooAbsTestStatistic::Configuration const& cfg, bool extended, RooDataHist::ErrorType etype) :
221 RooAbsOptTestStatistic(name,title,pdf,hdata,RooArgSet(), cfg),
222 _etype(etype), _funcMode(extended?ExtendedPdf:Pdf)
223{
224}
225
226
227////////////////////////////////////////////////////////////////////////////////
228/// Constructor of a chi2 for given p.d.f. with respect given binned
229/// dataset taking the observables specified in projDeps as projected
230/// observables. If cutRange is specified the calculation of the chi2
231/// is restricted to that named range. If addCoefRange is specified,
232/// the interpretation of fractions for all component RooAddPdfs that
233/// do not have a frozen range interpretation is set to chosen range
234/// name. If nCPU is greater than one the chi^2 calculation is
235/// parallelized over the specified number of processors. If
236/// interleave is true the partitioning of event over processors
237/// follows a (i % n == i_set) strategy rather than a bulk
238/// partitioning strategy which may result in unequal load balancing
239/// in binned datasets with many (adjacent) zero bins. If
240/// splitCutRange is true the cutRange is used to construct an
241/// individual cutRange for each RooSimultaneous index category state
242/// name cutRange_{indexStateName}.
243
244RooChi2Var::RooChi2Var(const char *name, const char *title, RooAbsReal& func, RooDataHist& hdata,
245 const RooArgSet& projDeps, RooChi2Var::FuncMode fmode,
248 RooAbsOptTestStatistic(name,title,func,hdata,projDeps,cfg),
249 _etype(etype), _funcMode(fmode)
250{
251}
252
253
254////////////////////////////////////////////////////////////////////////////////
255/// Copy constructor
256
257RooChi2Var::RooChi2Var(const RooChi2Var& other, const char* name) :
259 _etype(other._etype),
260 _funcMode(other._funcMode)
261{
262}
263
264
265////////////////////////////////////////////////////////////////////////////////
266/// Destructor
267
269{
270}
271
272
273////////////////////////////////////////////////////////////////////////////////
274/// Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize
275/// Throughout the calculation, we use Kahan's algorithm for summing to
276/// prevent loss of precision - this is a factor four more expensive than
277/// straight addition, but since evaluating the PDF is usually much more
278/// expensive than that, we tolerate the additional cost...
279
280Double_t RooChi2Var::evaluatePartition(std::size_t firstEvent, std::size_t lastEvent, std::size_t stepSize) const
281{
282
283 Double_t result(0), carry(0);
284
285 _dataClone->store()->recalculateCache( _projDeps, firstEvent, lastEvent, stepSize, kFALSE) ;
286
287
288 // Determine normalization factor depending on type of input function
289 Double_t normFactor(1) ;
290 switch (_funcMode) {
291 case Function: normFactor=1 ; break ;
292 case Pdf: normFactor = _dataClone->sumEntries() ; break ;
293 case ExtendedPdf: normFactor = ((RooAbsPdf*)_funcClone)->expectedEvents(_dataClone->get()) ; break ;
294 }
295
296 // Loop over bins of dataset
298 for (auto i=firstEvent ; i<lastEvent ; i+=stepSize) {
299
300 // get the data values for this event
301 hdata->get(i);
302
303 if (!hdata->valid()) continue;
304
305 const Double_t nData = hdata->weight() ;
306
307 const Double_t nPdf = _funcClone->getVal(_normSet) * normFactor * hdata->binVolume() ;
308
309 const Double_t eExt = nPdf-nData ;
310
311
312 Double_t eInt ;
314 Double_t eIntLo,eIntHi ;
315 hdata->weightError(eIntLo,eIntHi,_etype) ;
316 eInt = (eExt>0) ? eIntHi : eIntLo ;
317 } else {
318 eInt = sqrt(nPdf) ;
319 }
320
321 // Skip cases where pdf=0 and there is no data
322 if (0. == eInt * eInt && 0. == nData * nData && 0. == nPdf * nPdf) continue ;
323
324 // Return 0 if eInt=0, special handling in MINUIT will follow
325 if (0. == eInt * eInt) {
326 coutE(Eval) << "RooChi2Var::RooChi2Var(" << GetName() << ") INFINITY ERROR: bin " << i
327 << " has zero error" << endl ;
328 return 0.;
329 }
330
331// cout << "Chi2Var[" << i << "] nData = " << nData << " nPdf = " << nPdf << " errorExt = " << eExt << " errorInt = " << eInt << " contrib = " << eExt*eExt/(eInt*eInt) << endl ;
332
333 Double_t term = eExt*eExt/(eInt*eInt) ;
334 Double_t y = term - carry;
335 Double_t t = result + y;
336 carry = (t - result) - y;
337 result = t;
338 }
339
340 _evalCarry = carry;
341 return result ;
342}
void Class()
Definition: Class.C:29
#define coutE(a)
Definition: RooMsgService.h:33
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:364
char name[80]
Definition: TGX11.cxx:110
double sqrt(double)
virtual void recalculateCache(const RooArgSet *, Int_t, Int_t, Int_t, Bool_t)
virtual const RooArgSet * get() const
Definition: RooAbsData.h:125
RooAbsDataStore * store()
Definition: RooAbsData.h:101
virtual Double_t sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
RooAbsOptTestStatistic is the abstract base class for test statistics objects that evaluate a functio...
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
Double_t _evalCarry
Offset as KahanSum to avoid loss of precision.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:35
RooChi2Var implements a simple calculation from a binned dataset and a PDF.
Definition: RooChi2Var.h:25
RooDataHist::ErrorType _etype
Definition: RooChi2Var.h:70
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 Throughout the calcula...
Definition: RooChi2Var.cxx:280
FuncMode _funcMode
Definition: RooChi2Var.h:71
RooChi2Var(const char *name, const char *title, RooAbsReal &func, RooDataHist &data, 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 RooCmdArg &arg9=RooCmdArg::none())
RooChi2Var constructor.
Definition: RooChi2Var.cxx:118
static RooArgSet _emptySet
Definition: RooChi2Var.h:68
virtual ~RooChi2Var()
Destructor.
Definition: RooChi2Var.cxx:268
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:27
Class RooCmdConfig is a configurable parser for RooCmdArg named arguments.
Definition: RooCmdConfig.h:27
static std::string decodeStringOnTheFly(const char *callerID, const char *cmdArgName, Int_t intIdx, const char *defVal, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Static decoder function allows to retrieve string property from set of RooCmdArgs For use in base mem...
static Int_t decodeIntOnTheFly(const char *callerID, const char *cmdArgName, Int_t intIdx, Int_t defVal, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Static decoder function allows to retrieve integer property from set of RooCmdArgs For use in base me...
static double decodeDoubleOnTheFly(const char *callerID, const char *cmdArgName, int idx, double defVal, std::initializer_list< std::reference_wrapper< const RooCmdArg > > args)
Find a given double in a list of RooCmdArg.
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:44
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:111
void weightError(double &lo, double &hi, ErrorType etype=Poisson) const override
Return the asymmetric errors on the current weight.
Bool_t isNonPoissonWeighted() const override
Returns true if dataset contains entries with a non-integer weight.
bool valid(std::size_t i) const
Return true if bin i is considered valid within the current range definitions of all observables.
Definition: RooDataHist.h:120
double binVolume(std::size_t i) const
Return bin volume of i-th bin.
Definition: RooDataHist.h:117
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:83
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:445
Double_t y[n]
Definition: legend1.C:17
const UInt_t eInt[256]
@ Interleave
Definition: RooGlobalFunc.h:63
static constexpr double pc
std::string rangeName
Stores the configuration parameters for RooAbsTestStatistic.