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 
29 using 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.
39 RooHistConstraint::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(kTRUE)
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_t 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_t 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(kTRUE) ;
146 
147  vname = Form("%s_nominal_error_bin_%i",GetName(),i) ;
148  RooRealVar* vare = new RooRealVar(vname,vname,0,1000) ;
149 
150  Double_t 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(kTRUE) ;
156 
157  allVars.add(RooArgSet(*var,*vare)) ;
158  _nominal.add(*var) ;
159  // _nominalErr.add(*vare) ;
160 
161  ((RooRealVar*)_gamma.at(i))->setConstant(kFALSE) ;
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 
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 
ROOT::Math::Cephes::gamma
double gamma(double x)
Definition: SpecFuncCephes.cxx:339
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:226
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooAbsReal.h
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:176
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:68
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
Form
char * Form(const char *fmt,...)
TMath::Poisson
Double_t Poisson(Double_t x, Double_t par)
Compute the Poisson distribution function for (x,par).
Definition: TMath.cxx:564
coutE
#define coutE(a)
Definition: RooMsgService.h:33
coutW
#define coutW(a)
Definition: RooMsgService.h:32
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooDataHist::get
virtual const RooArgSet * get() const
Definition: RooDataHist.h:78
log
double log(double)
Int_t
int Int_t
Definition: RtypesCore.h:45
RooParamHistFunc::_dh
RooDataHist _dh
Definition: RooParamHistFunc.h:56
RooDataHist::weight
virtual Double_t weight() const
Definition: RooDataHist.h:105
RooHistConstraint::getLogVal
Double_t getLogVal(const RooArgSet *set=0) const
Return the log of the current value with given normalization An error message is printed if the argum...
Definition: RooHistConstraint.cxx:201
RooArgSet::add
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooAbsReal
Definition: RooAbsReal.h:61
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:88
RooParamHistFunc::_p
RooListProxy _p
Definition: RooParamHistFunc.h:55
RooHistConstraint::evaluate
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooHistConstraint.cxx:180
RooHistConstraint::_relParam
Bool_t _relParam
Definition: RooHistConstraint.h:33
RooParamHistFunc.h
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooHistConstraint::_gamma
RooListProxy _gamma
Definition: RooHistConstraint.h:31
RooAbsArg::RooArgSet
friend class RooArgSet
Definition: RooAbsArg.h:579
RooAbsCollection::size
Storage_t::size_type size() const
Definition: RooAbsCollection.h:154
RooListProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
Definition: RooListProxy.cxx:104
sqrt
double sqrt(double)
RooRealVar.h
RooHistConstraint
Definition: RooHistConstraint.h:19
sum
static long int sum(long int i)
Definition: Factory.cxx:2272
Double_t
double Double_t
Definition: RtypesCore.h:59
RooAbsCategory.h
RooAbsArg::addOwnedComponents
Bool_t addOwnedComponents(const RooArgSet &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2109
RooHistConstraint.h
RooDataHist::numEntries
virtual Int_t numEntries() const
Return the number of bins.
Definition: RooDataHist.cxx:1694
RooParamHistFunc
Definition: RooParamHistFunc.h:20
name
char name[80]
Definition: TGX11.cxx:110
RooAbsRealLValue::setConstant
void setConstant(Bool_t value=kTRUE)
Definition: RooAbsRealLValue.h:109
RooHistConstraint::RooHistConstraint
RooHistConstraint()
Definition: RooHistConstraint.h:21
ROOT::Math::lgamma
double lgamma(double x)
Calculates the logarithm of the gamma function.
Definition: SpecFuncMathCore.cxx:74
RooAbsPdf
Definition: RooAbsPdf.h:40
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooRealVar
Definition: RooRealVar.h:35
RooHistConstraint::_nominal
RooListProxy _nominal
Definition: RooHistConstraint.h:32
Riostream.h
pow
double pow(double, double)
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
TMath.h
RooArgSet
Definition: RooArgSet.h:28
int