Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 <RooHistConstraint.h>
20
21#include <RooParamHistFunc.h>
22#include <RooRealVar.h>
23
25
27
28////////////////////////////////////////////////////////////////////////////////
29/// Create a new RooHistConstraint.
30/// \param[in] name Name of the PDF. This is used to identify it in a likelihood model.
31/// \param[in] title Title for plotting etc.
32/// \param[in] phfSet Set of parametrised histogram functions (RooParamHistFunc).
33/// \param[in] threshold Threshold (bin content) up to which statistcal uncertainties are taken into account.
34RooHistConstraint::RooHistConstraint(const char *name, const char *title,
35 const RooArgSet& phfSet, int threshold) :
36 RooAbsPdf(name,title),
37 _gamma("gamma","gamma",this),
38 _nominal("nominal","nominal",this),
39 _relParam(true)
40{
41 // Implementing constraint on sum of RooParamHists
42 //
43 // Step 1 - Create new gamma parameters for sum
44 // Step 2 - Replace entries in gamma listproxy of components with new sum components
45 // Step 3 - Implement constraints in terms of gamma sum parameters
46
47
48 if (phfSet.size()==1) {
49
50 auto phf = dynamic_cast<RooParamHistFunc*>(phfSet.first()) ;
51
52 if (!phf) {
53 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
54 << ") ERROR: input object must be a RooParamHistFunc" << std::endl ;
55 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
56 }
57
58 // Now populate nominal with parameters
59 for (int i=0 ; i<phf->_dh.numEntries() ; i++) {
60 phf->_dh.get(i) ;
61 if (phf->_dh.weight()<threshold && phf->_dh.weight() != 0.) {
62 const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
63 auto var = std::make_unique<RooRealVar>(vname,vname,0,1.E30);
64 var->setVal(phf->_dh.weight()) ;
65 var->setConstant(true);
66
67 auto gam = static_cast<RooRealVar*>(phf->_p.at(i));
68 if (var->getVal() > 0.0) {
69 gam->setConstant(false);
70 }
71
72 _nominal.addOwned(std::move(var)) ;
73 _gamma.add(*gam) ;
74 }
75 }
76
77 return ;
78 }
79
80
81
82 int nbins(-1) ;
83 std::vector<RooParamHistFunc*> phvec ;
84 RooArgSet gammaSet ;
85 std::string bin0_name ;
86 for (const auto arg : phfSet) {
87
88 auto phfComp = dynamic_cast<RooParamHistFunc*>(arg) ;
89 if (phfComp) {
90 phvec.push_back(phfComp) ;
91 if (nbins==-1) {
92 nbins = phfComp->_p.getSize() ;
93 bin0_name = phfComp->_p.at(0)->GetName() ;
94 gammaSet.add(phfComp->_p) ;
95 } else {
96 if (phfComp->_p.getSize()!=nbins) {
97 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
98 << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have same #bins" << std::endl ;
99 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
100 }
101 if (bin0_name != phfComp->_p.at(0)->GetName()) {
102 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
103 << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have the same bin parameters.\n"
104 << "Previously found " << bin0_name << ", now found " << phfComp->_p.at(0)->GetName() << ".\n"
105 << "Check that the right RooParamHistFuncs have been passed to this RooHistConstraint." << std::endl;
106 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
107 }
108
109 }
110 } else {
111 coutW(InputArguments) << "RooHistConstraint::ctor(" << GetName()
112 << ") WARNING: ignoring input argument " << arg->GetName() << " which is not of type RooParamHistFunc" << std::endl;
113 }
114 }
115
116 _gamma.add(gammaSet) ;
117
118 // Now populate nominal and nominalErr with parameters
119 for (int i=0 ; i<nbins ; i++) {
120
121 double sumVal(0) ;
122 for (const auto phfunc : phvec) {
123 sumVal += phfunc->getNominal(i);
124 }
125
126 if (sumVal<threshold && sumVal != 0.) {
127
128 const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
129 auto var = std::make_unique<RooRealVar>(vname,vname,0,1000);
130
131 double sumVal2(0) ;
132 for(auto const& elem : phvec) {
133 sumVal2 += elem->getNominal(i) ;
134 }
135 var->setVal(sumVal2) ;
136 var->setConstant(true) ;
137
138 vname = Form("%s_nominal_error_bin_%i",GetName(),i) ;
139 //RooRealVar* vare = new RooRealVar(vname,vname,0,1000) ;
140
141 //double sumErr2(0) ;
142 //for(auto const& elem : phvec) {
143 //sumErr2 += std::pow(elem->getNominalError(i),2) ;
144 //}
145 //vare->setVal(sqrt(sumErr2)) ;
146 //vare->setConstant(true) ;
147
148 _nominal.addOwned(std::move(var));
149 // _nominalErr.add(*vare) ;
150
151 ((RooRealVar*)_gamma.at(i))->setConstant(false) ;
152
153 }
154 }
155}
156
157////////////////////////////////////////////////////////////////////////////////
158
160 RooAbsPdf(other,name),
161 _gamma("gamma",this,other._gamma),
162 _nominal("nominal",this,other._nominal),
163 _relParam(other._relParam)
164 {
165 }
166
167////////////////////////////////////////////////////////////////////////////////
168
170 {
171 double prod(1.0);
172
173 for (unsigned int i=0; i < _nominal.size(); ++i) {
174 const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
175 const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
176 double gamVal = gamma.getVal();
177 const int nomVal = static_cast<int>(nominal.getVal());
178
179 if (_relParam) {
180 gamVal *= nomVal;
181 }
182
183 if (gamVal>0) {
184 const double pois = ROOT::Math::poisson_pdf(nomVal, gamVal);
185 prod *= pois;
186 } else if (nomVal > 0) {
187 coutE(Eval) << "ERROR in RooHistConstraint: gam=0 and nom>0" << std::endl;
188 }
189 }
190
191 return prod;
192 }
193
194////////////////////////////////////////////////////////////////////////////////
195
196double RooHistConstraint::getLogVal(const RooArgSet* /*set*/) const
197{
198 double sum = 0.;
199 for (unsigned int i=0; i < _nominal.size(); ++i) {
200 const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
201 const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
202 double gamVal = gamma.getVal();
203 const int nomVal = static_cast<int>(nominal.getVal());
204
205 if (_relParam) {
206 gamVal *= nomVal;
207 }
208
209 if (gamVal>0) {
210 const double logPoisson = nomVal * log(gamVal) - gamVal - std::lgamma(nomVal + 1);
211 sum += logPoisson ;
212 } else if (nomVal > 0) {
213 coutE(Eval) << "ERROR in RooHistConstraint: gam=0 and nom>0" << std::endl;
214 }
215 }
216
217 return sum ;
218}
#define coutW(a)
#define coutE(a)
#define ClassImp(name)
Definition Rtypes.h:377
char name[80]
Definition TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2467
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:62
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:55
bool addOwned(RooAbsArg &var, bool silent=false) override
Overloaded RooCollection_t::addOwned() method insert object into owning set and registers object as s...
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 RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties.
double getLogVal(const RooArgSet *set=nullptr) 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
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
double poisson_pdf(unsigned int n, double mu)
Probability density function of the Poisson distribution.
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345