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