Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rf608_fitresultaspdf.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Likelihood and minimization: representing the parabolic approximation of the fit as a
5## multi-variate Gaussian on the parameters of the fitted pdf
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 and dataset
16# -----------------------------------------------
17
18# Observable
19x = ROOT.RooRealVar("x", "x", -20, 20)
20
21# Model (intentional strong correlations)
22mean = ROOT.RooRealVar("mean", "mean of g1 and g2", 0, -1, 1)
23sigma_g1 = ROOT.RooRealVar("sigma_g1", "width of g1", 2)
24g1 = ROOT.RooGaussian("g1", "g1", x, mean, sigma_g1)
25
26sigma_g2 = ROOT.RooRealVar("sigma_g2", "width of g2", 4, 3.0, 5.0)
27g2 = ROOT.RooGaussian("g2", "g2", x, mean, sigma_g2)
28
29frac = ROOT.RooRealVar("frac", "frac", 0.5, 0.0, 1.0)
30model = ROOT.RooAddPdf(
31 "model", "model", ROOT.RooArgList(
32 g1, g2), ROOT.RooArgList(frac))
33
34# Generate 1000 events
35data = model.generate(ROOT.RooArgSet(x), 1000)
36
37# Fit model to data
38# ----------------------------------
39
40r = model.fitTo(data, ROOT.RooFit.Save())
41
42# Create MV Gaussian pdf of fitted parameters
43# ------------------------------------------------------------------------------------
44
45parabPdf = r.createHessePdf(ROOT.RooArgSet(frac, mean, sigma_g2))
46
47# Some exercises with the parameter pdf
48# -----------------------------------------------------------------------------
49
50# Generate 100K points in the parameter space, from the MVGaussian pdf
51d = parabPdf.generate(ROOT.RooArgSet(mean, sigma_g2, frac), 100000)
52
53# Sample a 3-D histogram of the pdf to be visualized as an error
54# ellipsoid using the GLISO draw option
55hh_3d = parabPdf.createHistogram("mean,sigma_g2,frac", 25, 25, 25)
56hh_3d.SetFillColor(ROOT.kBlue)
57
58# Project 3D parameter pdf down to 3 permutations of two-dimensional pdfs
59# The integrations corresponding to these projections are performed analytically
60# by the MV Gaussian pdf
61pdf_sigmag2_frac = parabPdf.createProjection(ROOT.RooArgSet(mean))
62pdf_mean_frac = parabPdf.createProjection(ROOT.RooArgSet(sigma_g2))
63pdf_mean_sigmag2 = parabPdf.createProjection(ROOT.RooArgSet(frac))
64
65# Make 2D plots of the 3 two-dimensional pdf projections
66hh_sigmag2_frac = pdf_sigmag2_frac.createHistogram("sigma_g2,frac", 50, 50)
67hh_mean_frac = pdf_mean_frac.createHistogram("mean,frac", 50, 50)
68hh_mean_sigmag2 = pdf_mean_sigmag2.createHistogram("mean,sigma_g2", 50, 50)
69hh_mean_frac.SetLineColor(ROOT.kBlue)
70hh_sigmag2_frac.SetLineColor(ROOT.kBlue)
71hh_mean_sigmag2.SetLineColor(ROOT.kBlue)
72
73# Draw the 'sigar'
74ROOT.gStyle.SetCanvasPreferGL(True)
75ROOT.gStyle.SetPalette(1)
76c1 = ROOT.TCanvas("rf608_fitresultaspdf_1", "rf608_fitresultaspdf_1", 600, 600)
77hh_3d.Draw("gliso")
78
79c1.SaveAs("rf608_fitresultaspdf_1.png")
80
81# Draw the 2D projections of the 3D pdf
82c2 = ROOT.TCanvas("rf608_fitresultaspdf_2",
83 "rf608_fitresultaspdf_2", 900, 600)
84c2.Divide(3, 2)
85c2.cd(1)
86ROOT.gPad.SetLeftMargin(0.15)
87hh_mean_sigmag2.GetZaxis().SetTitleOffset(1.4)
88hh_mean_sigmag2.Draw("surf3")
89c2.cd(2)
90ROOT.gPad.SetLeftMargin(0.15)
91hh_sigmag2_frac.GetZaxis().SetTitleOffset(1.4)
92hh_sigmag2_frac.Draw("surf3")
93c2.cd(3)
94ROOT.gPad.SetLeftMargin(0.15)
95hh_mean_frac.GetZaxis().SetTitleOffset(1.4)
96hh_mean_frac.Draw("surf3")
97
98# Draw the distributions of parameter points sampled from the pdf
99tmp1 = d.createHistogram(mean, sigma_g2, 50, 50)
100tmp2 = d.createHistogram(sigma_g2, frac, 50, 50)
101tmp3 = d.createHistogram(mean, frac, 50, 50)
102
103c2.cd(4)
104ROOT.gPad.SetLeftMargin(0.15)
105tmp1.GetZaxis().SetTitleOffset(1.4)
106tmp1.Draw("lego3")
107c2.cd(5)
108ROOT.gPad.SetLeftMargin(0.15)
109tmp2.GetZaxis().SetTitleOffset(1.4)
110tmp2.Draw("lego3")
111c2.cd(6)
112ROOT.gPad.SetLeftMargin(0.15)
113tmp3.GetZaxis().SetTitleOffset(1.4)
114tmp3.Draw("lego3")
115
116c2.SaveAs("rf608_fitresultaspdf_2.png")