Logo ROOT  
Reference Guide
rf316_llratioplot.py
Go to the documentation of this file.
1 ## \file
2 ## \ingroup tutorial_roofit
3 ## \notebook
4 ##
5 ## Multidimensional models: using the likelihood ratio techique to construct a signal
6 ## enhanced one-dimensional projection of a multi-dimensional p.d.f.
7 ##
8 ## \macro_code
9 ##
10 ## \date February 2018
11 ## \authors Clemens Lange, Wouter Verkerke (C++ version)
12 
13 import ROOT
14 
15 
16 # Create 3D pdf and data
17 # -------------------------------------------
18 
19 # Create observables
20 x = ROOT.RooRealVar("x", "x", -5, 5)
21 y = ROOT.RooRealVar("y", "y", -5, 5)
22 z = ROOT.RooRealVar("z", "z", -5, 5)
23 
24 # Create signal pdf gauss(x)*gauss(y)*gauss(z)
25 gx = ROOT.RooGaussian(
26  "gx", "gx", x, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
27 gy = ROOT.RooGaussian(
28  "gy", "gy", y, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
29 gz = ROOT.RooGaussian(
30  "gz", "gz", z, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
31 sig = ROOT.RooProdPdf("sig", "sig", ROOT.RooArgList(gx, gy, gz))
32 
33 # Create background pdf poly(x)*poly(y)*poly(z)
34 px = ROOT.RooPolynomial("px", "px", x, ROOT.RooArgList(
35  ROOT.RooFit.RooConst(-0.1), ROOT.RooFit.RooConst(0.004)))
36 py = ROOT.RooPolynomial("py", "py", y, ROOT.RooArgList(
37  ROOT.RooFit.RooConst(0.1), ROOT.RooFit.RooConst(-0.004)))
38 pz = ROOT.RooPolynomial("pz", "pz", z)
39 bkg = ROOT.RooProdPdf("bkg", "bkg", ROOT.RooArgList(px, py, pz))
40 
41 # Create composite pdf sig+bkg
42 fsig = ROOT.RooRealVar("fsig", "signal fraction", 0.1, 0., 1.)
43 model = ROOT.RooAddPdf(
44  "model", "model", ROOT.RooArgList(
45  sig, bkg), ROOT.RooArgList(fsig))
46 
47 data = model.generate(ROOT.RooArgSet(x, y, z), 20000)
48 
49 # Project pdf and data on x
50 # -------------------------------------------------
51 
52 # Make plain projection of data and pdf on x observable
53 frame = x.frame(ROOT.RooFit.Title(
54  "Projection of 3D data and pdf on X"), ROOT.RooFit.Bins(40))
55 data.plotOn(frame)
56 model.plotOn(frame)
57 
58 # Define projected signal likelihood ratio
59 # ----------------------------------------------------------------------------------
60 
61 # Calculate projection of signal and total likelihood on (y,z) observables
62 # i.e. integrate signal and composite model over x
63 sigyz = sig.createProjection(ROOT.RooArgSet(x))
64 totyz = model.createProjection(ROOT.RooArgSet(x))
65 
66 # Construct the log of the signal / signal+background probability
67 llratio_func = ROOT.RooFormulaVar(
68  "llratio", "log10(@0)-log10(@1)", ROOT.RooArgList(sigyz, totyz))
69 
70 # Plot data with a LL ratio cut
71 # -------------------------------------------------------
72 
73 # Calculate the llratio value for each event in the dataset
74 data.addColumn(llratio_func)
75 
76 # Extract the subset of data with large signal likelihood
77 dataSel = data.reduce(ROOT.RooFit.Cut("llratio>0.7"))
78 
79 # Make plot frame
80 frame2 = x.frame(ROOT.RooFit.Title(
81  "Same projection on X with LLratio(y,z)>0.7"), ROOT.RooFit.Bins(40))
82 
83 # Plot select data on frame
84 dataSel.plotOn(frame2)
85 
86 # Make MC projection of pdf with same LL ratio cut
87 # ---------------------------------------------------------------------------------------------
88 
89 # Generate large number of events for MC integration of pdf projection
90 mcprojData = model.generate(ROOT.RooArgSet(x, y, z), 10000)
91 
92 # Calculate LL ratio for each generated event and select MC events with
93 # llratio)0.7
94 mcprojData.addColumn(llratio_func)
95 mcprojDataSel = mcprojData.reduce(ROOT.RooFit.Cut("llratio>0.7"))
96 
97 # Project model on x, projected observables (y,z) with Monte Carlo technique
98 # on set of events with the same llratio cut as was applied to data
99 model.plotOn(frame2, ROOT.RooFit.ProjWData(mcprojDataSel))
100 
101 c = ROOT.TCanvas("rf316_llratioplot", "rf316_llratioplot", 800, 400)
102 c.Divide(2)
103 c.cd(1)
104 ROOT.gPad.SetLeftMargin(0.15)
105 frame.GetYaxis().SetTitleOffset(1.4)
106 frame.Draw()
107 c.cd(2)
108 ROOT.gPad.SetLeftMargin(0.15)
109 frame2.GetYaxis().SetTitleOffset(1.4)
110 frame2.Draw()
111 c.SaveAs("rf316_llratioplot.png")