Logo ROOT   6.18/05
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/**
19// \class RooChi2Var
20// Class RooChi2Var implements a simple chi^2 calculation from a binned dataset
21// and a PDF. It calculates
22\f[
23 \chi^2 = \sum_{[\mathrm{bins}]} \left( \frac{(f_\mathrm{PDF} \cdot N_\mathrm{tot} / V_\mathrm{bin}) - N_\mathrm{bin}}{\mathrm{err}_\mathrm{bin}} \right)^2
24\f]
25// If no user-defined errors are defined for the dataset, poisson errors
26// are used. In extended PDF mode, N_tot is substituted with N_expected.
27//
28*/
29
30#include "RooFit.h"
31
32#include "RooChi2Var.h"
33#include "RooChi2Var.h"
34#include "RooDataHist.h"
35#include "RooAbsPdf.h"
36#include "RooCmdConfig.h"
37#include "RooMsgService.h"
38
39#include "Riostream.h"
40
41#include "RooRealVar.h"
42#include "RooAbsDataStore.h"
43
44
45using namespace std;
46
48;
49
51
52
53////////////////////////////////////////////////////////////////////////////////
54/// RooChi2Var constructor. Optional arguments are:
55/// \param[in] name Name of the PDF
56/// \param[in] title Title for plotting etc.
57/// \param[in] func Function
58/// \param[in] hdata Data histogram
59/// \param[in] argX Optional arguments according to table below.
60/// <table>
61/// <tr><th> Argument <th> Effect
62/// <tr><td>
63/// DataError() <td> Choose between Poisson errors and Sum-of-weights errors
64/// <tr><td>
65/// NumCPU() <td> Activate parallel processing feature
66/// <tr><td>
67/// Range() <td> Fit only selected region
68/// <tr><td>
69/// Verbose() <td> Verbose output of GOF framework
70
71RooChi2Var::RooChi2Var(const char *name, const char* title, RooAbsReal& func, RooDataHist& hdata,
72 const RooCmdArg& arg1,const RooCmdArg& arg2,const RooCmdArg& arg3,
73 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,
74 const RooCmdArg& arg7,const RooCmdArg& arg8,const RooCmdArg& arg9) :
75 RooAbsOptTestStatistic(name,title,func,hdata,_emptySet,
76 RooCmdConfig::decodeStringOnTheFly("RooChi2Var::RooChi2Var","RangeWithName",0,"",arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
77 0,
78 RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","NumCPU",0,1,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
80 RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","Verbose",0,1,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
81 0)
82{
83 RooCmdConfig pc("RooChi2Var::RooChi2Var") ;
84 pc.defineInt("etype","DataError",0,(Int_t)RooDataHist::Auto) ;
85 pc.defineInt("extended","Extended",0,kFALSE) ;
86 pc.allowUndefined() ;
87
88 pc.process(arg1) ; pc.process(arg2) ; pc.process(arg3) ;
89 pc.process(arg4) ; pc.process(arg5) ; pc.process(arg6) ;
90 pc.process(arg7) ; pc.process(arg8) ; pc.process(arg9) ;
91
92 if (func.IsA()->InheritsFrom(RooAbsPdf::Class())) {
93 _funcMode = pc.getInt("extended") ? ExtendedPdf : Pdf ;
94 } else {
96 }
97 _etype = (RooDataHist::ErrorType) pc.getInt("etype") ;
98
101 }
102
103}
104
105
106
107////////////////////////////////////////////////////////////////////////////////
108/// RooChi2Var constructor. Optional arguments taken
109///
110/// \param[in] name Name of the PDF
111/// \param[in] title Title for plotting etc.
112/// \param[in] pdf PDF to fit
113/// \param[in] hdata Data histogram
114/// \param[in] argX Optional arguments according to table below.
115/// <table>
116/// <tr><th> Argument <th> Effect
117/// <tr><td>
118/// Extended() <td> Include extended term in calculation
119/// <tr><td>
120/// DataError() <td> Choose between Poisson errors and Sum-of-weights errors
121/// <tr><td>
122/// NumCPU() <td> Activate parallel processing feature
123/// <tr><td>
124/// Range() <td> Fit only selected region
125/// <tr><td>
126/// SumCoefRange() <td> Set the range in which to interpret the coefficients of RooAddPdf components
127/// <tr><td>
128/// SplitRange() <td> Fit range is split by index catory of simultaneous PDF
129/// <tr><td>
130/// ConditionalObservables() <td> Define projected observables
131/// <tr><td>
132/// Verbose() <td> Verbose output of GOF framework
133
134RooChi2Var::RooChi2Var(const char *name, const char* title, RooAbsPdf& pdf, RooDataHist& hdata,
135 const RooCmdArg& arg1,const RooCmdArg& arg2,const RooCmdArg& arg3,
136 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,
137 const RooCmdArg& arg7,const RooCmdArg& arg8,const RooCmdArg& arg9) :
138 RooAbsOptTestStatistic(name,title,pdf,hdata,
139 *(const RooArgSet*)RooCmdConfig::decodeObjOnTheFly("RooChi2Var::RooChi2Var","ProjectedObservables",0,&_emptySet
140 ,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
141 RooCmdConfig::decodeStringOnTheFly("RooChi2Var::RooChi2Var","RangeWithName",0,"",arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
142 RooCmdConfig::decodeStringOnTheFly("RooChi2Var::RooChi2Var","AddCoefRange",0,"",arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
143 RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","NumCPU",0,1,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
145 RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","Verbose",0,1,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
146 RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","SplitRange",0,0,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9))
147{
148 RooCmdConfig pc("RooChi2Var::RooChi2Var") ;
149 pc.defineInt("extended","Extended",0,kFALSE) ;
150 pc.defineInt("etype","DataError",0,(Int_t)RooDataHist::Auto) ;
151 pc.allowUndefined() ;
152
153 pc.process(arg1) ; pc.process(arg2) ; pc.process(arg3) ;
154 pc.process(arg4) ; pc.process(arg5) ; pc.process(arg6) ;
155 pc.process(arg7) ; pc.process(arg8) ; pc.process(arg9) ;
156
157 _funcMode = pc.getInt("extended") ? ExtendedPdf : Pdf ;
158 _etype = (RooDataHist::ErrorType) pc.getInt("etype") ;
161 }
162}
163
164
165
166////////////////////////////////////////////////////////////////////////////////
167/// Constructor of a chi2 for given p.d.f. with respect given binned
168/// dataset. If cutRange is specified the calculation of the chi2 is
169/// restricted to that named range. If addCoefRange is specified, the
170/// interpretation of fractions for all component RooAddPdfs that do
171/// not have a frozen range interpretation is set to chosen range
172/// name. If nCPU is greater than one the chi^2 calculation is
173/// paralellized over the specified number of processors. If
174/// interleave is true the partitioning of event over processors
175/// follows a (i % n == i_set) strategy rather than a bulk
176/// partitioning strategy which may result in unequal load balancing
177/// in binned datasets with many (adjacent) zero bins. If
178/// splitCutRange is true the cutRange is used to construct an
179/// individual cutRange for each RooSimultaneous index category state
180/// name cutRange_{indexStateName}.
181
182RooChi2Var::RooChi2Var(const char *name, const char *title, RooAbsPdf& pdf, RooDataHist& hdata,
183 Bool_t extended, const char* cutRange, const char* addCoefRange,
184 Int_t nCPU, RooFit::MPSplit interleave, Bool_t verbose, Bool_t splitCutRange, RooDataHist::ErrorType etype) :
185 RooAbsOptTestStatistic(name,title,pdf,hdata,RooArgSet(),cutRange,addCoefRange,nCPU,interleave,verbose,splitCutRange),
186 _etype(etype), _funcMode(extended?ExtendedPdf:Pdf)
187{
188}
189
190
191
192////////////////////////////////////////////////////////////////////////////////
193/// Constructor of a chi2 for given p.d.f. with respect given binned
194/// dataset taking the observables specified in projDeps as projected
195/// observables. If cutRange is specified the calculation of the chi2
196/// is restricted to that named range. If addCoefRange is specified,
197/// the interpretation of fractions for all component RooAddPdfs that
198/// do not have a frozen range interpretation is set to chosen range
199/// name. If nCPU is greater than one the chi^2 calculation is
200/// paralellized over the specified number of processors. If
201/// interleave is true the partitioning of event over processors
202/// follows a (i % n == i_set) strategy rather than a bulk
203/// partitioning strategy which may result in unequal load balancing
204/// in binned datasets with many (adjacent) zero bins. If
205/// splitCutRange is true the cutRange is used to construct an
206/// individual cutRange for each RooSimultaneous index category state
207/// name cutRange_{indexStateName}.
208
209RooChi2Var::RooChi2Var(const char *name, const char *title, RooAbsReal& func, RooDataHist& hdata,
210 const RooArgSet& projDeps, RooChi2Var::FuncMode fmode, const char* cutRange, const char* addCoefRange,
211 Int_t nCPU, RooFit::MPSplit interleave, Bool_t verbose, Bool_t splitCutRange, RooDataHist::ErrorType etype) :
212 RooAbsOptTestStatistic(name,title,func,hdata,projDeps,cutRange,addCoefRange,nCPU,interleave,verbose,splitCutRange),
213 _etype(etype), _funcMode(fmode)
214{
215}
216
217
218
219////////////////////////////////////////////////////////////////////////////////
220/// Copy constructor
221
222RooChi2Var::RooChi2Var(const RooChi2Var& other, const char* name) :
224 _etype(other._etype),
225 _funcMode(other._funcMode)
226{
227}
228
229
230
231////////////////////////////////////////////////////////////////////////////////
232/// Destructor
233
235{
236}
237
238
239
240////////////////////////////////////////////////////////////////////////////////
241/// Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize
242
243Double_t RooChi2Var::evaluatePartition(Int_t firstEvent, Int_t lastEvent, Int_t stepSize) const
244{
245 // Throughout the calculation, we use Kahan's algorithm for summing to
246 // prevent loss of precision - this is a factor four more expensive than
247 // straight addition, but since evaluating the PDF is usually much more
248 // expensive than that, we tolerate the additional cost...
249 Int_t i ;
250 Double_t result(0), carry(0);
251
252 _dataClone->store()->recalculateCache( _projDeps, firstEvent, lastEvent, stepSize, kFALSE) ;
253
254
255 // Determine normalization factor depending on type of input function
256 Double_t normFactor(1) ;
257 switch (_funcMode) {
258 case Function: normFactor=1 ; break ;
259 case Pdf: normFactor = _dataClone->sumEntries() ; break ;
260 case ExtendedPdf: normFactor = ((RooAbsPdf*)_funcClone)->expectedEvents(_dataClone->get()) ; break ;
261 }
262
263 // Loop over bins of dataset
265 for (i=firstEvent ; i<lastEvent ; i+=stepSize) {
266
267 // get the data values for this event
268 hdata->get(i);
269
270 if (!hdata->valid()) continue;
271
272 const Double_t nData = hdata->weight() ;
273
274 const Double_t nPdf = _funcClone->getVal(_normSet) * normFactor * hdata->binVolume() ;
275
276 const Double_t eExt = nPdf-nData ;
277
278
279 Double_t eInt ;
281 Double_t eIntLo,eIntHi ;
282 hdata->weightError(eIntLo,eIntHi,_etype) ;
283 eInt = (eExt>0) ? eIntHi : eIntLo ;
284 } else {
285 eInt = sqrt(nPdf) ;
286 }
287
288 // Skip cases where pdf=0 and there is no data
289 if (0. == eInt * eInt && 0. == nData * nData && 0. == nPdf * nPdf) continue ;
290
291 // Return 0 if eInt=0, special handling in MINUIT will follow
292 if (0. == eInt * eInt) {
293 coutE(Eval) << "RooChi2Var::RooChi2Var(" << GetName() << ") INFINITY ERROR: bin " << i
294 << " has zero error" << endl ;
295 return 0.;
296 }
297
298// cout << "Chi2Var[" << i << "] nData = " << nData << " nPdf = " << nPdf << " errorExt = " << eExt << " errorInt = " << eInt << " contrib = " << eExt*eExt/(eInt*eInt) << endl ;
299
300 Double_t term = eExt*eExt/(eInt*eInt) ;
301 Double_t y = term - carry;
302 Double_t t = result + y;
303 carry = (t - result) - y;
304 result = t;
305 }
306
307 _evalCarry = carry;
308 return result ;
309}
310
311
312
void Class()
Definition: Class.C:29
#define coutE(a)
Definition: RooMsgService.h:34
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
double sqrt(double)
virtual void recalculateCache(const RooArgSet *, Int_t, Int_t, Int_t, Bool_t)
virtual const RooArgSet * get() const
Definition: RooAbsData.h:80
RooAbsDataStore * store()
Definition: RooAbsData.h:56
virtual Double_t sumEntries() const =0
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:53
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:81
Double_t _evalCarry
avoids loss of precision
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Class RooChi2Var implements a simple chi^2 calculation from a binned dataset and a PDF.
Definition: RooChi2Var.h:25
RooDataHist::ErrorType _etype
Definition: RooChi2Var.h:71
virtual Double_t evaluatePartition(Int_t firstEvent, Int_t lastEvent, Int_t stepSize) const
Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize.
Definition: RooChi2Var.cxx:243
FuncMode _funcMode
Definition: RooChi2Var.h:72
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:71
static RooArgSet _emptySet
Definition: RooChi2Var.h:69
virtual ~RooChi2Var()
Destructor.
Definition: RooChi2Var.cxx:234
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
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
virtual Double_t weight() const
Definition: RooDataHist.h:98
virtual void weightError(Double_t &lo, Double_t &hi, ErrorType etype=Poisson) const
Return the error on current weight.
virtual Bool_t isNonPoissonWeighted() const
Returns true if datasets contains entries with a non-integer weight.
virtual const RooArgSet * get() const
Definition: RooDataHist.h:79
virtual Bool_t valid() const
Return true if currently loaded coordinate is considered valid within the current range definitions o...
Double_t binVolume() const
Definition: RooDataHist.h:104
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:443
Double_t y[n]
Definition: legend1.C:17
const UInt_t eInt[256]
Template specialisation used in RooAbsArg:
@ Interleave
Definition: RooGlobalFunc.h:60
static constexpr double pc