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