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