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