-
Notifications
You must be signed in to change notification settings - Fork 18
/
rf607_fitresult.py
134 lines (102 loc) · 4.21 KB
/
rf607_fitresult.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
#####################################
#
# 'LIKELIHOOD AND MINIMIZATION' ROOT.RooFit tutorial macro #607
#
# Demonstration of options of the ROOT.RooFitResult class
#
#
#
# 07/2008 - Wouter Verkerke
#
# /
import ROOT
def rf607_fitresult():
# C r e a t e p d f , a t a
# --------------------------------
# Declare observable x
x = ROOT.RooRealVar("x", "x", 0, 10)
# Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and
# their parameters
mean = ROOT.RooRealVar("mean", "mean of gaussians", 5, -10, 10)
sigma1 = ROOT.RooRealVar("sigma1", "width of gaussians", 0.5, 0.1, 10)
sigma2 = ROOT.RooRealVar("sigma2", "width of gaussians", 1, 0.1, 10)
sig1 = ROOT.RooGaussian("sig1", "Signal component 1", x, mean, sigma1)
sig2 = ROOT.RooGaussian("sig2", "Signal component 2", x, mean, sigma2)
# Build Chebychev polynomial p.d.f.
a0 = ROOT.RooRealVar("a0", "a0", 0.5, 0., 1.)
a1 = ROOT.RooRealVar("a1", "a1", -0.2)
bkg = ROOT.RooChebychev("bkg", "Background", x, ROOT.RooArgList(a0, a1))
# Sum the signal components into a composite signal p.d.f.
sig1frac = ROOT.RooRealVar(
"sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.)
sig = ROOT.RooAddPdf(
"sig", "Signal", ROOT.RooArgList(sig1, sig2), ROOT.RooArgList(sig1frac))
# Sum the composite signal and background
bkgfrac = ROOT.RooRealVar("bkgfrac", "fraction of background", 0.5, 0., 1.)
model = ROOT.RooAddPdf(
"model", "g1+g2+a", ROOT.RooArgList(bkg, sig), ROOT.RooArgList(bkgfrac))
# Generate 1000 events
data = model.generate(ROOT.RooArgSet(x), 1000)
# F i t p d f t o d a t a , a v e f i t r e s u l t
# -------------------------------------------------------------
# Perform fit and save result
r = model.fitTo(data, ROOT.RooFit.Save())
# P r i n t f i t r e s u l t s
# ---------------------------------
# Summary printing: Basic info plus final values of floating fit parameters
r.Print()
# Verbose printing: Basic info, of constant parameters, and
# final values of floating parameters, correlations
r.Print("v")
# V i s u a l i z e c o r r e l a t i o n m a t r i x
# -------------------------------------------------------
# Construct 2D color plot of correlation matrix
ROOT.gStyle.SetOptStat(0)
ROOT.gStyle.SetPalette(1)
hcorr = r.correlationHist()
# Visualize ellipse corresponding to single correlation matrix element
frame = ROOT.RooPlot(sigma1, sig1frac, 0.45, 0.60, 0.65, 0.90)
frame.SetTitle("Covariance between sigma1 and sig1frac")
r.plotOn(frame, sigma1, sig1frac, "ME12ABHV")
# A c c e s s f i t r e s u l t i n f o r m a t i o n
# ---------------------------------------------------------
# Access basic information
print "EDM = ", r.edm()
print "-log(L) minimum = ", r.minNll()
# Access list of final fit parameter values
print "final value of floating parameters"
r.floatParsFinal().Print("s")
# Access correlation matrix elements
print "correlation between sig1frac and a0 is ", r.correlation(
sig1frac, a0)
print "correlation between bkgfrac and mean is ", r.correlation(
"bkgfrac", "mean")
# Extract covariance and correlation matrix as ROOT.TMatrixDSym
cor = r.correlationMatrix()
cov = r.covarianceMatrix()
# Print correlation, matrix
print "correlation matrix"
cor.Print()
print "covariance matrix"
cov.Print()
# P e r s i s t f i t r e s u l t i n r o o t f i l e
# -------------------------------------------------------------
# Open ROOT file save save result
f = ROOT.TFile("rf607_fitresult.root", "RECREATE")
r.Write("rf607")
f.Close()
# In a clean ROOT session retrieve the persisted fit result as follows:
# r = gDirectory.Get("rf607")
c = ROOT.TCanvas("rf607_fitresult", "rf607_fitresult", 800, 400)
c.Divide(2)
c.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
hcorr.GetYaxis().SetTitleOffset(1.4)
hcorr.Draw("colz")
c.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
frame.GetYaxis().SetTitleOffset(1.6)
frame.Draw()
c.SaveAs("rf607_fitresult.png")
if __name__ == "__main__":
rf607_fitresult()