Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooHistConstraint.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 *
4 * Copyright (c) 2023, CERN
5 *
6 * Redistribution and use in source and binary forms,
7 * with or without modification, are permitted according to the terms
8 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
9 */
10
11/** \class RooHistConstraint
12 * \ingroup Roofit
13 * The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties.
14 * Following the Barlow-Beeston method, it adds Poisson constraints for each bin that
15 * constrain the statistical uncertainty of the template histogram.
16 *
17 * It can therefore be used to estimate the Monte Carlo uncertainty of a fit.
18 *
19 * Check also the tutorial rf709_BarlowBeeston.C
20 *
21 */
22
23#include <RooHistConstraint.h>
24
25#include <RooParamHistFunc.h>
26#include <RooRealVar.h>
27
29
31
32////////////////////////////////////////////////////////////////////////////////
33/// Create a new RooHistConstraint.
34/// \param[in] name Name of the PDF. This is used to identify it in a likelihood model.
35/// \param[in] title Title for plotting etc.
36/// \param[in] phfSet Set of parametrised histogram functions (RooParamHistFunc).
37/// \param[in] threshold Threshold (bin content) up to which statistical uncertainties are taken into account.
38RooHistConstraint::RooHistConstraint(const char *name, const char *title,
39 const RooArgSet& phfSet, int threshold) :
40 RooAbsPdf(name,title),
41 _gamma("gamma","gamma",this),
42 _nominal("nominal","nominal",this),
43 _relParam(true)
44{
45 // Implementing constraint on sum of RooParamHists
46 //
47 // Step 1 - Create new gamma parameters for sum
48 // Step 2 - Replace entries in gamma listproxy of components with new sum components
49 // Step 3 - Implement constraints in terms of gamma sum parameters
50
51
52 if (phfSet.size()==1) {
53
54 auto phf = dynamic_cast<RooParamHistFunc*>(phfSet.first()) ;
55
56 if (!phf) {
57 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
58 << ") ERROR: input object must be a RooParamHistFunc" << std::endl ;
59 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
60 }
61
62 // Now populate nominal with parameters
63 for (int i=0 ; i<phf->_dh.numEntries() ; i++) {
64 phf->_dh.get(i) ;
65 if (phf->_dh.weight()<threshold && phf->_dh.weight() != 0.) {
66 const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
67 auto var = std::make_unique<RooRealVar>(vname,vname,0,1.E30);
68 var->setVal(phf->_dh.weight()) ;
69 var->setConstant(true);
70
71 auto gamma = static_cast<RooRealVar*>(phf->_p.at(i));
72 if (var->getVal() > 0.0) {
73 gamma->setConstant(false);
74 }
75
76 _nominal.addOwned(std::move(var)) ;
77 _gamma.add(*gamma) ;
78 }
79 }
80
81 return ;
82 }
83
84
85
86 int nbins(-1) ;
87 std::vector<RooParamHistFunc*> phvec ;
88 RooArgSet gammaSet ;
89 std::string bin0_name ;
90 for (const auto arg : phfSet) {
91
92 auto phfComp = dynamic_cast<RooParamHistFunc*>(arg) ;
93 if (phfComp) {
94 phvec.push_back(phfComp) ;
95 if (nbins==-1) {
96 nbins = phfComp->_p.size() ;
97 bin0_name = phfComp->_p.at(0)->GetName() ;
98 gammaSet.add(phfComp->_p) ;
99 } else {
100 if (int(phfComp->_p.size())!=nbins) {
101 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
102 << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have same #bins" << std::endl ;
103 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
104 }
105 if (bin0_name != phfComp->_p.at(0)->GetName()) {
106 coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
107 << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have the same bin parameters.\n"
108 << "Previously found " << bin0_name << ", now found " << phfComp->_p.at(0)->GetName() << ".\n"
109 << "Check that the right RooParamHistFuncs have been passed to this RooHistConstraint." << std::endl;
110 throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
111 }
112
113 }
114 } else {
115 coutW(InputArguments) << "RooHistConstraint::ctor(" << GetName()
116 << ") WARNING: ignoring input argument " << arg->GetName() << " which is not of type RooParamHistFunc" << std::endl;
117 }
118 }
119
120 _gamma.add(gammaSet) ;
121
122 // Now populate nominal and nominalErr with parameters
123 for (int i=0 ; i<nbins ; i++) {
124
125 double sumVal(0) ;
126 for (const auto phfunc : phvec) {
127 sumVal += phfunc->getNominal(i);
128 }
129
130 if (sumVal<threshold && sumVal != 0.) {
131
132 const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
133 auto var = std::make_unique<RooRealVar>(vname,vname,0,1000);
134
135 double sumVal2(0) ;
136 for(auto const& elem : phvec) {
137 sumVal2 += elem->getNominal(i) ;
138 }
139 var->setVal(sumVal2) ;
140 var->setConstant(true) ;
141
142 vname = Form("%s_nominal_error_bin_%i",GetName(),i) ;
143 //RooRealVar* vare = new RooRealVar(vname,vname,0,1000) ;
144
145 //double sumErr2(0) ;
146 //for(auto const& elem : phvec) {
147 //sumErr2 += std::pow(elem->getNominalError(i),2) ;
148 //}
149 //vare->setVal(sqrt(sumErr2)) ;
150 //vare->setConstant(true) ;
151
152 _nominal.addOwned(std::move(var));
153 // _nominalErr.add(*vare) ;
154
155 (static_cast<RooRealVar*>(_gamma.at(i)))->setConstant(false) ;
156
157 }
158 }
159}
160
161////////////////////////////////////////////////////////////////////////////////
162
164 RooAbsPdf(other,name),
165 _gamma("gamma",this,other._gamma),
166 _nominal("nominal",this,other._nominal),
167 _relParam(other._relParam)
168 {
169 }
170
171////////////////////////////////////////////////////////////////////////////////
172
174 {
175 double prod(1.0);
176
177 for (unsigned int i=0; i < _nominal.size(); ++i) {
178 const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
179 const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
180 double gammaVal = gamma.getVal();
181 const int nomVal = static_cast<int>(nominal.getVal());
182
183 if (_relParam) {
184 gammaVal *= nomVal;
185 }
186
187 if (gammaVal>0) {
188 const double pois = ROOT::Math::poisson_pdf(nomVal, gammaVal);
189 prod *= pois;
190 } else if (nomVal > 0) {
191 coutE(Eval) << "ERROR in RooHistConstraint: gamma=0 and nom>0" << std::endl;
192 }
193 }
194
195 return prod;
196 }
197
198////////////////////////////////////////////////////////////////////////////////
199
200double RooHistConstraint::getLogVal(const RooArgSet* /*set*/) const
201{
202 double sum = 0.;
203 for (unsigned int i=0; i < _nominal.size(); ++i) {
204 const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
205 const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
206 double gammaVal = gamma.getVal();
207 const int nomVal = static_cast<int>(nominal.getVal());
208
209 if (_relParam) {
210 gammaVal *= nomVal;
211 }
212
213 if (gammaVal>0) {
214 const double logPoisson = nomVal * log(gammaVal) - gammaVal - std::lgamma(nomVal + 1);
215 sum += logPoisson ;
216 } else if (nomVal > 0) {
217 coutE(Eval) << "ERROR in RooHistConstraint: gamma=0 and nom>0" << std::endl;
218 }
219 }
220
221 return sum ;
222}
#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:2489
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
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.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
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