Logo ROOT   master
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 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
static long int sum(long int i)
Definition: Factory.cxx:2274
#define coutE(a)
Definition: RooMsgService.h:33
Storage_t::size_type size() const
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
RooListProxy _nominal
virtual const RooArgSet * get() const
Definition: RooDataHist.h:79
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties...
int Int_t
Definition: RtypesCore.h:41
Bool_t addOwnedComponents(const RooArgSet &comps)
Take ownership of the contents of &#39;comps&#39;.
Definition: RooAbsArg.cxx:2112
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:32
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...
double sqrt(double)
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
double pow(double, double)
virtual Double_t weight() const
Definition: RooDataHist.h:106
friend class RooArgSet
Definition: RooAbsArg.h:550
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:232
double gamma(double x)
Int_t getSize() const
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
RooAbsArg * first() const
void setConstant(Bool_t value=kTRUE)
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:89
char * Form(const char *fmt,...)
A histogram function that assigns scale parameters to every bin.
Double_t Poisson(Double_t x, Double_t par)
const Bool_t kFALSE
Definition: RtypesCore.h:88
#define ClassImp(name)
Definition: Rtypes.h:365
virtual Int_t numEntries() const
Return the number of bins.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
double lgamma(double x)
Calculates the logarithm of the gamma function.
const Bool_t kTRUE
Definition: RtypesCore.h:87
char name[80]
Definition: TGX11.cxx:109
double log(double)