Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooParametricStepFunction.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitBabar *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * Aaron Roodman, Stanford Linear Accelerator Center, Stanford University *
7 * *
8 * Copyright (c) 2004, Stanford University. All rights reserved. *
9 *
10 * Redistribution and use in source and binary forms, *
11 * with or without modification, are permitted according to the terms *
12 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
13 *****************************************************************************/
14
15/** \class RooParametricStepFunction
16 \ingroup Roofit
17
18The Parametric Step Function PDF is a binned distribution whose parameters
19are the heights of each bin. This PDF was first used in BaBar's B0->pi0pi0
20paper BABAR Collaboration (B. Aubert et al.) Phys.Rev.Lett.91:241801,2003.
21
22This PDF may be used to describe oddly shaped distributions. It differs
23from a RooKeysPdf or a RooHistPdf in that a RooParametricStepFunction
24has free parameters. In particular, any statistical uncertainty in
25sample used to model this PDF may be understood with these free parameters;
26this is not possible with non-parametric PDFs.
27
28The RooParametricStepFunction has Nbins-1 free parameters. Note that
29the limits of the dependent variable must match the low and hi bin limits.
30
31Here is an example showing how to use the RooParametricStepFunction to fit toy
32data generated from a uniform distribution:
33
34~~~ {.cpp}
35// Define some constant parameters
36const int nBins = 10;
37const double xMin = 0.0;
38const double xMax = 10.0;
39const double binWidth = (xMax - xMin) / nBins;
40const std::size_t nEvents = 10000;
41
42// Fill the bin boundaries
43std::vector<double> binBoundaries(nBins + 1);
44
45for(int i = 0; i <= nBins; ++i) {
46 binBoundaries[i] = i * binWidth;
47}
48
49// The RooParametricStepFunction needs a TArrayD
50TArrayD binBoundariesTArr{int(binBoundaries.size()), binBoundaries.data()};
51
52RooRealVar x{"x", "x", xMin, xMax};
53
54// There are nBins-1 free coefficient parameters, whose sum must be <= 1.0.
55// We all set them to 0.1, such that the resulting step function pdf is
56// initially uniform.
57RooArgList coefList;
58for(int i = 0; i < nBins - 1; ++i) {
59 auto name = std::string("coef_") + std::to_string(i);
60 coefList.addOwned(std::make_unique<RooRealVar>(name.c_str(), name.c_str(), 0.1, 0.0, 1.0));
61}
62
63// Create the RooParametricStepFunction pdf
64RooParametricStepFunction pdf{"pdf", "pdf", x, coefList, binBoundariesTArr, int(nBins)};
65
66// Generate some toy data, following our uniform step function pdf
67std::unique_ptr<RooAbsData> data{pdf.generate(x, nEvents)};
68
69// Fit the step function to the toy data
70pdf.fitTo(*data);
71
72// Now we plot the data and the pdf, the latter should not be uniform
73// anymore because the coefficients were fit to the toy data
74RooPlot *xframe = x.frame(RooFit::Title("Fitting uniform toy data with a step function p.d.f."));
75data->plotOn(xframe);
76pdf.plotOn(xframe);
77xframe->Draw();
78~~~
79*/
80
82
83#include <RooCurve.h>
84#include <RooRealVar.h>
85
87
88////////////////////////////////////////////////////////////////////////////////
89/// Constructor
90
92 RooAbsReal& x, const RooArgList& coefList, TArrayD const& limits, Int_t nBins) :
93 RooAbsPdf(name, title),
94 _x("x", "Dependent", this, x),
95 _coefList("coefList","List of coefficients",this),
96 _nBins(nBins)
97{
98
99 // Check lowest order
100 if (_nBins<0) {
101 std::cout << "RooParametricStepFunction::ctor(" << GetName()
102 << ") WARNING: nBins must be >=0, setting value to 0" << std::endl ;
103 _nBins=0 ;
104 }
105
106 _coefList.addTyped<RooAbsReal>(coefList);
107
108 // Bin limits
109 limits.Copy(_limits);
110
111}
112
113////////////////////////////////////////////////////////////////////////////////
114/// Copy constructor
115
117 RooAbsPdf(other, name),
118 _x("x", this, other._x),
119 _coefList("coefList",this,other._coefList),
120 _nBins(other._nBins)
121{
122 (other._limits).Copy(_limits);
123}
124
125////////////////////////////////////////////////////////////////////////////////
126
127Int_t RooParametricStepFunction::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
128{
129 if (matchArgs(allVars, analVars, _x)) return 1;
130 return 0;
131}
132
133////////////////////////////////////////////////////////////////////////////////
134
135double RooParametricStepFunction::analyticalIntegral(Int_t /*code*/, const char* rangeName) const
136{
137 // Case without range is trivial: p.d.f is by construction normalized
138 if (!rangeName) {
139 return 1.0 ;
140 }
141
142 // Case with ranges, calculate integral explicitly
143 double xmin = _x.min(rangeName) ;
144 double xmax = _x.max(rangeName) ;
145
146 double sum=0 ;
147 Int_t i ;
148 for (i=1 ; i<=_nBins ; i++) {
149 double binVal = (i<_nBins) ? (static_cast<RooRealVar*>(_coefList.at(i-1))->getVal()) : lastBinValue() ;
150 if (_limits[i-1]>=xmin && _limits[i]<=xmax) {
151 // Bin fully in the integration domain
152 sum += (_limits[i]-_limits[i-1])*binVal ;
153 } else if (_limits[i-1]<xmin && _limits[i]>xmax) {
154 // Domain is fully contained in this bin
155 sum += (xmax-xmin)*binVal ;
156 // Exit here, this is the last bin to be processed by construction
157 return sum ;
158 } else if (_limits[i-1]<xmin && _limits[i]<=xmax && _limits[i]>xmin) {
159 // Lower domain boundary is in bin
160 sum += (_limits[i]-xmin)*binVal ;
161 } else if (_limits[i-1]>=xmin && _limits[i]>xmax && _limits[i-1]<xmax) {
162 sum += (xmax-_limits[i-1])*binVal ;
163 // Upper domain boundary is in bin
164 // Exit here, this is the last bin to be processed by construction
165 return sum ;
166 }
167 }
168
169 return sum;
170}
171
172////////////////////////////////////////////////////////////////////////////////
173
175{
176 double sum(0.);
177 double binSize(0.);
178 for (Int_t j=1;j<_nBins;j++){
179 RooRealVar* tmp = static_cast<RooRealVar*>(_coefList.at(j-1));
180 binSize = _limits[j] - _limits[j-1];
181 sum = sum + tmp->getVal()*binSize;
182 }
183 binSize = _limits[_nBins] - _limits[_nBins-1];
184 return (1.0 - sum)/binSize;
185}
186
187////////////////////////////////////////////////////////////////////////////////
188
190{
191 double value(0.);
192 if (_x >= _limits[0] && _x < _limits[_nBins]){
193
194 for (Int_t i=1;i<=_nBins;i++){
195 if (_x < _limits[i]){
196 // in Bin i-1 (starting with Bin 0)
197 if (i<_nBins) {
198 // not in last Bin
199 value = static_cast<RooRealVar*>(_coefList.at(i-1))->getVal();
200 break;
201 } else {
203 if (value<=0.0){
204 value = 0.000001;
205 // cout << "RooParametricStepFunction: sum of values gt 1.0 -- beware!!" <<endl;
206 }
207 break;
208 }
209 }
210 }
211
212 }
213 return value;
214
215}
216
217
218std::list<double>* RooParametricStepFunction::plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const
219{
220 if(obs.namePtr() != _x->namePtr()) {
221 return nullptr;
222 }
223
224 // Retrieve position of all bin boundaries
225 std::span<const double> boundaries{_limits.GetArray(), static_cast<std::size_t>(_limits.GetSize())};
226
227 // Use the helper function from RooCurve to make sure to get sampling hints
228 // that work with the RooFitPlotting.
229 return RooCurve::plotSamplingHintForBinBoundaries(boundaries, xlo, xhi);
230}
#define ClassImp(name)
Definition Rtypes.h:382
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
const TNamed * namePtr() const
De-duplicated pointer to this object's name.
Definition RooAbsArg.h:504
bool addTyped(const RooAbsCollection &list, bool silent=false)
Adds elements of a given RooAbsCollection to the container if they match the specified type.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
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().
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
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:24
static std::list< double > * plotSamplingHintForBinBoundaries(std::span< const double > boundaries, double xlo, double xhi)
Returns sampling hints for a histogram with given boundaries.
Definition RooCurve.cxx:893
The Parametric Step Function PDF is a binned distribution whose parameters are the heights of each bi...
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
std::list< double > * plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const override
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
Array of doubles (64 bits per element).
Definition TArrayD.h:27
void Copy(TArrayD &array) const
Definition TArrayD.h:42
const Double_t * GetArray() const
Definition TArrayD.h:43
Int_t GetSize() const
Definition TArray.h:47
void Copy(TObject &named) const override
Copy this to obj.
Definition TNamed.cxx:94
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Double_t x[n]
Definition legend1.C:17
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345