Logo ROOT  
Reference Guide
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 
11 import ROOT
12 
13 
14 # Setup example fit
15 # ---------------------------------------
16 
17 # Create sum of two Gaussians pdf with factory
18 x = ROOT.RooRealVar("x", "x", -10, 10)
19 
20 m = ROOT.RooRealVar("m", "m", 0, -10, 10)
21 s = ROOT.RooRealVar("s", "s", 2, 1, 50)
22 sig = ROOT.RooGaussian("sig", "sig", x, m, s)
23 
24 m2 = ROOT.RooRealVar("m2", "m2", -1, -10, 10)
25 s2 = ROOT.RooRealVar("s2", "s2", 6, 1, 50)
26 bkg = ROOT.RooGaussian("bkg", "bkg", x, m2, s2)
27 
28 fsig = ROOT.RooRealVar("fsig", "fsig", 0.33, 0, 1)
29 model = ROOT.RooAddPdf("model", "model", ROOT.RooArgList(
30  sig, bkg), ROOT.RooArgList(fsig))
31 
32 # Create binned dataset
33 x.setBins(25)
34 d = model.generateBinned(ROOT.RooArgSet(x), 1000)
35 
36 # Perform fit and save fit result
37 r = model.fitTo(d, ROOT.RooFit.Save())
38 
39 # Visualize fit error
40 # -------------------------------------
41 
42 # Make plot frame
43 frame = x.frame(ROOT.RooFit.Bins(40), ROOT.RooFit.Title(
44  "P.d.f with visualized 1-sigma error band"))
45 d.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 #
64 model.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
77 model.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..)
91 model.plotOn(
92  frame, ROOT.RooFit.VisualizeError(
93  r, 1), ROOT.RooFit.FillColor(
94  ROOT.kOrange), ROOT.RooFit.Components("bkg"))
95 model.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
110 model.plotOn(frame)
111 model.plotOn(frame, ROOT.RooFit.Components("bkg"),
112  ROOT.RooFit.LineStyle(ROOT.kDashed))
113 d.plotOn(frame)
114 frame.SetMinimum(0)
115 
116 # Visualize partial fit error
117 # ------------------------------------------------------
118 
119 # Make plot frame
120 frame2 = 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
135 model.plotOn(frame2, ROOT.RooFit.VisualizeError(
136  r, ROOT.RooArgSet(m, m2), 2), ROOT.RooFit.FillColor(ROOT.kCyan))
137 model.plotOn(frame2, ROOT.RooFit.Components("bkg"), ROOT.RooFit.VisualizeError(
138  r, ROOT.RooArgSet(m, m2), 2), ROOT.RooFit.FillColor(ROOT.kCyan))
139 
140 model.plotOn(frame2)
141 model.plotOn(frame2, ROOT.RooFit.Components("bkg"),
142  ROOT.RooFit.LineStyle(ROOT.kDashed))
143 frame2.SetMinimum(0)
144 
145 # Make plot frame
146 frame3 = 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
151 model.plotOn(frame3, ROOT.RooFit.VisualizeError(
152  r, ROOT.RooArgSet(s, s2), 2), ROOT.RooFit.FillColor(ROOT.kGreen))
153 model.plotOn(frame3, ROOT.RooFit.Components("bkg"), ROOT.RooFit.VisualizeError(
154  r, ROOT.RooArgSet(s, s2), 2), ROOT.RooFit.FillColor(ROOT.kGreen))
155 
156 model.plotOn(frame3)
157 model.plotOn(frame3, ROOT.RooFit.Components("bkg"),
158  ROOT.RooFit.LineStyle(ROOT.kDashed))
159 frame3.SetMinimum(0)
160 
161 # Make plot frame
162 frame4 = 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
167 model.plotOn(frame4, ROOT.RooFit.VisualizeError(
168  r, ROOT.RooArgSet(fsig), 2), ROOT.RooFit.FillColor(ROOT.kMagenta))
169 model.plotOn(frame4, ROOT.RooFit.Components("bkg"), ROOT.RooFit.VisualizeError(
170  r, ROOT.RooArgSet(fsig), 2), ROOT.RooFit.FillColor(ROOT.kMagenta))
171 
172 model.plotOn(frame4)
173 model.plotOn(frame4, ROOT.RooFit.Components("bkg"),
174  ROOT.RooFit.LineStyle(ROOT.kDashed))
175 frame4.SetMinimum(0)
176 
177 c = ROOT.TCanvas("rf610_visualerror", "rf610_visualerror", 800, 800)
178 c.Divide(2, 2)
179 c.cd(1)
180 ROOT.gPad.SetLeftMargin(0.15)
181 frame.GetYaxis().SetTitleOffset(1.4)
182 frame.Draw()
183 c.cd(2)
184 ROOT.gPad.SetLeftMargin(0.15)
185 frame2.GetYaxis().SetTitleOffset(1.6)
186 frame2.Draw()
187 c.cd(3)
188 ROOT.gPad.SetLeftMargin(0.15)
189 frame3.GetYaxis().SetTitleOffset(1.6)
190 frame3.Draw()
191 c.cd(4)
192 ROOT.gPad.SetLeftMargin(0.15)
193 frame4.GetYaxis().SetTitleOffset(1.6)
194 frame4.Draw()
195 
196 c.SaveAs("rf610_visualerror.png")