diff --git a/roofit/histfactory/src/ParamHistFunc.cxx b/roofit/histfactory/src/ParamHistFunc.cxx
index e492cf58982848b593ad75de04f445224dd7f811..19ec94514dc2693438ca6a995078c7b93180035a 100644
--- a/roofit/histfactory/src/ParamHistFunc.cxx
+++ b/roofit/histfactory/src/ParamHistFunc.cxx
@@ -124,8 +124,8 @@ ParamHistFunc::ParamHistFunc(const char* name, const char* title,
////////////////////////////////////////////////////////////////////////////////
-/// Create a function which returns binewise-values
-/// This class contains N RooRealVar's, one for each
+/// Create a function which returns bin-wise values.
+/// This class contains N RooRealVars, one for each
/// bin from the given RooRealVar.
///
/// The value of the function in the ith bin is
@@ -133,8 +133,7 @@ ParamHistFunc::ParamHistFunc(const char* name, const char* title,
///
/// F(i) = gamma_i * nominal(i)
///
-/// Where the nominal values are simply fixed
-/// numbers (default = 1.0 for all i)
+/// Where the nominal values are taken from the histogram.
ParamHistFunc::ParamHistFunc(const char* name, const char* title,
const RooArgList& vars, const RooArgList& paramSet,
const TH1* Hist ) :
@@ -214,7 +213,7 @@ ParamHistFunc::~ParamHistFunc()
////////////////////////////////////////////////////////////////////////////////
/// Get the index of the gamma parameter associated
-/// with the current bin
+/// with the current bin.
/// This number is the "RooDataSet" style index
/// and it must be because it uses the RooDataSet method directly
/// This is intended to be fed into the getParameter(Int_t) method:
@@ -460,9 +459,7 @@ RooArgList ParamHistFunc::createParamSet(RooWorkspace& w, const std::string& Pre
RooArgList params = ParamHistFunc::createParamSet( w, Prefix, vars );
- RooFIter paramIter = params.fwdIterator() ;
- RooAbsArg* comp ;
- while((comp = (RooAbsArg*) paramIter.next())) {
+ for (auto comp : params) {
RooRealVar* var = (RooRealVar*) comp;
@@ -493,7 +490,7 @@ RooArgList ParamHistFunc::createParamSet(const std::string& Prefix, Int_t numBin
if( gamma_max <= gamma_min ) {
- std::cout << "Warming: gamma_min <= gamma_max: Using default values (0, 10)" << std::endl;
+ std::cout << "Warning: gamma_min <= gamma_max: Using default values (0, 10)" << std::endl;
gamma_min = 0.0;
gamma_max = 10.0;
diff --git a/roofit/roofit/inc/RooHistConstraint.h b/roofit/roofit/inc/RooHistConstraint.h
index c199bc68c36823bec5c65f6b821e777ee46bde68..b0ef2ede05e566ae83f2b73982a728311dfc602f 100644
--- a/roofit/roofit/inc/RooHistConstraint.h
+++ b/roofit/roofit/inc/RooHistConstraint.h
@@ -28,18 +28,15 @@ public:
protected:
- Double_t logSum(Int_t i) const ;
-
RooListProxy _gamma ;
RooListProxy _nominal ;
- RooListProxy _nominalErr ;
Bool_t _relParam ;
Double_t evaluate() const ;
private:
- ClassDef(RooHistConstraint,1) // Your description goes here...
+ ClassDef(RooHistConstraint, 2)
};
#endif
diff --git a/roofit/roofit/src/RooHistConstraint.cxx b/roofit/roofit/src/RooHistConstraint.cxx
index 8d1e16123f18ea457dc56aaa79baed3db7dcbd7e..03f26eb4fc6b0989b58a05fcaeafd385774b6bb7 100644
--- a/roofit/roofit/src/RooHistConstraint.cxx
+++ b/roofit/roofit/src/RooHistConstraint.cxx
@@ -5,9 +5,16 @@
*****************************************************************************/
/** \class RooHistConstraint
- \ingroup Roofit
-
-**/
+ * \ingroup Roofit
+ * The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties.
+ * Following the Barlow-Beeston method, it adds Poisson constraints for each bin that
+ * constrain the statistical uncertainty of the template histogram.
+ *
+ * It can therefore be used to estimate the Monte Carlo uncertainty of a fit.
+ *
+ * Check also the tutorial rf709_BarlowBeeston.C
+ *
+ */
#include "Riostream.h"
@@ -24,12 +31,16 @@ using namespace std;
ClassImp(RooHistConstraint);
////////////////////////////////////////////////////////////////////////////////
-
-RooHistConstraint::RooHistConstraint(const char *name, const char *title, const RooArgSet& phfSet, Int_t threshold) :
+/// Create a new RooHistConstraint.
+/// \param[in] name Name of the PDF. This is used to identify it in a likelihood model.
+/// \param[in] title Title for plotting etc.
+/// \param[in] phfSet Set of parametrised histogram functions (RooParamHistFunc).
+/// \param[in] threshold Threshold (bin content) up to which statistcal uncertainties are taken into account.
+RooHistConstraint::RooHistConstraint(const char *name, const char *title,
+ const RooArgSet& phfSet, Int_t threshold) :
RooAbsPdf(name,title),
_gamma("gamma","gamma",this),
_nominal("nominal","nominal",this),
- _nominalErr("nominalErr","nominalErr",this),
_relParam(kTRUE)
{
// Implementing constraint on sum of RooParamHists
@@ -45,28 +56,30 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title, const
if (!phf) {
coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
- << ") ERROR: input object must be a RooParamHistFunc" << endl ;
- throw std::string("RoohistConstraint::ctor ERROR incongruent input arguments") ;
+ << ") ERROR: input object must be a RooParamHistFunc" << endl ;
+ throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
}
// Now populate nominal with parameters
RooArgSet allVars ;
for (Int_t i=0 ; i<phf->_dh.numEntries() ; i++) {
phf->_dh.get(i) ;
- if (phf->_dh.weight()<threshold) {
- const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
- RooRealVar* var = new RooRealVar(vname,vname,0,1000) ;
- var->setVal(phf->_dh.weight()) ;
- var->setConstant(kTRUE) ;
- allVars.add(*var) ;
- _nominal.add(*var) ;
- RooRealVar* gam = (RooRealVar*) phf->_p.at(i) ;
- if (var->getVal()>0) {
- gam->setConstant(kFALSE) ;
- }
- _gamma.add(*gam) ;
+ if (phf->_dh.weight()<threshold && phf->_dh.weight() != 0.) {
+ const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
+ RooRealVar* var = new RooRealVar(vname,vname,0,1.E30) ;
+ var->setVal(phf->_dh.weight()) ;
+ var->setConstant(true);
+ allVars.add(*var) ;
+ _nominal.add(*var) ;
+
+ RooRealVar* gam = (RooRealVar*) phf->_p.at(i) ;
+ if (var->getVal()>0) {
+ gam->setConstant(false);
+ }
+ _gamma.add(*gam) ;
}
}
+
addOwnedComponents(allVars) ;
return ;
@@ -74,37 +87,37 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title, const
- RooFIter phiter = phfSet.fwdIterator() ;
- RooAbsArg* arg ;
Int_t nbins(-1) ;
vector<RooParamHistFunc*> phvec ;
RooArgSet gammaSet ;
string bin0_name ;
- while((arg=phiter.next())) {
+ for (const auto arg : phfSet) {
RooParamHistFunc* phfComp = dynamic_cast<RooParamHistFunc*>(arg) ;
if (phfComp) {
phvec.push_back(phfComp) ;
if (nbins==-1) {
- nbins = phfComp->_p.getSize() ;
- bin0_name = phfComp->_p.at(0)->GetName() ;
- gammaSet.add(phfComp->_p) ;
+ nbins = phfComp->_p.getSize() ;
+ bin0_name = phfComp->_p.at(0)->GetName() ;
+ gammaSet.add(phfComp->_p) ;
} else {
- if (phfComp->_p.getSize()!=nbins) {
- coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
- << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have same #bins" << endl ;
- throw std::string("RoohistConstraint::ctor ERROR incongruent input arguments") ;
- }
- if (bin0_name != phfComp->_p.at(0)->GetName()) {
- coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
- << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have the same bin parameters" << endl ;
- throw std::string("RoohistConstraint::ctor ERROR incongruent input arguments") ;
- }
+ if (phfComp->_p.getSize()!=nbins) {
+ coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
+ << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have same #bins" << endl ;
+ throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
+ }
+ if (bin0_name != phfComp->_p.at(0)->GetName()) {
+ coutE(InputArguments) << "RooHistConstraint::ctor(" << GetName()
+ << ") ERROR: incongruent input arguments: all input RooParamHistFuncs should have the same bin parameters.\n"
+ << "Previously found " << bin0_name << ", now found " << phfComp->_p.at(0)->GetName() << ".\n"
+ << "Check that the right RooParamHistFuncs have been passed to this RooHistConstraint." << std::endl;
+ throw std::string("RooHistConstraint::ctor ERROR incongruent input arguments") ;
+ }
}
} else {
coutW(InputArguments) << "RooHistConstraint::ctor(" << GetName()
- << ") WARNING: ignoring input argument " << arg->GetName() << " which is not of type RooParamHistFunc" << endl;
+ << ") WARNING: ignoring input argument " << arg->GetName() << " which is not of type RooParamHistFunc" << endl;
}
}
@@ -115,18 +128,18 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title, const
for (Int_t i=0 ; i<nbins ; i++) {
Double_t sumVal(0) ;
- for (vector<RooParamHistFunc*>::iterator iter = phvec.begin() ; iter != phvec.end() ; ++iter) {
- sumVal += (*iter)->getNominal(i) ;
+ for (const auto phfunc : phvec) {
+ sumVal += phfunc->getNominal(i);
}
- if (sumVal<threshold) {
+ if (sumVal<threshold && sumVal != 0.) {
const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
RooRealVar* var = new RooRealVar(vname,vname,0,1000) ;
Double_t sumVal2(0) ;
for (vector<RooParamHistFunc*>::iterator iter = phvec.begin() ; iter != phvec.end() ; ++iter) {
- sumVal2 += (*iter)->getNominal(i) ;
+ sumVal2 += (*iter)->getNominal(i) ;
}
var->setVal(sumVal2) ;
var->setConstant(kTRUE) ;
@@ -136,14 +149,14 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title, const
Double_t sumErr2(0) ;
for (vector<RooParamHistFunc*>::iterator iter = phvec.begin() ; iter != phvec.end() ; ++iter) {
- sumErr2 += pow((*iter)->getNominalError(i),2) ;
+ sumErr2 += pow((*iter)->getNominalError(i),2) ;
}
vare->setVal(sqrt(sumErr2)) ;
vare->setConstant(kTRUE) ;
allVars.add(RooArgSet(*var,*vare)) ;
_nominal.add(*var) ;
- _nominalErr.add(*vare) ;
+ // _nominalErr.add(*vare) ;
((RooRealVar*)_gamma.at(i))->setConstant(kFALSE) ;
@@ -158,7 +171,6 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title, const
RooAbsPdf(other,name),
_gamma("gamma",this,other._gamma),
_nominal("nominal",this,other._nominal),
- _nominalErr("nominalErr",this,other._nominalErr),
_relParam(other._relParam)
{
}
@@ -167,69 +179,46 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title, const
Double_t RooHistConstraint::evaluate() const
{
- Double_t prod(1) ;
- RooFIter niter = _nominal.fwdIterator() ;
- RooFIter giter = _gamma.fwdIterator() ;
- RooAbsReal *gam, *nom ;
- while ((gam = (RooAbsReal*) giter.next())) {
- nom = (RooAbsReal*) niter.next() ;
- Double_t gamVal = gam->getVal() ;
- if (_relParam) gamVal *= nom->getVal() ;
- Double_t pois = TMath::Poisson(nom->getVal(),gamVal) ;
- prod *= pois ;
+ double prod(1);
+
+ for (unsigned int i=0; i < _nominal.size(); ++i) {
+ const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
+ const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
+
+ double gamVal = gamma.getVal();
+ if (_relParam)
+ gamVal *= nominal.getVal();
+
+ const double pois = TMath::Poisson(nominal.getVal(),gamVal);
+ prod *= pois;
}
- return prod ;
+
+ return prod;
}
////////////////////////////////////////////////////////////////////////////////
Double_t RooHistConstraint::getLogVal(const RooArgSet* /*set*/) const
{
- Double_t sum(0) ;
- RooFIter niter = _nominal.fwdIterator() ;
- RooFIter giter = _gamma.fwdIterator() ;
- RooAbsReal *gamv, *nomv ;
- while ((gamv = (RooAbsReal*) giter.next())) {
- nomv = (RooAbsReal*) niter.next() ;
- Double_t gam = gamv->getVal() ;
- Int_t nom = (Int_t) nomv->getVal() ;
- if (_relParam) gam *= nom ;
+ double sum = 0.;
+ for (unsigned int i=0; i < _nominal.size(); ++i) {
+ const auto& gamma = static_cast<const RooAbsReal&>(_gamma[i]);
+ const auto& nominal = static_cast<const RooAbsReal&>(_nominal[i]);
+ double gam = gamma.getVal();
+ Int_t nom = (Int_t) nominal.getVal();
+
+ if (_relParam)
+ gam *= nom;
+
if (gam>0) {
- Double_t logPoisson = nom*log(gam) - gam - logSum(nom) ;
+ const double logPoisson = nom * log(gam) - gam - std::lgamma(nom+1);
sum += logPoisson ;
} else if (nom>0) {
- cout << "ERROR gam=0 and nom>0" << endl ;
+ cerr << "ERROR in RooHistConstraint: gam=0 and nom>0" << endl ;
}
}
+
return sum ;
}
-////////////////////////////////////////////////////////////////////////////////
-Double_t RooHistConstraint::logSum(Int_t i) const
-{
- static Double_t* _lut = 0 ;
- if (!_lut) {
- _lut = new Double_t[5000] ;
- for (Int_t ii=0 ; ii<5000 ; ii++) _lut[ii] = 0 ;
-
- for (Int_t j=1 ; j<=5000 ; j++) {
- Double_t logj = log((Double_t)j) ;
- for (Int_t ii=j ; ii<=5000 ; ii++) {
- _lut[ii-1] += logj ;
- }
- }
- }
-
- if (i<5000) {
- return _lut[i-1] ;
- } else {
- Double_t ret = _lut[4999] ;
- cout << "logSum i=" << i << endl ;
- for (Int_t j=5000 ; j<=i ; j++) {
- ret += log((Double_t)j) ;
- }
- return ret ;
- }
-
-}
diff --git a/roofit/roofit/src/RooParamHistFunc.cxx b/roofit/roofit/src/RooParamHistFunc.cxx
index f3dca0100ff96883b68c59656fd3f20e66c72251..407d30d2bb4036d40c2f914f2e3f986a05d29a59 100644
--- a/roofit/roofit/src/RooParamHistFunc.cxx
+++ b/roofit/roofit/src/RooParamHistFunc.cxx
@@ -5,9 +5,19 @@
*****************************************************************************/
/** \class RooParamHistFunc
- \ingroup Roofit
-
-*/
+ * \ingroup Roofit
+ * A histogram function that assigns scale parameters to every bin. Instead of the bare bin contents,
+ * it therefore yields:
+ * \f[
+ * \gamma_{i} * \mathrm{bin}_i
+ * \f]
+ *
+ * The \f$ \gamma_i \f$ can therefore be used to parametrise statistical uncertainties of the histogram
+ * template. In conjuction with a constraint term, this can be used to implement the Barlow-Beeston method.
+ * The constraint can be implemented using RooHistConstraint.
+ *
+ * See also the tutorial rf709_BarlowBeeston.C
+ */
#include "Riostream.h"
#include "RooParamHistFunc.h"
@@ -37,6 +47,7 @@ RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataH
RooArgSet allVars ;
for (Int_t i=0 ; i<_dh.numEntries() ; i++) {
_dh.get(i) ;
+
const char* vname = Form("%s_gamma_bin_%i",GetName(),i) ;
RooRealVar* var = new RooRealVar(vname,vname,0,1000) ;
var->setVal(_relParam ? 1 : _dh.weight()) ;
@@ -237,22 +248,19 @@ Double_t RooParamHistFunc::analyticalIntegralWN(Int_t code, const RooArgSet* /*n
{
R__ASSERT(code==1) ;
- RooFIter iter = _p.fwdIterator() ;
- RooAbsReal* p ;
Double_t ret(0) ;
Int_t i(0) ;
- while((p=(RooAbsReal*)iter.next())) {
+ for (const auto param : _p) {
+ auto p = static_cast<const RooAbsReal*>(param);
Double_t bin = p->getVal() ;
if (_relParam) bin *= getNominal(i++) ;
ret += bin ;
}
// WVE fix this!!! Assume uniform binning for now!
- RooFIter xiter = _x.fwdIterator() ;
- RooAbsArg* obs ;
Double_t binV(1) ;
- while((obs=xiter.next())) {
- RooRealVar* xx = (RooRealVar*) obs ;
+ for (const auto obs : _x) {
+ auto xx = static_cast<const RooRealVar*>(obs);
binV *= (xx->getMax()-xx->getMin())/xx->numBins() ;
}
diff --git a/tutorials/roofit/rf706_histpdf.C b/tutorials/roofit/rf706_histpdf.C
index bf28d2725a3344c3e8fe024c241adf50780f7500..9c323ecfba3ddd33e230c2a6e1bcc2f0d98a5e3e 100644
--- a/tutorials/roofit/rf706_histpdf.C
+++ b/tutorials/roofit/rf706_histpdf.C
@@ -1,7 +1,7 @@
/// \file
/// \ingroup tutorial_roofit
/// \notebook -js
-/// Speecial p.d.f.'s: histogram based p.d.f.s and functions
+/// Special p.d.f.'s: histogram-based p.d.f.s and functions
///
/// \macro_image
/// \macro_output
diff --git a/tutorials/roofit/rf709_BarlowBeeston.C b/tutorials/roofit/rf709_BarlowBeeston.C
new file mode 100644
index 0000000000000000000000000000000000000000..61c1ce69f9accc8a9d992ba4b2de5c2b9438db7a
--- /dev/null
+++ b/tutorials/roofit/rf709_BarlowBeeston.C
@@ -0,0 +1,233 @@
+/// \file rf709_BarlowBeeston.C
+/// \ingroup tutorial_roofit
+/// \notebook -js
+/// Implementing the Barlow-Beeston method for taking into account the statistical uncertainty of a Monte-Carlo fit template.
+/// \macro_image
+/// \macro_output
+/// \macro_code
+/// \author 06/2019 - Stephan Hageboeck, CERN
+/// Based on a demo by Wouter Verkerke
+
+#include "RooRealVar.h"
+#include "RooGaussian.h"
+#include "RooUniform.h"
+#include "RooDataSet.h"
+#include "RooDataHist.h"
+#include "RooHistFunc.h"
+#include "RooRealSumPdf.h"
+#include "RooParamHistFunc.h"
+#include "RooHistConstraint.h"
+#include "RooProdPdf.h"
+#include "RooPlot.h"
+
+#include "TCanvas.h"
+#include "TPaveText.h"
+
+#include <iostream>
+#include <memory>
+
+using namespace RooFit;
+
+void rf709_BarlowBeeston()
+{
+ // First, construct a likelihood model with a Gaussian signal on top of a uniform background
+ RooRealVar x("x", "x", -20, 20);
+ x.setBins(25);
+
+ RooRealVar meanG("meanG", "meanG", 1, -10, 10);
+ RooRealVar sigG("sigG", "sigG", 1.5, -10, 10);
+ RooGaussian g("g", "Gauss", x, meanG, sigG);
+ RooUniform u("u", "Uniform", x);
+
+
+ // Generate the data to be fitted
+ std::unique_ptr<RooDataSet> sigData(g.generate(x, 50));
+ std::unique_ptr<RooDataSet> bkgData(u.generate(x, 1000));
+
+ RooDataSet sumData("sumData", "Gauss + Uniform", x);
+ sumData.append(*sigData);
+ sumData.append(*bkgData);
+
+
+ // Make histogram templates for signal and background.
+ // Let's take a signal distribution with low statistics and a more accurate
+ // background distribution.
+ // Normally, these come from Monte Carlo simulations, but we will just generate them.
+ std::unique_ptr<RooDataHist> dh_sig( g.generateBinned(x, 50) );
+ std::unique_ptr<RooDataHist> dh_bkg( u.generateBinned(x, 10000) );
+
+
+ // ***** Case 0 - 'Rigid templates' *****
+
+ // Construct histogram shapes for signal and background
+ RooHistFunc p_h_sig("p_h_sig","p_h_sig",x,*dh_sig);
+ RooHistFunc p_h_bkg("p_h_bkg","p_h_bkg",x,*dh_bkg);
+
+ // Construct scale factors for adding the two distributions
+ RooRealVar Asig0("Asig","Asig",1,0.01,5000);
+ RooRealVar Abkg0("Abkg","Abkg",1,0.01,5000);
+
+ // Construct the sum model
+ RooRealSumPdf model0("model0","model0",
+ RooArgList(p_h_sig,p_h_bkg),
+ RooArgList(Asig0,Abkg0),
+ true);
+
+
+
+ // ***** Case 1 - 'Barlow Beeston' *****
+
+ // Construct parameterized histogram shapes for signal and background
+ RooParamHistFunc p_ph_sig1("p_ph_sig","p_ph_sig",*dh_sig);
+ RooParamHistFunc p_ph_bkg1("p_ph_bkg","p_ph_bkg",*dh_bkg);
+
+ RooRealVar Asig1("Asig","Asig",1,0.01,5000);
+ RooRealVar Abkg1("Abkg","Abkg",1,0.01,5000);
+
+ // Construct the sum of these
+ RooRealSumPdf model_tmp("sp_ph", "sp_ph",
+ RooArgList(p_ph_sig1,p_ph_bkg1),
+ RooArgList(Asig1,Abkg1),
+ true);
+
+ // Construct the subsidiary poisson measurements constraining the histogram parameters
+ // These ensure that the bin contents of the histograms are only allowed to vary within
+ // the statistical uncertainty of the Monte Carlo.
+ RooHistConstraint hc_sig("hc_sig","hc_sig",p_ph_sig1);
+ RooHistConstraint hc_bkg("hc_bkg","hc_bkg",p_ph_bkg1);
+
+ // Construct the joint model with template PDFs and constraints
+ RooProdPdf model1("model1","model1",RooArgSet(hc_sig,hc_bkg),Conditional(model_tmp,x));
+
+
+
+ // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) *****
+
+ // Construct the histogram shapes, using the same parameters for signal and background
+ // This requires passing the first histogram to the second, so that their common parameters
+ // can be re-used.
+ // The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`.
+ // This allows bin 0 to fluctuate up and down.
+ // Then, the SAME parameters are connected to the background histogram, so the bins flucutate
+ // synchronously. This reduces the number of parameters.
+ RooParamHistFunc p_ph_sig2("p_ph_sig2", "p_ph_sig2", *dh_sig);
+ RooParamHistFunc p_ph_bkg2("p_ph_bkg2", "p_ph_bkg2", *dh_bkg, p_ph_sig2, true);
+
+ RooRealVar Asig2("Asig","Asig",1,0.01,5000);
+ RooRealVar Abkg2("Abkg","Abkg",1,0.01,5000);
+
+ // As before, construct the sum of signal2 and background2
+ RooRealSumPdf model2_tmp("sp_ph","sp_ph",
+ RooArgList(p_ph_sig2,p_ph_bkg2),
+ RooArgList(Asig2,Abkg2),
+ true);
+
+ // Construct the subsidiary poisson measurements constraining the statistical fluctuations
+ RooHistConstraint hc_sigbkg("hc_sigbkg","hc_sigbkg",RooArgSet(p_ph_sig2,p_ph_bkg2));
+
+ // Construct the joint model
+ RooProdPdf model2("model2","model2",hc_sigbkg, RooFit::Conditional(model2_tmp,x));
+
+
+
+ // ************ Fit all models to data and plot *********************
+
+ auto result0 = model0.fitTo(sumData, PrintLevel(0), Save());
+ auto result1 = model1.fitTo(sumData, PrintLevel(0), Save());
+ auto result2 = model2.fitTo(sumData, PrintLevel(0), Save());
+
+
+ TCanvas* can = new TCanvas("can", "", 1500, 600);
+ can->Divide(3,1);
+
+ TPaveText pt(-19.5, 1, -2, 25);
+ pt.SetFillStyle(0);
+ pt.SetBorderSize(0);
+
+
+ can->cd(1);
+ auto frame = x.frame(Title("No template uncertainties"));
+ // Plot data to enable automatic determination of model0 normalisation:
+ sumData.plotOn(frame);
+ model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0));
+ // Plot data again to show it on top of model0 error bands:
+ sumData.plotOn(frame);
+ // Plot model components
+ model0.plotOn(frame, LineColor(kBlue));
+ model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure));
+ model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed));
+ model0.paramOn(frame);
+
+ sigData->plotOn(frame, MarkerColor(kBlue));
+ frame->Draw();
+
+ for (auto text : {
+ "No template uncertainties",
+ "are taken into account.",
+ "This leads to low errors",
+ "for the parameters A, since",
+ "the only source of errors",
+ "are the data statistics."}) {
+ pt.AddText(text);
+ }
+ pt.DrawClone();
+
+
+ can->cd(2);
+ frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately"));
+ sumData.plotOn(frame);
+ model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1));
+ // Plot data again to show it on top of error bands:
+ sumData.plotOn(frame);
+ model1.plotOn(frame, LineColor(kBlue));
+ model1.plotOn(frame, Components(p_ph_sig1), LineColor(kAzure));
+ model1.plotOn(frame, Components(p_ph_bkg1), LineColor(kRed));
+ model1.paramOn(frame, Parameters(RooArgSet(Asig1, Abkg1)));
+
+ sigData->plotOn(frame, MarkerColor(kBlue));
+ frame->Draw();
+
+ pt.Clear();
+ for (auto text : {
+ "With gamma parameters, the",
+ "signal & background templates",
+ "can adapt to the data.",
+ "Note how the blue signal",
+ "template changes its shape.",
+ "This leads to higher errors",
+ "of the scale parameters A."}) {
+ pt.AddText(text);
+ }
+ pt.DrawClone();
+
+ can->cd(3);
+ frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)"));
+ sumData.plotOn(frame);
+ model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2));
+ // Plot data again to show it on top of model0 error bands:
+ sumData.plotOn(frame);
+ model2.plotOn(frame, LineColor(kBlue));
+ model2.plotOn(frame, Components(p_ph_sig2), LineColor(kAzure));
+ model2.plotOn(frame, Components(p_ph_bkg2), LineColor(kRed));
+ model2.paramOn(frame, Parameters(RooArgSet(Asig2, Abkg2)));
+
+ sigData->plotOn(frame, MarkerColor(kBlue));
+ frame->Draw();
+
+ pt.Clear();
+ for (auto text : {
+ "When signal and background",
+ "template share one gamma para-",
+ "meter per bin, they adapt less.",
+ "The errors of the A parameters",
+ "also shrink slightly."}) {
+ pt.AddText(text);
+ }
+ pt.DrawClone();
+
+
+ std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl;
+ std::cout << "Asig [BB ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl;
+ std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl;
+
+}