-
Notifications
You must be signed in to change notification settings - Fork 18
/
rf608_fitresultaspdf.py
120 lines (95 loc) · 4.21 KB
/
rf608_fitresultaspdf.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
#####################################
#
# 'LIKELIHOOD AND MINIMIZATION' ROOT.RooFit tutorial macro #608
#
# Representing the parabolic approximation of the fit as
# a multi-variate Gaussian on the parameters of the fitted p.d.f.
#
#
# 07/2008 - Wouter Verkerke
#
# /
import ROOT
def rf608_fitresultaspdf():
# C r e a t e m o d e l a n d d a t a s e t
# -----------------------------------------------
# Observable
x = ROOT.RooRealVar("x", "x", -20, 20)
# Model (intentional strong correlations)
mean = ROOT.RooRealVar("mean", "mean of g1 and g2", 0, -1, 1)
sigma_g1 = ROOT.RooRealVar("sigma_g1", "width of g1", 2)
g1 = ROOT.RooGaussian("g1", "g1", x, mean, sigma_g1)
sigma_g2 = ROOT.RooRealVar("sigma_g2", "width of g2", 4, 3.0, 5.0)
g2 = ROOT.RooGaussian("g2", "g2", x, mean, sigma_g2)
frac = ROOT.RooRealVar("frac", "frac", 0.5, 0.0, 1.0)
model = ROOT.RooAddPdf("model", "model", ROOT.RooArgList(g1, g2), ROOT.RooArgList(frac))
# Generate 1000 events
data = model.generate(ROOT.RooArgSet(x), 1000)
# F i t m o d e l t o d a t a
# ----------------------------------
r = model.fitTo(data, ROOT.RooFit.Save())
# C r e a t e M V G a u s s i a n p d f o f f i t t e d p a r a m e t e r s
# ------------------------------------------------------------------------------------
parabPdf = r.createHessePdf(ROOT.RooArgSet(frac, mean, sigma_g2))
# S o m e e x e c e r c i s e s w i t h t h e p a r a m e t e r p d f
# -----------------------------------------------------------------------------
# Generate 100K points in the parameter space, from the MVGaussian p.d.f.
d = parabPdf.generate(ROOT.RooArgSet(mean, sigma_g2, frac), 100000)
# Sample a 3-D histogram of the p.d.f. to be visualized as an error
# ellipsoid using the GLISO draw option
hh_3d = parabPdf.createHistogram("mean,sigma_g2,frac", 25, 25, 25)
hh_3d.SetFillColor(ROOT.kBlue)
# Project 3D parameter p.d.f. down to 3 permutations of two-dimensional p.d.f.s
# ROOT.The integrations corresponding to these projections are performed analytically
# by the MV Gaussian p.d.f.
pdf_sigmag2_frac = parabPdf.createProjection(ROOT.RooArgSet(mean))
pdf_mean_frac = parabPdf.createProjection(ROOT.RooArgSet(sigma_g2))
pdf_mean_sigmag2 = parabPdf.createProjection(ROOT.RooArgSet(frac))
# Make 2D plots of the 3 two-dimensional p.d.f. projections
hh_sigmag2_frac = pdf_sigmag2_frac.createHistogram("sigma_g2,frac", 50, 50)
hh_mean_frac = pdf_mean_frac.createHistogram("mean,frac", 50, 50)
hh_mean_sigmag2 = pdf_mean_sigmag2.createHistogram("mean,sigma_g2", 50, 50)
hh_mean_frac.SetLineColor(ROOT.kBlue)
hh_sigmag2_frac.SetLineColor(ROOT.kBlue)
hh_mean_sigmag2.SetLineColor(ROOT.kBlue)
# Draw the 'sigar'
ROOT.gStyle.SetCanvasPreferGL(True)
ROOT.gStyle.SetPalette(1)
c1 = ROOT.TCanvas("rf608_fitresultaspdf_1", "rf608_fitresultaspdf_1", 600, 600)
hh_3d.Draw("gliso")
c1.SaveAs("rf608_fitresultaspdf_1.png")
# Draw the 2D projections of the 3D p.d.f.
c2 = ROOT.TCanvas("rf608_fitresultaspdf_2",
"rf608_fitresultaspdf_2", 900, 600)
c2.Divide(3, 2)
c2.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
hh_mean_sigmag2.GetZaxis().SetTitleOffset(1.4)
hh_mean_sigmag2.Draw("surf3")
c2.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
hh_sigmag2_frac.GetZaxis().SetTitleOffset(1.4)
hh_sigmag2_frac.Draw("surf3")
c2.cd(3)
ROOT.gPad.SetLeftMargin(0.15)
hh_mean_frac.GetZaxis().SetTitleOffset(1.4)
hh_mean_frac.Draw("surf3")
# Draw the distributions of parameter points sampled from the p.d.f.
tmp1 = d.createHistogram(mean, sigma_g2, 50, 50)
tmp2 = d.createHistogram(sigma_g2, frac, 50, 50)
tmp3 = d.createHistogram(mean, frac, 50, 50)
c2.cd(4)
ROOT.gPad.SetLeftMargin(0.15)
tmp1.GetZaxis().SetTitleOffset(1.4)
tmp1.Draw("lego3")
c2.cd(5)
ROOT.gPad.SetLeftMargin(0.15)
tmp2.GetZaxis().SetTitleOffset(1.4)
tmp2.Draw("lego3")
c2.cd(6)
ROOT.gPad.SetLeftMargin(0.15)
tmp3.GetZaxis().SetTitleOffset(1.4)
tmp3.Draw("lego3")
c2.SaveAs("rf608_fitresultaspdf_2.png")
if __name__ == "__main__":
rf608_fitresultaspdf()