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