ROOT   Reference Guide
Searching...
No Matches
rf801_mcstudy.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Validation and MC studies: toy Monte Carlo study that perform cycles of event generation and fitting
5##
6## \macro_code
7##
8## \date February 2018
9## \authors Clemens Lange, Wouter Verkerke (C++ version)
10
11import ROOT
12
13
14# Create model
15# -----------------------
16
17# Declare observable x
18x = ROOT.RooRealVar("x", "x", 0, 10)
19x.setBins(40)
20
21# Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and
22# their parameters
23mean = ROOT.RooRealVar("mean", "mean of gaussians", 5, 0, 10)
24sigma1 = ROOT.RooRealVar("sigma1", "width of gaussians", 0.5)
25sigma2 = ROOT.RooRealVar("sigma2", "width of gaussians", 1)
26
27sig1 = ROOT.RooGaussian("sig1", "Signal component 1", x, mean, sigma1)
28sig2 = ROOT.RooGaussian("sig2", "Signal component 2", x, mean, sigma2)
29
30# Build Chebychev polynomial pdf
31a0 = ROOT.RooRealVar("a0", "a0", 0.5, 0., 1.)
32a1 = ROOT.RooRealVar("a1", "a1", -0.2, -1, 1.)
33bkg = ROOT.RooChebychev("bkg", "Background", x, ROOT.RooArgList(a0, a1))
34
35# Sum the signal components into a composite signal pdf
36sig1frac = ROOT.RooRealVar(
37 "sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.)
39 "sig", "Signal", ROOT.RooArgList(sig1, sig2), ROOT.RooArgList(sig1frac))
40
41# Sum the composite signal and background
42nbkg = ROOT.RooRealVar(
43 "nbkg", "number of background events, ", 150, 0, 1000)
44nsig = ROOT.RooRealVar("nsig", "number of signal events", 150, 0, 1000)
46 "model", "g1+g2+a", ROOT.RooArgList(bkg, sig), ROOT.RooArgList(nbkg, nsig))
47
48# Create manager
49# ---------------------------
50
51# Instantiate ROOT.RooMCStudy manager on model with x as observable and given choice of fit options
52#
53# The Silence() option kills all messages below the PROGRESS level, only a single message
54# per sample executed, any error message that occur during fitting
55#
56# The Extended() option has two effects:
57# 1) The extended ML term is included in the likelihood and
58# 2) A poisson fluctuation is introduced on the number of generated events
59#
60# The FitOptions() given here are passed to the fitting stage of each toy experiment.
61# If Save() is specified, fit result of each experiment is saved by the manager
62#
63# A Binned() option is added in self example to bin the data between generation and fitting
64# to speed up the study at the expemse of some precision
65
66mcstudy = ROOT.RooMCStudy(
67 model,
68 ROOT.RooArgSet(x),
69 ROOT.RooFit.Binned(
70 ROOT.kTRUE),
71 ROOT.RooFit.Silence(),
72 ROOT.RooFit.Extended(),
73 ROOT.RooFit.FitOptions(
74 ROOT.RooFit.Save(
75 ROOT.kTRUE),
76 ROOT.RooFit.PrintEvalErrors(0)))
77
78# Generate and fit events
79# ---------------------------------------------
80
81# Generate and fit 1000 samples of Poisson(nExpected) events
82mcstudy.generateAndFit(1000)
83
84# Explore results of study
85# ------------------------------------------------
86
87# Make plots of the distributions of mean, error on mean and the pull of
88# mean
89frame1 = mcstudy.plotParam(mean, ROOT.RooFit.Bins(40))
90frame2 = mcstudy.plotError(mean, ROOT.RooFit.Bins(40))
91frame3 = mcstudy.plotPull(mean, ROOT.RooFit.Bins(
92 40), ROOT.RooFit.FitGauss(ROOT.kTRUE))
93
94# Plot distribution of minimized likelihood
95frame4 = mcstudy.plotNLL(ROOT.RooFit.Bins(40))
96
97# Make some histograms from the parameter dataset
98hh_cor_a0_s1f = ROOT.RooAbsData.createHistogram(
99 mcstudy.fitParDataSet(), "hh", a1, ROOT.RooFit.YVar(sig1frac))
100hh_cor_a0_a1 = ROOT.RooAbsData.createHistogram(mcstudy.fitParDataSet(),
101 "hh", a0, ROOT.RooFit.YVar(a1))
102
103# Access some of the saved fit results from individual toys
104corrHist000 = mcstudy.fitResult(0).correlationHist("c000")
105corrHist127 = mcstudy.fitResult(127).correlationHist("c127")
106corrHist953 = mcstudy.fitResult(953).correlationHist("c953")
107
108# Draw all plots on a canvas
109ROOT.gStyle.SetPalette(1)
110ROOT.gStyle.SetOptStat(0)
111c = ROOT.TCanvas("rf801_mcstudy", "rf801_mcstudy", 900, 900)
112c.Divide(3, 3)
113c.cd(1)
115frame1.GetYaxis().SetTitleOffset(1.4)
116frame1.Draw()
117c.cd(2)
119frame2.GetYaxis().SetTitleOffset(1.4)
120frame2.Draw()
121c.cd(3)
123frame3.GetYaxis().SetTitleOffset(1.4)
124frame3.Draw()
125c.cd(4)
127frame4.GetYaxis().SetTitleOffset(1.4)
128frame4.Draw()
129c.cd(5)
131hh_cor_a0_s1f.GetYaxis().SetTitleOffset(1.4)
132hh_cor_a0_s1f.Draw("box")
133c.cd(6)
135hh_cor_a0_a1.GetYaxis().SetTitleOffset(1.4)
136hh_cor_a0_a1.Draw("box")
137c.cd(7)
139corrHist000.GetYaxis().SetTitleOffset(1.4)
140corrHist000.Draw("colz")
141c.cd(8)