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