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