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