Logo ROOT  
Reference Guide
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 
12 import ROOT
13 
14 
15 # Create model and dataset
16 # -----------------------------------------------
17 
18 # Observable
19 x = ROOT.RooRealVar("x", "x", -20, 20)
20 
21 # Model (intentional strong correlations)
22 mean = ROOT.RooRealVar("mean", "mean of g1 and g2", 0, -1, 1)
23 sigma_g1 = ROOT.RooRealVar("sigma_g1", "width of g1", 2)
24 g1 = ROOT.RooGaussian("g1", "g1", x, mean, sigma_g1)
25 
26 sigma_g2 = ROOT.RooRealVar("sigma_g2", "width of g2", 4, 3.0, 5.0)
27 g2 = ROOT.RooGaussian("g2", "g2", x, mean, sigma_g2)
28 
29 frac = ROOT.RooRealVar("frac", "frac", 0.5, 0.0, 1.0)
30 model = ROOT.RooAddPdf(
31  "model", "model", ROOT.RooArgList(
32  g1, g2), ROOT.RooArgList(frac))
33 
34 # Generate 1000 events
35 data = model.generate(ROOT.RooArgSet(x), 1000)
36 
37 # Fit model to data
38 # ----------------------------------
39 
40 r = model.fitTo(data, ROOT.RooFit.Save())
41 
42 # Create MV Gaussian pdf of fitted parameters
43 # ------------------------------------------------------------------------------------
44 
45 parabPdf = 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
51 d = 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
55 hh_3d = parabPdf.createHistogram("mean,sigma_g2,frac", 25, 25, 25)
56 hh_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
61 pdf_sigmag2_frac = parabPdf.createProjection(ROOT.RooArgSet(mean))
62 pdf_mean_frac = parabPdf.createProjection(ROOT.RooArgSet(sigma_g2))
63 pdf_mean_sigmag2 = parabPdf.createProjection(ROOT.RooArgSet(frac))
64 
65 # Make 2D plots of the 3 two-dimensional pdf projections
66 hh_sigmag2_frac = pdf_sigmag2_frac.createHistogram("sigma_g2,frac", 50, 50)
67 hh_mean_frac = pdf_mean_frac.createHistogram("mean,frac", 50, 50)
68 hh_mean_sigmag2 = pdf_mean_sigmag2.createHistogram("mean,sigma_g2", 50, 50)
69 hh_mean_frac.SetLineColor(ROOT.kBlue)
70 hh_sigmag2_frac.SetLineColor(ROOT.kBlue)
71 hh_mean_sigmag2.SetLineColor(ROOT.kBlue)
72 
73 # Draw the 'sigar'
74 ROOT.gStyle.SetCanvasPreferGL(True)
75 ROOT.gStyle.SetPalette(1)
76 c1 = ROOT.TCanvas("rf608_fitresultaspdf_1", "rf608_fitresultaspdf_1", 600, 600)
77 hh_3d.Draw("gliso")
78 
79 c1.SaveAs("rf608_fitresultaspdf_1.png")
80 
81 # Draw the 2D projections of the 3D pdf
82 c2 = ROOT.TCanvas("rf608_fitresultaspdf_2",
83  "rf608_fitresultaspdf_2", 900, 600)
84 c2.Divide(3, 2)
85 c2.cd(1)
86 ROOT.gPad.SetLeftMargin(0.15)
87 hh_mean_sigmag2.GetZaxis().SetTitleOffset(1.4)
88 hh_mean_sigmag2.Draw("surf3")
89 c2.cd(2)
90 ROOT.gPad.SetLeftMargin(0.15)
91 hh_sigmag2_frac.GetZaxis().SetTitleOffset(1.4)
92 hh_sigmag2_frac.Draw("surf3")
93 c2.cd(3)
94 ROOT.gPad.SetLeftMargin(0.15)
95 hh_mean_frac.GetZaxis().SetTitleOffset(1.4)
96 hh_mean_frac.Draw("surf3")
97 
98 # Draw the distributions of parameter points sampled from the pdf
99 tmp1 = d.createHistogram(mean, sigma_g2, 50, 50)
100 tmp2 = d.createHistogram(sigma_g2, frac, 50, 50)
101 tmp3 = d.createHistogram(mean, frac, 50, 50)
102 
103 c2.cd(4)
104 ROOT.gPad.SetLeftMargin(0.15)
105 tmp1.GetZaxis().SetTitleOffset(1.4)
106 tmp1.Draw("lego3")
107 c2.cd(5)
108 ROOT.gPad.SetLeftMargin(0.15)
109 tmp2.GetZaxis().SetTitleOffset(1.4)
110 tmp2.Draw("lego3")
111 c2.cd(6)
112 ROOT.gPad.SetLeftMargin(0.15)
113 tmp3.GetZaxis().SetTitleOffset(1.4)
114 tmp3.Draw("lego3")
115 
116 c2.SaveAs("rf608_fitresultaspdf_2.png")