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 "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
66using namespace std;
67
68namespace {
69 template<class ...Args>
70 RooAbsTestStatistic::Configuration makeRooAbsTestStatisticCfgForFunc(Args const& ... args) {
72 cfg.rangeName = RooCmdConfig::decodeStringOnTheFly("RooChi2Var::RooChi2Var","RangeWithName",0,"",args...);
73 cfg.nCPU = RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","NumCPU",0,1,args...);
75 cfg.verbose = static_cast<bool>(RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","Verbose",0,1,args...));
76 cfg.cloneInputData = false;
77 cfg.integrateOverBinsPrecision = RooCmdConfig::decodeDoubleOnTheFly("RooChi2Var::RooChi2Var", "IntegrateBins", 0, -1., {args...});
78 return cfg;
79 }
80
81 template<class ...Args>
82 RooAbsTestStatistic::Configuration makeRooAbsTestStatisticCfgForPdf(Args const& ... args) {
83 auto cfg = makeRooAbsTestStatisticCfgForFunc(args...);
84 cfg.addCoefRangeName = RooCmdConfig::decodeStringOnTheFly("RooChi2Var::RooChi2Var","AddCoefRange",0,"",args...);
85 cfg.splitCutRange = static_cast<bool>(RooCmdConfig::decodeIntOnTheFly("RooChi2Var::RooChi2Var","SplitRange",0,0,args...));
86 return cfg;
87 }
88}
89
91;
92
94
95
96////////////////////////////////////////////////////////////////////////////////
97/// RooChi2Var constructor. Optional arguments are:
98/// \param[in] name Name of the PDF
99/// \param[in] title Title for plotting etc.
100/// \param[in] func Function
101/// \param[in] hdata Data histogram
102/// \param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9 Optional arguments according to table below.
103/// <table>
104/// <tr><th> Type of CmdArg <th> Effect on \f$ \chi^2 \f$
105/// <tr><td>
106/// <tr><td> `DataError()` <td> Choose between:
107/// - RooAbsData::Expected: Expected Poisson error (\f$ \sqrt{n_\text{expected}} \f$ from the PDF).
108/// - RooAbsData::SumW2: The observed error from the square root of the sum of weights,
109/// i.e., symmetric errors calculated with the standard deviation of a Poisson distribution.
110/// - RooAbsData::Poisson: Asymmetric errors from the central 68 % interval around a Poisson distribution with mean \f$ n_\text{observed} \f$.
111/// If for a given bin \f$ n_\text{expected} \f$ is lower than the \f$ n_\text{observed} \f$, the lower uncertainty is taken
112/// (e.g., the difference between the mean and the 16 % quantile).
113/// If \f$ n_\text{expected} \f$ is higher than \f$ n_\text{observed} \f$, the higher uncertainty is taken
114/// (e.g., the difference between the 84 % quantile and the mean).
115/// - RooAbsData::Auto (default): RooAbsData::Expected for unweighted data, RooAbsData::SumW2 for weighted data.
116/// <tr><td>
117/// `Extended()` <td> Use expected number of events of an extended p.d.f as normalization
118/// <tr><td>
119/// NumCPU() <td> Activate parallel processing feature
120/// <tr><td>
121/// Range() <td> Calculate \f$ \chi^2 \f$ only in selected region
122/// <tr><td>
123/// Verbose() <td> Verbose output of GOF framework
124/// <tr><td>
125/// IntegrateBins() <td> Integrate PDF within each bin. This sets the desired precision. Only useful for binned fits.
126/// <tr><td> `SumCoefRange()` <td> Set the range in which to interpret the coefficients of RooAddPdf components
127/// <tr><td> `SplitRange()` <td> Fit ranges used in different categories get named after the category.
128/// Using `Range("range"), SplitRange()` as switches, different ranges could be set like this:
129/// ```
130/// myVariable.setRange("range_pi0", 135, 210);
131/// myVariable.setRange("range_gamma", 50, 210);
132/// ```
133/// <tr><td> `ConditionalObservables(Args_t &&... argsOrArgSet)` <td> Define projected observables.
134/// Arguments can either be multiple RooRealVar or a single RooArgSet containing them.
135///
136/// </table>
137
138RooChi2Var::RooChi2Var(const char *name, const char* title, RooAbsReal& func, RooDataHist& hdata,
139 const RooCmdArg& arg1,const RooCmdArg& arg2,const RooCmdArg& arg3,
140 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,
141 const RooCmdArg& arg7,const RooCmdArg& arg8,const RooCmdArg& arg9) :
142 RooAbsOptTestStatistic(name,title,func,hdata,_emptySet,
143 makeRooAbsTestStatisticCfgForFunc(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9))
144{
145 RooCmdConfig pc("RooChi2Var::RooChi2Var") ;
146 pc.defineInt("etype","DataError",0,(Int_t)RooDataHist::Auto) ;
147 pc.defineInt("extended","Extended",0,false) ;
148 pc.allowUndefined() ;
149
150 pc.process(arg1) ; pc.process(arg2) ; pc.process(arg3) ;
151 pc.process(arg4) ; pc.process(arg5) ; pc.process(arg6) ;
152 pc.process(arg7) ; pc.process(arg8) ; pc.process(arg9) ;
153
154 if (func.IsA()->InheritsFrom(RooAbsPdf::Class())) {
155 _funcMode = pc.getInt("extended") ? ExtendedPdf : Pdf ;
156 } else {
158 }
159 _etype = (RooDataHist::ErrorType) pc.getInt("etype") ;
160
163 }
164
165}
166
167
168////////////////////////////////////////////////////////////////////////////////
169/// RooChi2Var constructor. Optional arguments taken
170///
171/// \param[in] name Name of the PDF
172/// \param[in] title Title for plotting etc.
173/// \param[in] pdf PDF to fit
174/// \param[in] hdata Data histogram
175/// \param[in] arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9 Optional arguments according to table below.
176/// <table>
177/// <tr><th> Argument <th> Effect
178/// <tr><td>
179/// Extended() <td> Include extended term in calculation
180/// <tr><td>
181/// DataError() <td> Choose between Poisson errors and Sum-of-weights errors
182/// <tr><td>
183/// NumCPU() <td> Activate parallel processing feature
184/// <tr><td>
185/// Range() <td> Fit only selected region
186/// <tr><td>
187/// SumCoefRange() <td> Set the range in which to interpret the coefficients of RooAddPdf components
188/// <tr><td>
189/// SplitRange() <td> Fit range is split by index category of simultaneous PDF
190/// <tr><td>
191/// ConditionalObservables() <td> Define projected observables
192/// <tr><td>
193/// Verbose() <td> Verbose output of GOF framework
194/// <tr><td>
195/// IntegrateBins() <td> Integrate PDF within each bin. This sets the desired precision.
196
197RooChi2Var::RooChi2Var(const char *name, const char* title, RooAbsPdf& pdf, RooDataHist& hdata,
198 const RooCmdArg& arg1,const RooCmdArg& arg2,const RooCmdArg& arg3,
199 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,
200 const RooCmdArg& arg7,const RooCmdArg& arg8,const RooCmdArg& arg9) :
201 RooAbsOptTestStatistic(name,title,pdf,hdata,
202 *RooCmdConfig::decodeSetOnTheFly("RooChi2Var::RooChi2Var","ProjectedObservables",0,&_emptySet,
203 arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9),
204 makeRooAbsTestStatisticCfgForPdf(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9))
205{
206 RooCmdConfig pc("RooChi2Var::RooChi2Var") ;
207 pc.defineInt("extended","Extended",0,false) ;
208 pc.defineInt("etype","DataError",0,(Int_t)RooDataHist::Auto) ;
209 pc.allowUndefined() ;
210
211 pc.process(arg1) ; pc.process(arg2) ; pc.process(arg3) ;
212 pc.process(arg4) ; pc.process(arg5) ; pc.process(arg6) ;
213 pc.process(arg7) ; pc.process(arg8) ; pc.process(arg9) ;
214
215 _funcMode = pc.getInt("extended") ? ExtendedPdf : Pdf ;
216 _etype = (RooDataHist::ErrorType) pc.getInt("etype") ;
219 }
220}
221
222
223////////////////////////////////////////////////////////////////////////////////
224/// Copy constructor
225
226RooChi2Var::RooChi2Var(const RooChi2Var& other, const char* name) :
228 _etype(other._etype),
229 _funcMode(other._funcMode)
230{
231}
232
233
234////////////////////////////////////////////////////////////////////////////////
235/// Calculate chi^2 in partition from firstEvent to lastEvent using given stepSize
236/// Throughout the calculation, we use Kahan's algorithm for summing to
237/// prevent loss of precision - this is a factor four more expensive than
238/// straight addition, but since evaluating the PDF is usually much more
239/// expensive than that, we tolerate the additional cost...
240
241double RooChi2Var::evaluatePartition(std::size_t firstEvent, std::size_t lastEvent, std::size_t stepSize) const
242{
243
244 double result(0), carry(0);
245
246 _dataClone->store()->recalculateCache( _projDeps, firstEvent, lastEvent, stepSize, false) ;
247
248
249 // Determine normalization factor depending on type of input function
250 double normFactor(1) ;
251 switch (_funcMode) {
252 case Function: normFactor=1 ; break ;
253 case Pdf: normFactor = _dataClone->sumEntries() ; break ;
254 case ExtendedPdf: normFactor = ((RooAbsPdf*)_funcClone)->expectedEvents(_dataClone->get()) ; break ;
255 }
256
257 // Loop over bins of dataset
259 for (auto i=firstEvent ; i<lastEvent ; i+=stepSize) {
260
261 // get the data values for this event
262 hdata->get(i);
263
264 const double nData = hdata->weight() ;
265
266 const double nPdf = _funcClone->getVal(_normSet) * normFactor * hdata->binVolume() ;
267
268 const double eExt = nPdf-nData ;
269
270
271 double eInt ;
273 double eIntLo,eIntHi ;
274 hdata->weightError(eIntLo,eIntHi,_etype) ;
275 eInt = (eExt>0) ? eIntHi : eIntLo ;
276 } else {
277 eInt = sqrt(nPdf) ;
278 }
279
280 // Skip cases where pdf=0 and there is no data
281 if (0. == eInt * eInt && 0. == nData * nData && 0. == nPdf * nPdf) continue ;
282
283 // Return 0 if eInt=0, special handling in MINUIT will follow
284 if (0. == eInt * eInt) {
285 coutE(Eval) << "RooChi2Var::RooChi2Var(" << GetName() << ") INFINITY ERROR: bin " << i
286 << " has zero error" << endl ;
287 return 0.;
288 }
289
290// cout << "Chi2Var[" << i << "] nData = " << nData << " nPdf = " << nPdf << " errorExt = " << eExt << " errorInt = " << eInt << " contrib = " << eExt*eExt/(eInt*eInt) << endl ;
291
292 double term = eExt*eExt/(eInt*eInt) ;
293 double y = term - carry;
294 double t = result + y;
295 carry = (t - result) - y;
296 result = t;
297 }
298
299 _evalCarry = carry;
300 return result ;
301}
#define coutE(a)
Definition: RooMsgService.h:37
#define ClassImp(name)
Definition: Rtypes.h:375
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
virtual void recalculateCache(const RooArgSet *, Int_t, Int_t, Int_t, bool)
virtual double sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
virtual const RooArgSet * get() const
Definition: RooAbsData.h:106
RooAbsDataStore * store()
Definition: RooAbsData.h:82
RooAbsOptTestStatistic is the abstract base class for test statistics objects that evaluate a functio...
RooAbsReal * _funcClone
Pointer to internal clone of input function.
RooArgSet * _normSet
Pointer to set with observables used for normalization.
RooAbsData * _dataClone
Pointer to internal clone if input data.
RooArgSet * _projDeps
Set of projected observable.
static TClass * Class()
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
TClass * IsA() const override
Definition: RooAbsReal.h:581
double _evalCarry
! carry of Kahan sum in evaluatePartition
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooChi2Var implements a simple calculation from a binned dataset and a PDF.
Definition: RooChi2Var.h:25
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 Throughout the calcula...
Definition: RooChi2Var.cxx:241
RooDataHist::ErrorType _etype
Error type store in associated RooDataHist.
Definition: RooChi2Var.h:69
RooChi2Var(const char *name, const char *title, RooAbsReal &func, RooDataHist &data, const RooCmdArg &arg1=RooCmdArg::none(), 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:138
FuncMode _funcMode
Function, P.d.f. or extended p.d.f?
Definition: RooChi2Var.h:70
static RooArgSet _emptySet
Supports named argument constructor.
Definition: RooChi2Var.h:67
RooCmdArg is a named container for two doubles, two integers two object points and three string point...
Definition: RooCmdArg.h:26
Class RooCmdConfig is a configurable parser for RooCmdArg named arguments.
Definition: RooCmdConfig.h:31
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.
static std::string decodeStringOnTheFly(const char *callerID, const char *cmdArgName, Int_t intIdx, const char *defVal, Args_t &&...args)
Static decoder function allows to retrieve string property from set of RooCmdArgs For use in base mem...
Definition: RooCmdConfig.h:204
static Int_t decodeIntOnTheFly(const char *callerID, const char *cmdArgName, Int_t intIdx, Int_t defVal, Args_t &&...args)
Static decoder function allows to retrieve integer property from set of RooCmdArgs For use in base me...
Definition: RooCmdConfig.h:188
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:104
void weightError(double &lo, double &hi, ErrorType etype=Poisson) const override
Return the asymmetric errors on the current weight.
bool isNonPoissonWeighted() const override
Returns true if dataset contains entries with a non-integer weight.
double binVolume(std::size_t i) const
Return bin volume of i-th bin.
Definition: RooDataHist.h:112
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:76
Bool_t InheritsFrom(const char *cl) const override
Return kTRUE if this class inherits from a class with name "classname".
Definition: TClass.cxx:4863
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Double_t y[n]
Definition: legend1.C:17
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
const UInt_t eInt[256]
@ Interleave
Definition: RooGlobalFunc.h:64
static constexpr double pc
std::string rangeName
Stores the configuration parameters for RooAbsTestStatistic.