Logo ROOT  
Reference Guide
rf315_projectpdf.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Multidimensional models: marginizalization of multi-dimensional pdfs through integration
5##
6## \macro_code
7##
8## \date February 2018
9## \authors Clemens Lange, Wouter Verkerke (C++ version)
10
11import ROOT
12
13
14# Create pdf m(x,y) = gx(x|y) * g(y)
15# --------------------------------------------------------------
16
17# Increase default precision of numeric integration
18# as self exercise has high sensitivity to numeric integration precision
19ROOT.RooAbsPdf.defaultIntegratorConfig().setEpsRel(1e-8)
20ROOT.RooAbsPdf.defaultIntegratorConfig().setEpsAbs(1e-8)
21
22# Create observables
23x = ROOT.RooRealVar("x", "x", -5, 5)
24y = ROOT.RooRealVar("y", "y", -2, 2)
25
26# Create function f(y) = a0 + a1*y
27a0 = ROOT.RooRealVar("a0", "a0", 0)
28a1 = ROOT.RooRealVar("a1", "a1", -1.5, -3, 1)
29fy = ROOT.RooPolyVar("fy", "fy", y, ROOT.RooArgList(a0, a1))
30
31# Create gaussx(x,f(y),sx)
32sigmax = ROOT.RooRealVar("sigmax", "width of gaussian", 0.5)
33gaussx = ROOT.RooGaussian(
34 "gaussx", "Gaussian in x with shifting mean in y", x, fy, sigmax)
35
36# Create gaussy(y,0,2)
37gaussy = ROOT.RooGaussian(
38 "gaussy",
39 "Gaussian in y",
40 y,
41 ROOT.RooFit.RooConst(0),
42 ROOT.RooFit.RooConst(2))
43
44# Create gaussx(x,sx|y) * gaussy(y)
45model = ROOT.RooProdPdf(
46 "model",
47 "gaussx(x|y)*gaussy(y)",
48 ROOT.RooArgSet(gaussy),
49 ROOT.RooFit.Conditional(
50 ROOT.RooArgSet(gaussx),
51 ROOT.RooArgSet(x)))
52
53# Marginalize m(x,y) to m(x)
54# ----------------------------------------------------
55
56# modelx(x) = Int model(x,y) dy
57modelx = model.createProjection(ROOT.RooArgSet(y))
58
59# Use marginalized pdf as regular 1D pdf
60# -----------------------------------------------
61
62# Sample 1000 events from modelx
63data = modelx.generateBinned(ROOT.RooArgSet(x), 1000)
64
65# Fit modelx to toy data
66modelx.fitTo(data, ROOT.RooFit.Verbose())
67
68# Plot modelx over data
69frame = x.frame(40)
70data.plotOn(frame)
71modelx.plotOn(frame)
72
73# Make 2D histogram of model(x,y)
74hh = model.createHistogram("x,y")
75hh.SetLineColor(ROOT.kBlue)
76
77c = ROOT.TCanvas("rf315_projectpdf", "rf315_projectpdf", 800, 400)
78c.Divide(2)
79c.cd(1)
80ROOT.gPad.SetLeftMargin(0.15)
81frame.GetYaxis().SetTitleOffset(1.4)
82frame.Draw()
83c.cd(2)
84ROOT.gPad.SetLeftMargin(0.20)
85hh.GetZaxis().SetTitleOffset(2.5)
86hh.Draw("surf")
87c.SaveAs("rf315_projectpdf.png")