ROOT   Reference Guide
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_image
8## \macro_code
9## \macro_output
10##
11## \date February 2018
12## \authors Clemens Lange, Wouter Verkerke (C++ version)
13
14import ROOT
15
16
17# Create model and dataset
18# -----------------------------------------------
19
20# Observable
21x = ROOT.RooRealVar("x", "x", -20, 20)
22
23# Model (intentional strong correlations)
24mean = ROOT.RooRealVar("mean", "mean of g1 and g2", 0, -1, 1)
25sigma_g1 = ROOT.RooRealVar("sigma_g1", "width of g1", 2)
26g1 = ROOT.RooGaussian("g1", "g1", x, mean, sigma_g1)
27
28sigma_g2 = ROOT.RooRealVar("sigma_g2", "width of g2", 4, 3.0, 5.0)
29g2 = ROOT.RooGaussian("g2", "g2", x, mean, sigma_g2)
30
31frac = ROOT.RooRealVar("frac", "frac", 0.5, 0.0, 1.0)
32model = ROOT.RooAddPdf("model", "model", [g1, g2], [frac])
33
34# Generate 1000 events
35data = model.generate({x}, 1000)
36
37# Fit model to data
38# ----------------------------------
39
40r = model.fitTo(data, Save=True, PrintLevel=-1)
41
42# Create MV Gaussian pdf of fitted parameters
43# ------------------------------------------------------------------------------------
44
45parabPdf = r.createHessePdf({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({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({mean})
62pdf_mean_frac = parabPdf.createProjection({sigma_g2})
63pdf_mean_sigmag2 = parabPdf.createProjection({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", "rf608_fitresultaspdf_2", 900, 600)
83c2.Divide(3, 2)
84c2.cd(1)
86hh_mean_sigmag2.GetZaxis().SetTitleOffset(1.4)
87hh_mean_sigmag2.Draw("surf3")
88c2.cd(2)
90hh_sigmag2_frac.GetZaxis().SetTitleOffset(1.4)
91hh_sigmag2_frac.Draw("surf3")
92c2.cd(3)
94hh_mean_frac.GetZaxis().SetTitleOffset(1.4)
95hh_mean_frac.Draw("surf3")
96
97# Draw the distributions of parameter points sampled from the pdf
98tmp1 = d.createHistogram(mean, sigma_g2, 50, 50)
99tmp2 = d.createHistogram(sigma_g2, frac, 50, 50)
100tmp3 = d.createHistogram(mean, frac, 50, 50)
101
102c2.cd(4)