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