-
Notifications
You must be signed in to change notification settings - Fork 18
/
rf510_wsnamedsets.py
141 lines (107 loc) · 5.19 KB
/
rf510_wsnamedsets.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
135
136
137
138
139
140
141
# /
#
# 'ORGANIZATION AND SIMULTANEOUS FITS' ROOT.RooFit tutorial macro #510
#
# Working with named parameter sets and parameter snapshots in
# workspaces
#
# 04/2009 - Wouter Verkerke
#
# /
import ROOT
def rf510_wsnamedsets():
# C r e a t e m o d e l a n d d a t a s e t
# -----------------------------------------------
w = ROOT.RooWorkspace("w")
fillWorkspace(w)
# Exploit convention encoded in named set "parameters" and "observables"
# to use workspace contents w/o need for introspected
model = w.pdf("model")
# Generate data from p.d.f. in given observables
data = model.generate(w.set("observables"), 1000)
# Fit model to data
model.fitTo(data)
# Plot fitted model and data on frame of first (only) observable
frame = (w.set("observables").first()).frame()
data.plotOn(frame)
model.plotOn(frame)
# Overlay plot with model with reference parameters as stored in snapshots
w.loadSnapshot("reference_fit")
model.plotOn(frame, ROOT.RooFit.LineColor(ROOT.kRed))
w.loadSnapshot("reference_fit_bkgonly")
model.plotOn(frame, ROOT.RooFit.LineColor(ROOT.kRed),
ROOT.RooFit.LineStyle(ROOT.kDashed))
# Draw the frame on the canvas
c = ROOT.TCanvas("rf510_wsnamedsets", "rf503_wsnamedsets", 600, 600)
ROOT.gPad.SetLeftMargin(0.15)
frame.GetYaxis().SetTitleOffset(1.4)
frame.Draw()
c.SaveAs("rf510_wsnamedsets.png")
# Print workspace contents
w.Print()
# Workspace will remain in memory after macro finishes
ROOT.gDirectory.Add(w)
def fillWorkspace(w):
# C r e a t e m o d e l
# -----------------------
# 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, 0, 10)
sigma1 = ROOT.RooRealVar("sigma1", "width of gaussians", 0.5)
sigma2 = ROOT.RooRealVar("sigma2", "width of gaussians", 1)
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, 0., 1.)
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))
# Import model into p.d.f.
getattr(w, 'import')(model)
# E n c o d e d e f i n i t i o n o f p a r a m e t e r s i n w o r k s p a c e
# ---------------------------------------------------------------------------------------
# Define named sets "parameters" and "observables", list which variables should be considered
# parameters and observables by the users convention
#
# Variables appearing in sets _must_ live in the workspace already, the autoImport flag
# of defineSet must be set to import them on the fly. Named sets contain only references
# to the original variables, the value of observables in named sets already
# reflect their 'current' value
params = model.getParameters(ROOT.RooArgSet(x))
w.defineSet("parameters", params)
w.defineSet("observables", ROOT.RooArgSet(x))
# E n c o d e r e f e r e n c e v a l u e f o r p a r a m e t e r s i n w o r k s p a c e
# ---------------------------------------------------------------------------------------------------
# Define a parameter 'snapshot' in the p.d.f.
# Unlike a named set, parameter snapshot stores an independent set of values for
# a given set of variables in the workspace. ROOT.The values can be stored and reloaded
# into the workspace variable objects using the loadSnapshot() and saveSnapshot()
# methods. A snapshot saves the value of each variable, errors that are stored
# with it as well as the 'Constant' flag that is used in fits to determine if a
# parameter is kept fixed or not.
# Do a dummy fit to a (supposedly) reference dataset here and store the results
# of that fit into a snapshot
refData = model.generate(ROOT.RooArgSet(x), 10000)
model.fitTo(refData, ROOT.RooFit.PrintLevel(-1))
# ROOT.The kTRUE flag imports the values of the objects in (*params) into the workspace
# If not set, present values of the workspace parameters objects are stored
w.saveSnapshot("reference_fit", params, ROOT.kTRUE)
# Make another fit with the signal componentforced to zero
# and save those parameters too
bkgfrac.setVal(1)
bkgfrac.setConstant(ROOT.kTRUE)
bkgfrac.removeError()
model.fitTo(refData, ROOT.RooFit.PrintLevel(-1))
w.saveSnapshot("reference_fit_bkgonly", params, ROOT.kTRUE)
if __name__ == "__main__":
rf510_wsnamedsets()