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