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 
11 import 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
19 ROOT.RooAbsPdf.defaultIntegratorConfig().setEpsRel(1e-8)
20 ROOT.RooAbsPdf.defaultIntegratorConfig().setEpsAbs(1e-8)
21 
22 # Create observables
23 x = ROOT.RooRealVar("x", "x", -5, 5)
24 y = ROOT.RooRealVar("y", "y", -2, 2)
25 
26 # Create function f(y) = a0 + a1*y
27 a0 = ROOT.RooRealVar("a0", "a0", 0)
28 a1 = ROOT.RooRealVar("a1", "a1", -1.5, -3, 1)
29 fy = ROOT.RooPolyVar("fy", "fy", y, ROOT.RooArgList(a0, a1))
30 
31 # Create gaussx(x,f(y),sx)
32 sigmax = ROOT.RooRealVar("sigmax", "width of gaussian", 0.5)
33 gaussx = ROOT.RooGaussian(
34  "gaussx", "Gaussian in x with shifting mean in y", x, fy, sigmax)
35 
36 # Create gaussy(y,0,2)
37 gaussy = 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)
45 model = 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
57 modelx = model.createProjection(ROOT.RooArgSet(y))
58 
59 # Use marginalized pdf as regular 1D pdf
60 # -----------------------------------------------
61 
62 # Sample 1000 events from modelx
63 data = modelx.generateBinned(ROOT.RooArgSet(x), 1000)
64 
65 # Fit modelx to toy data
66 modelx.fitTo(data, ROOT.RooFit.Verbose())
67 
68 # Plot modelx over data
69 frame = x.frame(40)
70 data.plotOn(frame)
71 modelx.plotOn(frame)
72 
73 # Make 2D histogram of model(x,y)
74 hh = model.createHistogram("x,y")
75 hh.SetLineColor(ROOT.kBlue)
76 
77 c = ROOT.TCanvas("rf315_projectpdf", "rf315_projectpdf", 800, 400)
78 c.Divide(2)
79 c.cd(1)
80 ROOT.gPad.SetLeftMargin(0.15)
81 frame.GetYaxis().SetTitleOffset(1.4)
82 frame.Draw()
83 c.cd(2)
84 ROOT.gPad.SetLeftMargin(0.20)
85 hh.GetZaxis().SetTitleOffset(2.5)
86 hh.Draw("surf")
87 c.SaveAs("rf315_projectpdf.png")