Logo ROOT  
Reference Guide
RooHistConstraint.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * *
4 * This code was autogenerated by RooClassFactory *
5 *****************************************************************************/
6
7/** \class RooHistConstraint
8 * \ingroup Roofit
9 * The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties.
10 * Following the Barlow-Beeston method, it adds Poisson constraints for each bin that
11 * constrain the statistical uncertainty of the template histogram.
12 *
13 * It can therefore be used to estimate the Monte Carlo uncertainty of a fit.
14 *
15 * Check also the tutorial rf709_BarlowBeeston.C
16 *
17 */
18
19#include "Riostream.h"
20
21#include "RooHistConstraint.h"
22#include "RooAbsReal.h"
23#include "RooAbsCategory.h"
24#include "RooParamHistFunc.h"
25#include "RooRealVar.h"
26#include <math.h>
27#include "TMath.h"
28
29using namespace std;
30
32
33////////////////////////////////////////////////////////////////////////////////
34/// Create a new RooHistConstraint.
35/// \param[in] name Name of the PDF. This is used to identify it in a likelihood model.
36/// \param[in] title Title for plotting etc.
37/// \param[in] phfSet Set of parametrised histogram functions (RooParamHistFunc).
38/// \param[in] threshold Threshold (bin content) up to which statistcal uncertainties are taken into account.
39RooHistConstraint::RooHistConstraint(const char *name, const char *title,
40 const RooArgSet& phfSet, Int_t threshold) :
41 RooAbsPdf(name,title),
42 _gamma("gamma","gamma",this),
43 _nominal("nominal","nominal",this),
44 _relParam(true)
45{
46 // Implementing constraint on sum of RooParamHists
47 //
48 // Step 1 - Create new gamma parameters for sum
49 // Step 2 - Replace entries in gamma listproxy of components with new sum components
50 // Step 3 - Implement constraints in terms of gamma sum parameters
51
52
53 if (phfSet.getSize()==1) {
54
55 RooParamHistFunc* phf = dynamic_cast<RooParamHistFunc*>(phfSet.first()) ;
56
57 if (!phf) {
58 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
59 << ") ERROR: input object must be a RooParamHistFunc" << endl ;
60 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
61 }
62
63 // Now populate nominal with parameters
64 RooArgSet allVars ;
65 for (Int_t i=0 ; i<phf->_dh.numEntries() ; i++) {
66 phf->_dh.get(i) ;
67 if (phf->_dh.weight()<threshold && phf->_dh.weight() != 0.) {
68 const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
69 RooRealVar* var = new RooRealVar(vname,vname,0,1.E30) ;
70 var->setVal(phf->_dh.weight()) ;
71 var->setConstant(true);
72 allVars.add(*var) ;
73 _nominal.add(*var) ;
74
75 RooRealVar* gam = (RooRealVar*) phf->_p.at(i) ;
76 if (var->getVal()>0) {
77 gam->setConstant(false);
78 }
79 _gamma.add(*gam) ;
80 }
81 }
82
83 addOwnedComponents(allVars) ;
84
85 return ;
86 }
87
88
89
90 Int_t nbins(-1) ;
91 vector<RooParamHistFunc*> phvec ;
92 RooArgSet gammaSet ;
93 string bin0_name ;
94 for (const auto arg : phfSet) {
95
96 RooParamHistFunc* phfComp = dynamic_cast<RooParamHistFunc*>(arg) ;
97 if (phfComp) {
98 phvec.push_back(phfComp) ;
99 if (nbins==-1) {
100 nbins = phfComp->_p.getSize() ;
101 bin0_name = phfComp->_p.at(0)->GetName() ;
102 gammaSet.add(phfComp->_p) ;
103 } else {
104 if (phfComp->_p.getSize()!=nbins) {
105 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
106 << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have same #bins" << endl ;
107 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
108 }
109 if (bin0_name != phfComp->_p.at(0)->GetName()) {
110 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
111 << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have the same bin parameters.\n"
112 << "Previously found " << bin0_name << ", now found " << phfComp->_p.at(0)->GetName() << ".\n"
113 << "Check that the right RooParamHistFuncs have been passed to this RooHistConstraint." << std::endl;
114 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
115 }
116
117 }
118 } else {
119 coutW(InputArguments) << "RooHistConstraint::ctor(" << GetName()
120 << ") WARNING: ignoring input argument " << arg->GetName() << " which is not of type RooParamHistFunc" << endl;
121 }
122 }
123
124 _gamma.add(gammaSet) ;
125
126 // Now populate nominal and nominalErr with parameters
127 RooArgSet allVars ;
128 for (Int_t i=0 ; i<nbins ; i++) {
129
130 double sumVal(0) ;
131 for (const auto phfunc : phvec) {
132 sumVal += phfunc->getNominal(i);
133 }
134
135 if (sumVal<threshold && sumVal != 0.) {
136
137 const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
138 RooRealVar* var = new RooRealVar(vname,vname,0,1000) ;
139
140 double sumVal2(0) ;
141 for (vector<RooParamHistFunc*>::iterator iter = phvec.begin() ; iter != phvec.end() ; ++iter) {
142 sumVal2 += (*iter)->getNominal(i) ;
143 }
144 var->setVal(sumVal2) ;
145 var->setConstant(true) ;
146
147 vname = Form("%s_nominal_error_bin_%i",GetName(),i) ;
148 RooRealVar* vare = new RooRealVar(vname,vname,0,1000) ;
149
150 double sumErr2(0) ;
151 for (vector<RooParamHistFunc*>::iterator iter = phvec.begin() ; iter != phvec.end() ; ++iter) {
152 sumErr2 += pow((*iter)->getNominalError(i),2) ;
153 }
154 vare->setVal(sqrt(sumErr2)) ;
155 vare->setConstant(true) ;
156
157 allVars.add(RooArgSet(*var,*vare)) ;
158 _nominal.add(*var) ;
159 // _nominalErr.add(*vare) ;
160
161 ((RooRealVar*)_gamma.at(i))->setConstant(false) ;
162
163 }
164 }
165 addOwnedComponents(allVars) ;
166}
167
168////////////////////////////////////////////////////////////////////////////////
169
171 RooAbsPdf(other,name),
172 _gamma("gamma",this,other._gamma),
173 _nominal("nominal",this,other._nominal),
174 _relParam(other._relParam)
175 {
176 }
177
178////////////////////////////////////////////////////////////////////////////////
179
181 {
182 double prod(1);
183
184 for (unsigned int i=0; i < _nominal.size(); ++i) {
185 const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
186 const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
187
188 double gamVal = gamma.getVal();
189 if (_relParam)
190 gamVal *= nominal.getVal();
191
192 const double pois = TMath::Poisson(nominal.getVal(),gamVal);
193 prod *= pois;
194 }
195
196 return prod;
197 }
198
199////////////////////////////////////////////////////////////////////////////////
200
201double RooHistConstraint::getLogVal(const RooArgSet* /*set*/) const
202{
203 double sum = 0.;
204 for (unsigned int i=0; i < _nominal.size(); ++i) {
205 const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
206 const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
207 double gam = gamma.getVal();
208 Int_t nom = (Int_t) nominal.getVal();
209
210 if (_relParam)
211 gam *= nom;
212
213 if (gam>0) {
214 const double logPoisson = nom * log(gam) - gam - std::lgamma(nom+1);
215 sum += logPoisson ;
216 } else if (nom>0) {
217 cerr << "ERROR in RooHistConstraint: gam=0 and nom>0" << endl ;
218 }
219 }
220
221 return sum ;
222}
223
224
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
int Int_t
Definition: RtypesCore.h:45
#define ClassImp(name)
Definition: Rtypes.h:375
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2447
friend class RooArgSet
Definition: RooAbsArg.h:645
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2185
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
void setConstant(bool value=true)
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
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...
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:110
Int_t numEntries() const override
Return the number of bins.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:82
The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties.
RooListProxy _nominal
double getLogVal(const RooArgSet *set=0) const override
Return the log of the current value with given normalization An error message is printed if the argum...
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
A histogram function that assigns scale parameters to every bin.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:281
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
double lgamma(double x)
Calculates the logarithm of the gamma function.
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1753
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1748
double gamma(double x)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
@ InputArguments
Definition: RooGlobalFunc.h:64
Double_t Poisson(Double_t x, Double_t par)
Computes the Poisson distribution function for (x,par).
Definition: TMath.cxx:569
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345