Logo ROOT  
Reference Guide
rf709_BarlowBeeston.py
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## \date June 2021
13## \author Harshal Shende, Stephan Hageboeck (C++ version)
14
15
16import ROOT
17
18# First, construct a likelihood model with a Gaussian signal on top of a uniform background
19x = ROOT.RooRealVar("x", "x", -20, 20)
20x.setBins(25)
21
22meanG = ROOT.RooRealVar("meanG", "meanG", 1, -10, 10)
23sigG = ROOT.RooRealVar("sigG", "sigG", 1.5, -10, 10)
24g = ROOT.RooGaussian("g", "Gauss", x, meanG, sigG)
25u = ROOT.RooUniform("u", "Uniform", x)
26
27
28# Generate the data to be fitted
29sigData = g.generate(x, 50)
30bkgData = u.generate(x, 1000)
31
32sumData = ROOT.RooDataSet("sumData", "Gauss + Uniform", x)
33sumData.append(sigData)
34sumData.append(bkgData)
35
36
37# Make histogram templates for signal and background.
38# Let's take a signal distribution with low statistics and a more accurate
39# background distribution.
40# Normally, these come from Monte Carlo simulations, but we will just generate them.
41dh_sig = g.generateBinned(x, 50)
42dh_bkg = u.generateBinned(x, 10000)
43
44
45# Case 0 - 'Rigid templates'
46
47# Construct histogram shapes for signal and background
48p_h_sig = ROOT.RooHistFunc("p_h_sig", "p_h_sig", x, dh_sig)
49p_h_bkg = ROOT.RooHistFunc("p_h_bkg", "p_h_bkg", x, dh_bkg)
50
51# Construct scale factors for adding the two distributions
52Asig0 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000)
53Abkg0 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000)
54
55# Construct the sum model
56model0 = ROOT.RooRealSumPdf("model0", "model0", [p_h_sig, p_h_bkg], [Asig0, Abkg0], True)
57
58
59# Case 1 - 'Barlow Beeston'
60
61# Construct parameterized histogram shapes for signal and background
62p_ph_sig1 = ROOT.RooParamHistFunc("p_ph_sig", "p_ph_sig", dh_sig)
63p_ph_bkg1 = ROOT.RooParamHistFunc("p_ph_bkg", "p_ph_bkg", dh_bkg)
64
65Asig1 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000)
66Abkg1 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000)
67
68# Construct the sum of these
69model_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_ph_sig1, p_ph_bkg1], [Asig1, Abkg1], True)
70
71# Construct the subsidiary poisson measurements constraining the histogram parameters
72# These ensure that the bin contents of the histograms are only allowed to vary within
73# the statistical uncertainty of the Monte Carlo.
74hc_sig = ROOT.RooHistConstraint("hc_sig", "hc_sig", p_ph_sig1)
75hc_bkg = ROOT.RooHistConstraint("hc_bkg", "hc_bkg", p_ph_bkg1)
76
77# Construct the joint model with template PDFs and constraints
78model1 = ROOT.RooProdPdf("model1", "model1", {hc_sig, hc_bkg}, Conditional=(model_tmp, x))
79
80
81# Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples)
82
83# Construct the histogram shapes, using the same parameters for signal and background
84# This requires passing the first histogram to the second, so that their common parameters
85# can be re-used.
86# The first ParamHistFunc will create one parameter per bin, such as `p_ph_sig2_gamma_bin_0`.
87# This allows bin 0 to fluctuate up and down.
88# Then, the SAME parameters are connected to the background histogram, so the bins flucutate
89# synchronously. This reduces the number of parameters.
90p_ph_sig2 = ROOT.RooParamHistFunc("p_ph_sig2", "p_ph_sig2", dh_sig)
91p_ph_bkg2 = ROOT.RooParamHistFunc("p_ph_bkg2", "p_ph_bkg2", dh_bkg, p_ph_sig2, True)
92
93Asig2 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000)
94Abkg2 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000)
95
96# As before, construct the sum of signal2 and background2
97model2_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_ph_sig2, p_ph_bkg2], [Asig2, Abkg2], True)
98
99# Construct the subsidiary poisson measurements constraining the statistical fluctuations
100hc_sigbkg = ROOT.RooHistConstraint("hc_sigbkg", "hc_sigbkg", {p_ph_sig2, p_ph_bkg2})
101
102# Construct the joint model
103model2 = ROOT.RooProdPdf("model2", "model2", hc_sigbkg, Conditional=(model2_tmp, x))
104
105
106# ************ Fit all models to data and plot *********************
107
108result0 = model0.fitTo(sumData, PrintLevel=0, Save=True)
109result1 = model1.fitTo(sumData, PrintLevel=0, Save=True)
110result2 = model2.fitTo(sumData, PrintLevel=0, Save=True)
111
112
113can = ROOT.TCanvas("can", "", 1500, 600)
114can.Divide(3, 1)
115
116pt = ROOT.TPaveText(-19.5, 1, -2, 25)
117pt.SetFillStyle(0)
118pt.SetBorderSize(0)
119
120
121can.cd(1)
122frame = x.frame(Title="No template uncertainties")
123# Plot data to enable automatic determination of model0 normalisation:
124sumData.plotOn(frame)
125model0.plotOn(frame, LineColor="b", VisualizeError=result0)
126# Plot data again to show it on top of model0 error bands:
127sumData.plotOn(frame)
128# Plot model components
129model0.plotOn(frame, LineColor="b")
130p_ph_sig_set = {p_h_sig}
131p_ph_bkg_set = {p_h_bkg}
132model0.plotOn(frame, Components=p_ph_sig_set, LineColor="kAzure")
133model0.plotOn(frame, Components=p_ph_bkg_set, LineColor="r")
134model0.paramOn(frame)
135
136sigData.plotOn(frame, MarkerColor="b")
137frame.Draw()
138
139pt_text1 = [
140 "No template uncertainties",
141 "are taken into account.",
142 "This leads to low errors",
143 "for the parameters A, since",
144 "the only source of errors",
145 "are the data statistics.",
146]
147for text in pt_text1:
148 pt.AddText(text)
149
150pt.DrawClone()
151
152
153can.cd(2)
154frame = x.frame(Title="Barlow Beeston for Sig & Bkg separately")
155sumData.plotOn(frame)
156model1.plotOn(frame, LineColor="b", VisualizeError=result1)
157# Plot data again to show it on top of error bands:
158sumData.plotOn(frame)
159model1.plotOn(frame, LineColor="b")
160p_ph_sig1_set = {p_ph_sig1}
161p_ph_bkg1_set = {p_ph_bkg1}
162model1.plotOn(frame, Components=p_ph_sig1_set, LineColor="kAzure")
163model1.plotOn(frame, Components=p_ph_bkg1_set, LineColor="r")
164model1.paramOn(frame, Parameters={Asig1, Abkg1})
165
166sigData.plotOn(frame, MarkerColor="b")
167frame.Draw()
168
169pt.Clear()
170pt_text2 = [
171 "With gamma parameters, the",
172 "signal & background templates",
173 "can adapt to the data.",
174 "Note how the blue signal",
175 "template changes its shape.",
176 "This leads to higher errors",
177 "of the scale parameters A.",
178]
179
180for text in pt_text2:
181 pt.AddText(text)
182
183pt.DrawClone()
184
185can.cd(3)
186frame = x.frame(Title="Barlow Beeston light for (Sig+Bkg)")
187sumData.plotOn(frame)
188model2.plotOn(frame, LineColor="b", VisualizeError=result2)
189# Plot data again to show it on top of model0 error bands:
190sumData.plotOn(frame)
191model2.plotOn(frame, LineColor="b")
192p_ph_sig2_set = {p_ph_sig2}
193p_ph_bkg2_set = {p_ph_bkg2}
194model2.plotOn(frame, Components=p_ph_sig2_set, LineColor="kAzure")
195model2.plotOn(frame, Components=p_ph_bkg2_set, LineColor="r")
196model2.paramOn(frame, Parameters={Asig2, Abkg2})
197
198sigData.plotOn(frame, MarkerColor="b")
199frame.Draw()
200
201pt.Clear()
202pt_text3 = [
203 "When signal and background",
204 "template share one gamma para-",
205 "meter per bin, they adapt less.",
206 "The errors of the A parameters",
207 "also shrink slightly.",
208]
209for text in pt_text3:
210 pt.AddText(text)
211pt.DrawClone()
212
213
214print("Asig [normal ] = {} +/- {}".format(Asig0.getVal(), Asig0.getError()))
215print("Asig [BB ] = {} +/- {}".format(Asig1.getVal(), Asig1.getError()))
216print("Asig [BBlight] = {} +/- {}".format(Asig2.getVal(), Asig2.getError()))
217
218can.SaveAs("rf709_BarlowBeeston.png")