Logo ROOT  
Reference Guide
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
81#include "Riostream.h"
82#include "TArrayD.h"
83#include <math.h>
84
86#include "RooAbsReal.h"
87#include "RooRealVar.h"
88#include "RooArgList.h"
89
90#include "TError.h"
91
93
94////////////////////////////////////////////////////////////////////////////////
95/// Constructor
96
98 RooAbsReal& x, const RooArgList& coefList, TArrayD& limits, Int_t nBins) :
99 RooAbsPdf(name, title),
100 _x("x", "Dependent", this, x),
101 _coefList("coefList","List of coefficients",this),
102 _nBins(nBins)
103{
104
105 // Check lowest order
106 if (_nBins<0) {
107 std::cout << "RooParametricStepFunction::ctor(" << GetName()
108 << ") WARNING: nBins must be >=0, setting value to 0" << std::endl ;
109 _nBins=0 ;
110 }
111
112 for (auto *coef : coefList) {
113 if (!dynamic_cast<RooAbsReal*>(coef)) {
114 std::cout << "RooParametricStepFunction::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
115 << " is not of type RooAbsReal" << std::endl ;
116 R__ASSERT(0) ;
117 }
118 _coefList.add(*coef) ;
119 }
120
121 // Bin limits
122 limits.Copy(_limits);
123
124}
125
126////////////////////////////////////////////////////////////////////////////////
127/// Copy constructor
128
130 RooAbsPdf(other, name),
131 _x("x", this, other._x),
132 _coefList("coefList",this,other._coefList),
133 _nBins(other._nBins)
134{
135 (other._limits).Copy(_limits);
136}
137
138////////////////////////////////////////////////////////////////////////////////
139
140Int_t RooParametricStepFunction::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
141{
142 if (matchArgs(allVars, analVars, _x)) return 1;
143 return 0;
144}
145
146////////////////////////////////////////////////////////////////////////////////
147
148double RooParametricStepFunction::analyticalIntegral(Int_t code, const char* rangeName) const
149{
150 R__ASSERT(code==1) ;
151
152 // Case without range is trivial: p.d.f is by construction normalized
153 if (!rangeName) {
154 return 1.0 ;
155 }
156
157 // Case with ranges, calculate integral explicitly
158 double xmin = _x.min(rangeName) ;
159 double xmax = _x.max(rangeName) ;
160
161 double sum=0 ;
162 Int_t i ;
163 for (i=1 ; i<=_nBins ; i++) {
164 double binVal = (i<_nBins) ? (static_cast<RooRealVar*>(_coefList.at(i-1))->getVal()) : lastBinValue() ;
165 if (_limits[i-1]>=xmin && _limits[i]<=xmax) {
166 // Bin fully in the integration domain
167 sum += (_limits[i]-_limits[i-1])*binVal ;
168 } else if (_limits[i-1]<xmin && _limits[i]>xmax) {
169 // Domain is fully contained in this bin
170 sum += (xmax-xmin)*binVal ;
171 // Exit here, this is the last bin to be processed by construction
172 return sum ;
173 } else if (_limits[i-1]<xmin && _limits[i]<=xmax && _limits[i]>xmin) {
174 // Lower domain boundary is in bin
175 sum += (_limits[i]-xmin)*binVal ;
176 } else if (_limits[i-1]>=xmin && _limits[i]>xmax && _limits[i-1]<xmax) {
177 sum += (xmax-_limits[i-1])*binVal ;
178 // Upper domain boundary is in bin
179 // Exit here, this is the last bin to be processed by construction
180 return sum ;
181 }
182 }
183
184 return sum;
185}
186
187////////////////////////////////////////////////////////////////////////////////
188
190{
191 double sum(0.);
192 double binSize(0.);
193 for (Int_t j=1;j<_nBins;j++){
194 RooRealVar* tmp = (RooRealVar*) _coefList.at(j-1);
195 binSize = _limits[j] - _limits[j-1];
196 sum = sum + tmp->getVal()*binSize;
197 }
198 binSize = _limits[_nBins] - _limits[_nBins-1];
199 return (1.0 - sum)/binSize;
200}
201
202////////////////////////////////////////////////////////////////////////////////
203
205{
206 double value(0.);
207 if (_x >= _limits[0] && _x < _limits[_nBins]){
208
209 for (Int_t i=1;i<=_nBins;i++){
210 if (_x < _limits[i]){
211 // in Bin i-1 (starting with Bin 0)
212 if (i<_nBins) {
213 // not in last Bin
214 RooRealVar* tmp = (RooRealVar*) _coefList.at(i-1);
215 value = tmp->getVal();
216 break;
217 } else {
218 // in last Bin
219 double sum(0.);
220 double binSize(0.);
221 for (Int_t j=1;j<_nBins;j++){
222 RooRealVar* tmp = (RooRealVar*) _coefList.at(j-1);
223 binSize = _limits[j] - _limits[j-1];
224 sum = sum + tmp->getVal()*binSize;
225 }
226 binSize = _limits[_nBins] - _limits[_nBins-1];
227 value = (1.0 - sum)/binSize;
228 if (value<=0.0){
229 value = 0.000001;
230 // cout << "RooParametricStepFunction: sum of values gt 1.0 -- beware!!" <<endl;
231 }
232 break;
233 }
234 }
235 }
236
237 }
238 return value;
239
240}
#define ClassImp(name)
Definition: Rtypes.h:375
#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 value
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:94
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:57
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...
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.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
double max(const char *rname=0) const
Query upper 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
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