-
Notifications
You must be signed in to change notification settings - Fork 18
/
Copy pathrf801_mcstudy.py
151 lines (126 loc) · 5.34 KB
/
rf801_mcstudy.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
142
143
144
145
146
147
148
149
150
151
# /
#
# 'VALIDATION AND MC STUDIES' ROOT.RooFit tutorial macro #801
#
# A ROOT.Toy Monte Carlo study that perform cycles of
# event generation and fittting
#
#
# /
import ROOT
def rf801_mcstudy():
# C r e a t e m o d e l
# -----------------------
# Declare observable x
x = ROOT.RooRealVar("x", "x", 0, 10)
x.setBins(40)
# 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, -1, 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
nbkg = ROOT.RooRealVar(
"nbkg", "number of background events, ", 150, 0, 1000)
nsig = ROOT.RooRealVar("nsig", "number of signal events", 150, 0, 1000)
model = ROOT.RooAddPdf(
"model", "g1+g2+a", ROOT.RooArgList(bkg, sig), ROOT.RooArgList(nbkg, nsig))
# C r e a t e m a n a g e r
# ---------------------------
# Instantiate ROOT.RooMCStudy manager on model with x as observable and given choice of fit options
#
# ROOT.The Silence() option kills all messages below the PROGRESS level, only a single message
# per sample executed, any error message that occur during fitting
#
# ROOT.The Extended() option has two effects:
# 1) ROOT.The extended ML term is included in the likelihood and
# 2) A poisson fluctuation is introduced on the number of generated events
#
# ROOT.The FitOptions() given here are passed to the fitting stage of each toy experiment.
# If Save() is specified, fit result of each experiment is saved by the manager
#
# A Binned() option is added in self example to bin the data between generation and fitting
# to speed up the study at the expemse of some precision
mcstudy = ROOT.RooMCStudy(model, ROOT.RooArgSet(x), ROOT.RooFit.Binned(ROOT.kTRUE), ROOT.RooFit.Silence(), ROOT.RooFit.Extended(),
ROOT.RooFit.FitOptions(ROOT.RooFit.Save(ROOT.kTRUE), ROOT.RooFit.PrintEvalErrors(0)))
# G e n e r a t e a n d f i t e v e n t s
# ---------------------------------------------
# Generate and fit 1000 samples of Poisson(nExpected) events
mcstudy.generateAndFit(1000)
# E x p l o r e r e s u l t s o f s t u d y
# ------------------------------------------------
# Make plots of the distributions of mean, error on mean and the pull of
# mean
frame1 = mcstudy.plotParam(mean, ROOT.RooFit.Bins(40))
frame2 = mcstudy.plotError(mean, ROOT.RooFit.Bins(40))
frame3 = mcstudy.plotPull(mean, ROOT.RooFit.Bins(
40), ROOT.RooFit.FitGauss(ROOT.kTRUE))
# Plot distribution of minimized likelihood
frame4 = mcstudy.plotNLL(ROOT.RooFit.Bins(40))
# Make some histograms from the parameter dataset
hh_cor_a0_s1f = ROOT.RooAbsData.createHistogram(mcstudy.fitParDataSet(),
"hh", a1, ROOT.RooFit.YVar(sig1frac))
hh_cor_a0_a1 = ROOT.RooAbsData.createHistogram(mcstudy.fitParDataSet(),
"hh", a0, ROOT.RooFit.YVar(a1))
# Access some of the saved fit results from individual toys
corrHist000 = mcstudy.fitResult(0).correlationHist("c000")
corrHist127 = mcstudy.fitResult(127).correlationHist("c127")
corrHist953 = mcstudy.fitResult(953).correlationHist("c953")
# Draw all plots on a canvas
ROOT.gStyle.SetPalette(1)
ROOT.gStyle.SetOptStat(0)
c = ROOT.TCanvas("rf801_mcstudy", "rf801_mcstudy", 900, 900)
c.Divide(3, 3)
c.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
frame1.GetYaxis().SetTitleOffset(1.4)
frame1.Draw()
c.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
frame2.GetYaxis().SetTitleOffset(1.4)
frame2.Draw()
c.cd(3)
ROOT.gPad.SetLeftMargin(0.15)
frame3.GetYaxis().SetTitleOffset(1.4)
frame3.Draw()
c.cd(4)
ROOT.gPad.SetLeftMargin(0.15)
frame4.GetYaxis().SetTitleOffset(1.4)
frame4.Draw()
c.cd(5)
ROOT.gPad.SetLeftMargin(0.15)
hh_cor_a0_s1f.GetYaxis().SetTitleOffset(1.4)
hh_cor_a0_s1f.Draw("box")
c.cd(6)
ROOT.gPad.SetLeftMargin(0.15)
hh_cor_a0_a1.GetYaxis().SetTitleOffset(1.4)
hh_cor_a0_a1.Draw("box")
c.cd(7)
ROOT.gPad.SetLeftMargin(0.15)
corrHist000.GetYaxis().SetTitleOffset(1.4)
corrHist000.Draw("colz")
c.cd(8)
ROOT.gPad.SetLeftMargin(0.15)
corrHist127.GetYaxis().SetTitleOffset(1.4)
corrHist127.Draw("colz")
c.cd(9)
ROOT.gPad.SetLeftMargin(0.15)
corrHist953.GetYaxis().SetTitleOffset(1.4)
corrHist953.Draw("colz")
c.SaveAs("rf801_mcstudy.png")
# Make ROOT.RooMCStudy object available on command line after
# macro finishes
ROOT.gDirectory.Add(mcstudy)
if __name__ == "__main__":
rf801_mcstudy()