Logo ROOT  
Reference Guide
rf709_BarlowBeeston.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roofit
3/// \notebook -js
4///
5/// Implementing the Barlow-Beeston method for taking into account the statistical
6/// uncertainty of a Monte-Carlo fit template.
7///
8/// \macro_image
9/// \macro_output
10/// \macro_code
11///
12/// Based on a demo by Wouter Verkerke
13//
14/// \date 06/2019
15/// \author Stephan Hageboeck, CERN
16
17#include "RooRealVar.h"
18#include "RooGaussian.h"
19#include "RooUniform.h"
20#include "RooDataSet.h"
21#include "RooDataHist.h"
22#include "RooHistFunc.h"
23#include "RooRealSumPdf.h"
24#include "RooParamHistFunc.h"
25#include "RooHistConstraint.h"
26#include "RooProdPdf.h"
27#include "RooPlot.h"
28
29#include "TCanvas.h"
30#include "TPaveText.h"
31
32#include <iostream>
33#include <memory>
34
35using namespace RooFit;
36
37void rf709_BarlowBeeston()
38{
39 // First, construct a likelihood model with a Gaussian signal on top of a uniform background
40 RooRealVar x("x", "x", -20, 20);
41 x.setBins(25);
42
43 RooRealVar meanG("meanG", "meanG", 1, -10, 10);
44 RooRealVar sigG("sigG", "sigG", 1.5, -10, 10);
45 RooGaussian g("g", "Gauss", x, meanG, sigG);
46 RooUniform u("u", "Uniform", x);
47
48
49 // Generate the data to be fitted
50 std::unique_ptr<RooDataSet> sigData(g.generate(x, 50));
51 std::unique_ptr<RooDataSet> bkgData(u.generate(x, 1000));
52
53 RooDataSet sumData("sumData", "Gauss + Uniform", x);
54 sumData.append(*sigData);
55 sumData.append(*bkgData);
56
57
58 // Make histogram templates for signal and background.
59 // Let's take a signal distribution with low statistics and a more accurate
60 // background distribution.
61 // Normally, these come from Monte Carlo simulations, but we will just generate them.
62 std::unique_ptr<RooDataHist> dh_sig( g.generateBinned(x, 50) );
63 std::unique_ptr<RooDataHist> dh_bkg( u.generateBinned(x, 10000) );
64
65
66 // ***** Case 0 - 'Rigid templates' *****
67
68 // Construct histogram shapes for signal and background
69 RooHistFunc p_h_sig("p_h_sig","p_h_sig",x,*dh_sig);
70 RooHistFunc p_h_bkg("p_h_bkg","p_h_bkg",x,*dh_bkg);
71
72 // Construct scale factors for adding the two distributions
73 RooRealVar Asig0("Asig","Asig",1,0.01,5000);
74 RooRealVar Abkg0("Abkg","Abkg",1,0.01,5000);
75
76 // Construct the sum model
77 RooRealSumPdf model0("model0","model0",
78 RooArgList(p_h_sig,p_h_bkg),
79 RooArgList(Asig0,Abkg0),
80 true);
81
82
83
84 // ***** Case 1 - 'Barlow Beeston' *****
85
86 // Construct parameterized histogram shapes for signal and background
87 RooParamHistFunc p_ph_sig1("p_ph_sig","p_ph_sig",*dh_sig);
88 RooParamHistFunc p_ph_bkg1("p_ph_bkg","p_ph_bkg",*dh_bkg);
89
90 RooRealVar Asig1("Asig","Asig",1,0.01,5000);
91 RooRealVar Abkg1("Abkg","Abkg",1,0.01,5000);
92
93 // Construct the sum of these
94 RooRealSumPdf model_tmp("sp_ph", "sp_ph",
95 RooArgList(p_ph_sig1,p_ph_bkg1),
96 RooArgList(Asig1,Abkg1),
97 true);
98
99 // Construct the subsidiary poisson measurements constraining the histogram parameters
100 // These ensure that the bin contents of the histograms are only allowed to vary within
101 // the statistical uncertainty of the Monte Carlo.
102 RooHistConstraint hc_sig("hc_sig","hc_sig",p_ph_sig1);
103 RooHistConstraint hc_bkg("hc_bkg","hc_bkg",p_ph_bkg1);
104
105 // Construct the joint model with template PDFs and constraints
106 RooProdPdf model1("model1","model1",RooArgSet(hc_sig,hc_bkg),Conditional(model_tmp,x));
107
108
109
110 // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) *****
111
112 // Construct the histogram shapes, using the same parameters for signal and background
113 // This requires passing the first histogram to the second, so that their common parameters
114 // can be re-used.
115 // The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`.
116 // This allows bin 0 to fluctuate up and down.
117 // Then, the SAME parameters are connected to the background histogram, so the bins flucutate
118 // synchronously. This reduces the number of parameters.
119 RooParamHistFunc p_ph_sig2("p_ph_sig2", "p_ph_sig2", *dh_sig);
120 RooParamHistFunc p_ph_bkg2("p_ph_bkg2", "p_ph_bkg2", *dh_bkg, p_ph_sig2, true);
121
122 RooRealVar Asig2("Asig","Asig",1,0.01,5000);
123 RooRealVar Abkg2("Abkg","Abkg",1,0.01,5000);
124
125 // As before, construct the sum of signal2 and background2
126 RooRealSumPdf model2_tmp("sp_ph","sp_ph",
127 RooArgList(p_ph_sig2,p_ph_bkg2),
128 RooArgList(Asig2,Abkg2),
129 true);
130
131 // Construct the subsidiary poisson measurements constraining the statistical fluctuations
132 RooHistConstraint hc_sigbkg("hc_sigbkg","hc_sigbkg",RooArgSet(p_ph_sig2,p_ph_bkg2));
133
134 // Construct the joint model
135 RooProdPdf model2("model2","model2",hc_sigbkg, RooFit::Conditional(model2_tmp,x));
136
137
138
139 // ************ Fit all models to data and plot *********************
140
141 auto result0 = model0.fitTo(sumData, PrintLevel(0), Save());
142 auto result1 = model1.fitTo(sumData, PrintLevel(0), Save());
143 auto result2 = model2.fitTo(sumData, PrintLevel(0), Save());
144
145
146 TCanvas* can = new TCanvas("can", "", 1500, 600);
147 can->Divide(3,1);
148
149 TPaveText pt(-19.5, 1, -2, 25);
150 pt.SetFillStyle(0);
151 pt.SetBorderSize(0);
152
153
154 can->cd(1);
155 auto frame = x.frame(Title("No template uncertainties"));
156 // Plot data to enable automatic determination of model0 normalisation:
157 sumData.plotOn(frame);
158 model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0));
159 // Plot data again to show it on top of model0 error bands:
160 sumData.plotOn(frame);
161 // Plot model components
162 model0.plotOn(frame, LineColor(kBlue));
163 model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure));
164 model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed));
165 model0.paramOn(frame);
166
167 sigData->plotOn(frame, MarkerColor(kBlue));
168 frame->Draw();
169
170 for (auto text : {
171 "No template uncertainties",
172 "are taken into account.",
173 "This leads to low errors",
174 "for the parameters A, since",
175 "the only source of errors",
176 "are the data statistics."}) {
177 pt.AddText(text);
178 }
179 pt.DrawClone();
180
181
182 can->cd(2);
183 frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately"));
184 sumData.plotOn(frame);
185 model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1));
186 // Plot data again to show it on top of error bands:
187 sumData.plotOn(frame);
188 model1.plotOn(frame, LineColor(kBlue));
189 model1.plotOn(frame, Components(p_ph_sig1), LineColor(kAzure));
190 model1.plotOn(frame, Components(p_ph_bkg1), LineColor(kRed));
191 model1.paramOn(frame, Parameters(RooArgSet(Asig1, Abkg1)));
192
193 sigData->plotOn(frame, MarkerColor(kBlue));
194 frame->Draw();
195
196 pt.Clear();
197 for (auto text : {
198 "With gamma parameters, the",
199 "signal & background templates",
200 "can adapt to the data.",
201 "Note how the blue signal",
202 "template changes its shape.",
203 "This leads to higher errors",
204 "of the scale parameters A."}) {
205 pt.AddText(text);
206 }
207 pt.DrawClone();
208
209 can->cd(3);
210 frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)"));
211 sumData.plotOn(frame);
212 model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2));
213 // Plot data again to show it on top of model0 error bands:
214 sumData.plotOn(frame);
215 model2.plotOn(frame, LineColor(kBlue));
216 model2.plotOn(frame, Components(p_ph_sig2), LineColor(kAzure));
217 model2.plotOn(frame, Components(p_ph_bkg2), LineColor(kRed));
218 model2.paramOn(frame, Parameters(RooArgSet(Asig2, Abkg2)));
219
220 sigData->plotOn(frame, MarkerColor(kBlue));
221 frame->Draw();
222
223 pt.Clear();
224 for (auto text : {
225 "When signal and background",
226 "template share one gamma para-",
227 "meter per bin, they adapt less.",
228 "The errors of the A parameters",
229 "also shrink slightly."}) {
230 pt.AddText(text);
231 }
232 pt.DrawClone();
233
234
235 std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl;
236 std::cout << "Asig [BB ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl;
237 std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl;
238
239}
#define g(i)
Definition: RSha256.hxx:105
@ kRed
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
@ kAzure
Definition: Rtypes.h:65
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
Plain Gaussian p.d.f.
Definition: RooGaussian.h:25
The RooHistConstraint implements constraint terms for a binned PDF with statistical uncertainties.
RooHistFunc implements a real-valued function sampled from a multidimensional histogram.
Definition: RooHistFunc.h:29
A histogram function that assigns scale parameters to every bin.
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
The class RooRealSumPdf implements a PDF constructed from a sum of functions:
Definition: RooRealSumPdf.h:24
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
Flat p.d.f.
Definition: RooUniform.h:24
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
The Canvas class.
Definition: TCanvas.h:27
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:701
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad for instance with: gROOT->SetSelectedPad(gPad...
Definition: TObject.cxx:219
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1165
A Pave (see TPave) with text, lines or/and boxes inside.
Definition: TPaveText.h:21
virtual TText * AddText(Double_t x1, Double_t y1, const char *label)
Add a new Text line to this pavetext at given coordinates.
Definition: TPaveText.cxx:182
virtual void Clear(Option_t *option="")
Clear all lines in this pavetext.
Definition: TPaveText.cxx:208
virtual void SetBorderSize(Int_t bordersize=4)
Definition: TPave.h:73
TPaveText * pt
TText * text
RooCmdArg Parameters(const RooArgSet &params)
RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, Bool_t depsAreCond=kFALSE)
RooCmdArg Save(Bool_t flag=kTRUE)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg VisualizeError(const RooDataSet &paramData, Double_t Z=1)
RooCmdArg MarkerColor(Color_t color)
RooCmdArg Components(const RooArgSet &compSet)
RooCmdArg LineColor(Color_t color)
Double_t x[n]
Definition: legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
const char * Title
Definition: TXMLSetup.cxx:67