-
Notifications
You must be signed in to change notification settings - Fork 18
/
rf610_visualerror.py
182 lines (152 loc) · 7.15 KB
/
rf610_visualerror.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
# /
#
# 'LIKELIHOOD AND MINIMIZATION' ROOT.RooFit tutorial macro #610
#
# Visualization of errors from a covariance matrix
#
#
#
# 04/2009 - Wouter Verkerke
#
# /
import ROOT
def rf610_visualerror():
# S e t u p e x a m p l e f i t
# ---------------------------------------
# Create sum of two Gaussians p.d.f. with factory
x = ROOT.RooRealVar("x", "x", -10, 10)
m = ROOT.RooRealVar("m", "m", 0, -10, 10)
s = ROOT.RooRealVar("s", "s", 2, 1, 50)
sig = ROOT.RooGaussian("sig", "sig", x, m, s)
m2 = ROOT.RooRealVar("m2", "m2", -1, -10, 10)
s2 = ROOT.RooRealVar("s2", "s2", 6, 1, 50)
bkg = ROOT.RooGaussian("bkg", "bkg", x, m2, s2)
fsig = ROOT.RooRealVar("fsig", "fsig", 0.33, 0, 1)
model = ROOT.RooAddPdf("model", "model", ROOT.RooArgList(
sig, bkg), ROOT.RooArgList(fsig))
# Create binned dataset
x.setBins(25)
d = model.generateBinned(ROOT.RooArgSet(x), 1000)
# Perform fit and save fit result
r = model.fitTo(d, ROOT.RooFit.Save())
# V i s u a l i z e f i t e r r o r
# -------------------------------------
# Make plot frame
frame = x.frame(ROOT.RooFit.Bins(40), ROOT.RooFit.Title(
"P.d.f with visualized 1-sigma error band"))
d.plotOn(frame)
# Visualize 1-sigma error encoded in fit result 'r' as orange band using linear error propagation
# ROOT.This results in an error band that is by construction symmetric
#
# ROOT.The linear error is calculated as
# error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)
#
# where F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2,
#
# with f(x) = the plotted curve
# 'da' = error taken from the fit result
# Corr(a,a') = the correlation matrix from the fit result
# Z = requested significance 'Z sigma band'
#
# ROOT.The linear method is fast (required 2*N evaluations of the curve, N is the number of parameters),
# but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and
# Gaussian approximations made
#
model.plotOn(frame, ROOT.RooFit.VisualizeError(
r, 1), ROOT.RooFit.FillColor(ROOT.kOrange))
# Calculate error using sampling method and visualize as dashed red line.
#
# In self method a number of curves is calculated with variations of the parameter values, sampled
# from a multi-variate Gaussian p.d.f. that is constructed from the fit results covariance matrix.
# ROOT.The error(x) is determined by calculating a central interval that capture N% of the variations
# for each valye of x, N% is controlled by Z (i.e. Z=1 gives N=68%). ROOT.The number of sampling curves
# is chosen to be such that at least 100 curves are expected to be outside the N% interval, is minimally
# 100 (e.g. Z=1.Ncurve=356, Z=2.Ncurve=2156)) Intervals from the sampling method can be asymmetric,
# and may perform better in the presence of strong correlations, may take
# (much) longer to calculate
model.plotOn(frame, ROOT.RooFit.VisualizeError(r, 1, ROOT.kFALSE), ROOT.RooFit.DrawOption(
"L"), ROOT.RooFit.LineWidth(2), ROOT.RooFit.LineColor(ROOT.kRed))
# Perform the same type of error visualization on the background component only.
# ROOT.The VisualizeError() option can generally applied to _any_ kind of
# plot (components, asymmetries, etc..)
model.plotOn(frame, ROOT.RooFit.VisualizeError(r, 1),
ROOT.RooFit.FillColor(ROOT.kOrange), ROOT.RooFit.Components("bkg"))
model.plotOn(frame, ROOT.RooFit.VisualizeError(r, 1, ROOT.kFALSE), ROOT.RooFit.DrawOption("L"), ROOT.RooFit.LineWidth(
2), ROOT.RooFit.LineColor(ROOT.kRed), ROOT.RooFit.Components("bkg"), ROOT.RooFit.LineStyle(ROOT.kDashed))
# Overlay central value
model.plotOn(frame)
model.plotOn(frame, ROOT.RooFit.Components("bkg"),
ROOT.RooFit.LineStyle(ROOT.kDashed))
d.plotOn(frame)
frame.SetMinimum(0)
# V i s u a l i z e p a r t i a l f i t e r r o r
# ------------------------------------------------------
# Make plot frame
frame2 = x.frame(ROOT.RooFit.Bins(40), ROOT.RooFit.Title(
"Visualization of 2-sigma partial error from (m,m2)"))
# Visualize partial error. For partial error visualization the covariance matrix is first reduced as follows
# ___ -1
# Vred = V22 = V11 - V12 * V22 * V21
#
# Where V11,V12,V21, represent a block decomposition of the covariance matrix into observables that
# are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), V22bar
# is the Shur complement of V22, as shown above
#
# (Note that Vred is _not_ a simple sub-matrix of V)
# Propagate partial error due to shape parameters (m,m2) using linear and
# sampling method
model.plotOn(frame2, ROOT.RooFit.VisualizeError(
r, ROOT.RooArgSet(m, m2), 2), ROOT.RooFit.FillColor(ROOT.kCyan))
model.plotOn(frame2, ROOT.RooFit.Components("bkg"), ROOT.RooFit.VisualizeError(
r, ROOT.RooArgSet(m, m2), 2), ROOT.RooFit.FillColor(ROOT.kCyan))
model.plotOn(frame2)
model.plotOn(frame2, ROOT.RooFit.Components("bkg"),
ROOT.RooFit.LineStyle(ROOT.kDashed))
frame2.SetMinimum(0)
# Make plot frame
frame3 = x.frame(ROOT.RooFit.Bins(40), ROOT.RooFit.Title(
"Visualization of 2-sigma partial error from (s,s2)"))
# Propagate partial error due to yield parameter using linear and sampling
# method
model.plotOn(frame3, ROOT.RooFit.VisualizeError(
r, ROOT.RooArgSet(s, s2), 2), ROOT.RooFit.FillColor(ROOT.kGreen))
model.plotOn(frame3, ROOT.RooFit.Components("bkg"), ROOT.RooFit.VisualizeError(
r, ROOT.RooArgSet(s, s2), 2), ROOT.RooFit.FillColor(ROOT.kGreen))
model.plotOn(frame3)
model.plotOn(frame3, ROOT.RooFit.Components("bkg"),
ROOT.RooFit.LineStyle(ROOT.kDashed))
frame3.SetMinimum(0)
# Make plot frame
frame4 = x.frame(ROOT.RooFit.Bins(40), ROOT.RooFit.Title(
"Visualization of 2-sigma partial error from fsig"))
# Propagate partial error due to yield parameter using linear and sampling
# method
model.plotOn(frame4, ROOT.RooFit.VisualizeError(
r, ROOT.RooArgSet(fsig), 2), ROOT.RooFit.FillColor(ROOT.kMagenta))
model.plotOn(frame4, ROOT.RooFit.Components("bkg"), ROOT.RooFit.VisualizeError(
r, ROOT.RooArgSet(fsig), 2), ROOT.RooFit.FillColor(ROOT.kMagenta))
model.plotOn(frame4)
model.plotOn(frame4, ROOT.RooFit.Components("bkg"),
ROOT.RooFit.LineStyle(ROOT.kDashed))
frame4.SetMinimum(0)
c = ROOT.TCanvas("rf610_visualerror", "rf610_visualerror", 800, 800)
c.Divide(2, 2)
c.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
frame.GetYaxis().SetTitleOffset(1.4)
frame.Draw()
c.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
frame2.GetYaxis().SetTitleOffset(1.6)
frame2.Draw()
c.cd(3)
ROOT.gPad.SetLeftMargin(0.15)
frame3.GetYaxis().SetTitleOffset(1.6)
frame3.Draw()
c.cd(4)
ROOT.gPad.SetLeftMargin(0.15)
frame4.GetYaxis().SetTitleOffset(1.6)
frame4.Draw()
c.SaveAs("rf610_visualerror.png")
if __name__ == "__main__":
rf610_visualerror()