Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooParamHistFunc.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 *
4 * Copyright (c) 2023, CERN
5 *
6 * Redistribution and use in source and binary forms,
7 * with or without modification, are permitted according to the terms
8 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
9 */
10
11/** \class RooParamHistFunc
12 * \ingroup Roofit
13 * A histogram function that assigns scale parameters to every bin. Instead of the bare bin contents,
14 * it therefore yields:
15 * \f[
16 * \gamma_{i} * \mathrm{bin}_i
17 * \f]
18 *
19 * The \f$ \gamma_i \f$ can therefore be used to parametrise statistical uncertainties of the histogram
20 * template. In conjunction with a constraint term, this can be used to implement the Barlow-Beeston method.
21 * The constraint can be implemented using RooHistConstraint.
22 *
23 * See also the tutorial rf709_BarlowBeeston.C
24 */
25
26#include "Riostream.h"
27#include "RooParamHistFunc.h"
28#include "RooAbsCategory.h"
29#include "RooRealVar.h"
30#include "RooFitImplHelpers.h"
31#include <cmath>
32#include "TMath.h"
33
34
35using namespace std ;
36
38
39////////////////////////////////////////////////////////////////////////////////
40/// Populate x with observables
41
42RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataHist& dh, bool paramRelative) :
43 RooAbsReal(name,title),
44 _x("x","x",this),
45 _p("p","p",this),
46 _dh(dh),
47 _relParam(paramRelative)
48{
49 _x.add(*_dh.get()) ;
50
51 // Now populate p with parameters
52 RooArgSet allVars ;
53 for (Int_t i=0 ; i<_dh.numEntries() ; i++) {
54 _dh.get(i) ;
55
56 const char* vname = Form("%s_gamma_bin_%i",GetName(),i) ;
57 RooRealVar* var = new RooRealVar(vname,vname,0,1000) ;
58 var->setVal(_relParam ? 1 : _dh.weight()) ;
59 var->setError(_relParam ? 1 / sqrt(_dh.weight()) : sqrt(_dh.weight()));
60 var->setConstant(true) ;
61 allVars.add(*var) ;
62 _p.add(*var) ;
63 }
64 addOwnedComponents(allVars) ;
65}
66
67////////////////////////////////////////////////////////////////////////////////
68/// Populate x with observables
69
70RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, const RooAbsArg& /*x*/, RooDataHist& dh, bool paramRelative) :
71 RooAbsReal(name,title),
72 _x("x","x",this),
73 _p("p","p",this),
74 _dh(dh),
75 _relParam(paramRelative)
76{
77 _x.add(*_dh.get()) ;
78
79 // Now populate p with parameters
80 RooArgSet allVars ;
81 for (Int_t i=0 ; i<_dh.numEntries() ; i++) {
82 _dh.get(i) ;
83 const char* vname = Form("%s_gamma_bin_%i",GetName(),i) ;
84 RooRealVar* var = new RooRealVar(vname,vname,0,1000) ;
85 var->setVal(_relParam ? 1 : _dh.weight()) ;
86 var->setError(_relParam ? 1 / sqrt(_dh.weight()) : sqrt(_dh.weight()));
87 var->setConstant(true) ;
88 allVars.add(*var) ;
89 _p.add(*var) ;
90 }
91 addOwnedComponents(allVars) ;
92}
93
94////////////////////////////////////////////////////////////////////////////////
95
96RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataHist& dh, const RooParamHistFunc& paramSource, bool paramRelative) :
97 RooAbsReal(name,title),
98 _x("x","x",this),
99 _p("p","p",this),
100 _dh(dh),
101 _relParam(paramRelative)
102{
103 // Populate x with observables
104 _x.add(*_dh.get()) ;
105
106 // Now populate p with existing parameters
107 _p.add(paramSource._p) ;
108}
109
110////////////////////////////////////////////////////////////////////////////////
111
113 RooAbsReal(other,name),
114 _x("x",this,other._x),
115 _p("p",this,other._p),
116 _dh(other._dh),
117 _relParam(other._relParam)
118{
119}
120
121////////////////////////////////////////////////////////////////////////////////
122
124{
125 Int_t idx = ((RooDataHist&)_dh).getIndex(_x,true) ;
126 double ret = (static_cast<RooAbsReal*>(_p.at(idx)))->getVal() ;
127 if (_relParam) {
128 double nom = getNominal(idx) ;
129 ret *= nom ;
130 }
131 return ret ;
132}
133
135{
136 std::string const &idx = _dh.calculateTreeIndexForCodeSquash(this, ctx, _x);
137 std::string arrName = ctx.buildArg(_p);
138 std::string result = arrName + "[" + idx + "]";
139 if (_relParam) {
140 // get weight[idx] * binv[idx]. Here we get the bin volume for the first element as we assume the distribution to
141 // be binned uniformly.
142 std::string const &weightName = _dh.declWeightArrayForCodeSquash(this, ctx, false);
143 std::string nominalVal = weightName + "[" + idx + "] * " + std::to_string(_dh.binVolume(0));
144 result += " * " + nominalVal;
145 }
146 ctx.addResult(this, result);
147}
148
149////////////////////////////////////////////////////////////////////////////////
150
152{
153 return (static_cast<RooAbsReal&>(_p[ibin])).getVal() ;
154}
155
156////////////////////////////////////////////////////////////////////////////////
157
158void RooParamHistFunc::setActual(Int_t ibin, double newVal)
159{
160 (static_cast<RooRealVar&>(_p[ibin])).setVal(newVal) ;
161}
162
163////////////////////////////////////////////////////////////////////////////////
164
166{
167 _dh.get(ibin) ;
168 return _dh.weight() ;
169}
170
171////////////////////////////////////////////////////////////////////////////////
172
174{
175 _dh.get(ibin) ;
176 return _dh.weightError() ;
177}
178
179////////////////////////////////////////////////////////////////////////////////
180/// Return sampling hint for making curves of (projections) of this function
181/// as the recursive division strategy of RooCurve cannot deal efficiently
182/// with the vertical lines that occur in a non-interpolated histogram
183
184list<double>* RooParamHistFunc::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const
185{
186 // Check that observable is in dataset, if not no hint is generated
187 RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(_dh.get()->find(obs.GetName())) ;
188 if (!lvarg) {
189 return nullptr ;
190 }
191
192 // Retrieve position of all bin boundaries
193 const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr);
194 double* boundaries = binning->array() ;
195
196 list<double>* hint = new list<double> ;
197
198 // Widen range slightly
199 xlo = xlo - 0.01*(xhi-xlo) ;
200 xhi = xhi + 0.01*(xhi-xlo) ;
201
202 double delta = (xhi-xlo)*1e-8 ;
203
204 // Construct array with pairs of points positioned epsilon to the left and
205 // right of the bin boundaries
206 for (Int_t i=0 ; i<binning->numBoundaries() ; i++) {
207 if (boundaries[i]>=xlo && boundaries[i]<=xhi) {
208 hint->push_back(boundaries[i]-delta) ;
209 hint->push_back(boundaries[i]+delta) ;
210 }
211 }
212
213 return hint ;
214}
215
216////////////////////////////////////////////////////////////////////////////////
217/// Return sampling hint for making curves of (projections) of this function
218/// as the recursive division strategy of RooCurve cannot deal efficiently
219/// with the vertical lines that occur in a non-interpolated histogram
220
221std::list<double>* RooParamHistFunc::binBoundaries(RooAbsRealLValue& obs, double xlo, double xhi) const
222{
223 // Check that observable is in dataset, if not no hint is generated
224 RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(_dh.get()->find(obs.GetName())) ;
225 if (!lvarg) {
226 return nullptr ;
227 }
228
229 // Retrieve position of all bin boundaries
230 const RooAbsBinning* binning = lvarg->getBinningPtr(nullptr);
231 double* boundaries = binning->array() ;
232
233 list<double>* hint = new list<double> ;
234
235 // Construct array with pairs of points positioned epsilon to the left and
236 // right of the bin boundaries
237 for (Int_t i=0 ; i<binning->numBoundaries() ; i++) {
238 if (boundaries[i]>=xlo && boundaries[i]<=xhi) {
239 hint->push_back(boundaries[i]) ;
240 }
241 }
242
243 return hint ;
244}
245
246////////////////////////////////////////////////////////////////////////////////
247/// Advertise that all integrals can be handled internally.
248
250 const RooArgSet* /*normSet*/, const char* /*rangeName*/) const
251{
252 // Simplest scenario, integrate over all dependents
253 std::unique_ptr<RooAbsCollection> allVarsCommon{allVars.selectCommon(_x)};
254 bool intAllObs = (allVarsCommon->size()==_x.size()) ;
255 if (intAllObs && matchArgs(allVars,analVars,_x)) {
256 return 1 ;
257 }
258
259 return 0 ;
260}
261
262////////////////////////////////////////////////////////////////////////////////
263/// Implement analytical integrations by doing appropriate weighting from component integrals
264/// functions to integrators of components
265
266double RooParamHistFunc::analyticalIntegralWN(Int_t code, const RooArgSet* /*normSet2*/,const char* rangeName) const
267{
268 // Supports only the scenario of integration over all dependents
269 R__ASSERT(code==1) ;
270
271 // The logic for summing over the histogram is borrowed from RooHistPdf with some differences:
272 //
273 // - a lambda function is used to inject the parameters for bin scaling into the RooDataHist::sum method
274 //
275 // - for simplicity, there is no check for the possibility of full-range integration with another overload of
276 // RooDataHist::sum
277 std::map<const RooAbsArg*, std::pair<double, double> > ranges;
278 for (const auto obs : _x) {
279 ranges[obs] = RooHelpers::getRangeOrBinningInterval(obs, rangeName);
280 }
281
282 auto getBinScale = [&](int iBin){ return static_cast<const RooAbsReal&>(_p[iBin]).getVal(); };
283
284 RooArgSet sliceSet{};
285 return const_cast<RooDataHist&>(_dh).sum(_x, sliceSet, true, false, ranges, getBinScale);
286}
#define e(i)
Definition RSha256.hxx:103
#define ClassImp(name)
Definition Rtypes.h:377
#define R__ASSERT(e)
Definition TError.h:118
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
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
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Abstract base class for RooRealVar binning definitions.
virtual Int_t numBoundaries() const =0
virtual double * array() const =0
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
bool selectCommon(const RooAbsCollection &refColl, RooAbsCollection &outColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Abstract base class for objects that are lvalues, i.e.
virtual const RooAbsBinning * getBinningPtr(const char *rangeName) const =0
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
void setConstant(bool value=true)
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:39
std::string calculateTreeIndexForCodeSquash(RooAbsArg const *klass, RooFit::Detail::CodeSquashContext &ctx, const RooAbsCollection &coords, bool reverse=false) const
double weight(std::size_t i) const
Return weight of i-th bin.
void weightError(double &lo, double &hi, ErrorType etype=Poisson) const override
Return the asymmetric errors on the current weight.
std::string declWeightArrayForCodeSquash(RooAbsArg const *klass, RooFit::Detail::CodeSquashContext &ctx, bool correctForBinSize) const
double binVolume(std::size_t i) const
Return bin volume of i-th bin.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition RooDataHist.h:76
A class to maintain the context for squashing of RooFit models into code.
void addResult(RooAbsArg const *key, std::string const &value)
A function to save an expression that includes/depends on the result of the input node.
std::string buildArg(RooAbsCollection const &x)
Function to save a RooListProxy as an array in the squashed code.
A histogram function that assigns scale parameters to every bin.
void setActual(Int_t ibin, double newVal)
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
double getNominalError(Int_t ibin) const
double getActual(Int_t ibin)
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Implement analytical integrations by doing appropriate weighting from component integrals functions t...
double getNominal(Int_t ibin) const
void translate(RooFit::Detail::CodeSquashContext &ctx) const override
This function defines a translation for each RooAbsReal based object that can be used to express the ...
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Return sampling hint for making curves of (projections) of this function as the recursive division st...
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &analVars, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Advertise that all integrals can be handled internally.
std::list< double > * plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const override
Return sampling hint for making curves of (projections) of this function as the recursive division st...
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
void setError(double value)
Definition RooRealVar.h:60
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
std::pair< double, double > getRangeOrBinningInterval(RooAbsArg const *arg, const char *rangeName)
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345