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