Logo ROOT  
Reference Guide
rf702_efficiencyfit_2D.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Special pdf's: unbinned maximum likelihood fit of an efficiency eff(x) function
5## to a dataset D(x,cut), cut is a category encoding a selection whose efficiency as function
6## of x should be described by eff(x)
7##
8## \macro_code
9##
10## \date February 2018
11## \authors Clemens Lange, Wouter Verkerke (C++ version)
12
13import ROOT
14
15
16flat = ROOT.kFALSE
17# Construct efficiency function e(x,y)
18# -----------------------------------------------------------------------
19
20# Declare variables x,mean, with associated name, title, value and allowed
21# range
22x = ROOT.RooRealVar("x", "x", -10, 10)
23y = ROOT.RooRealVar("y", "y", -10, 10)
24
25# Efficiency function eff(x;a,b)
26ax = ROOT.RooRealVar("ax", "ay", 0.6, 0, 1)
27bx = ROOT.RooRealVar("bx", "by", 5)
28cx = ROOT.RooRealVar("cx", "cy", -1, -10, 10)
29
30ay = ROOT.RooRealVar("ay", "ay", 0.2, 0, 1)
31by = ROOT.RooRealVar("by", "by", 5)
32cy = ROOT.RooRealVar("cy", "cy", -1, -10, 10)
33
34effFunc = ROOT.RooFormulaVar(
35 "effFunc",
36 "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))",
37 ROOT.RooArgList(
38 ax,
39 bx,
40 cx,
41 x,
42 ay,
43 by,
44 cy,
45 y))
46
47# Acceptance state cut (1 or 0)
48cut = ROOT.RooCategory("cut", "cutr")
49cut.defineType("accept", 1)
50cut.defineType("reject", 0)
51
52# Construct conditional efficiency pdf E(cut|x,y)
53# ---------------------------------------------------------------------------------------------
54
55# Construct efficiency pdf eff(cut|x)
56effPdf = ROOT.RooEfficiency("effPdf", "effPdf", effFunc, cut, "accept")
57
58# Generate data(x,y,cut) from a toy model
59# -------------------------------------------------------------------------------
60
61# Construct global shape pdf shape(x) and product model(x,cut) = eff(cut|x)*shape(x)
62# (These are _only_ needed to generate some toy MC here to be used later)
63shapePdfX = ROOT.RooPolynomial(
64 "shapePdfX", "shapePdfX", x, ROOT.RooArgList(
65 ROOT.RooFit.RooConst(
66 0 if flat else -0.095)))
67shapePdfY = ROOT.RooPolynomial(
68 "shapePdfY", "shapePdfY", y, ROOT.RooArgList(
69 ROOT.RooFit.RooConst(
70 0 if flat else +0.095)))
71shapePdf = ROOT.RooProdPdf(
72 "shapePdf",
73 "shapePdf",
74 ROOT.RooArgList(
75 shapePdfX,
76 shapePdfY))
77model = ROOT.RooProdPdf(
78 "model",
79 "model",
80 ROOT.RooArgSet(shapePdf),
81 ROOT.RooFit.Conditional(
82 ROOT.RooArgSet(effPdf),
83 ROOT.RooArgSet(cut)))
84
85# Generate some toy data from model
86data = model.generate(ROOT.RooArgSet(x, y, cut), 10000)
87
88# Fit conditional efficiency pdf to data
89# --------------------------------------------------------------------------
90
91# Fit conditional efficiency pdf to data
92effPdf.fitTo(data, ROOT.RooFit.ConditionalObservables(ROOT.RooArgSet(x, y)))
93
94# Plot fitted, data efficiency
95# --------------------------------------------------------
96
97# Make 2D histograms of all data, data and efficiency function
98hh_data_all = ROOT.RooAbsData.createHistogram(
99 data, "hh_data_all", x, ROOT.RooFit.Binning(8), ROOT.RooFit.YVar(
100 y, ROOT.RooFit.Binning(8)))
101hh_data_sel = ROOT.RooAbsData.createHistogram(
102 data, "hh_data_sel", x, ROOT.RooFit.Binning(8), ROOT.RooFit.YVar(
103 y, ROOT.RooFit.Binning(8)), ROOT.RooFit.Cut("cut==cut::accept"))
104hh_eff = effFunc.createHistogram(
105 "hh_eff", x, ROOT.RooFit.Binning(
106 50), ROOT.RooFit.YVar(y, ROOT.RooFit.Binning(50)))
107
108# Some adjustsment for good visualization
109hh_data_all.SetMinimum(0)
110hh_data_sel.SetMinimum(0)
111hh_eff.SetMinimum(0)
112hh_eff.SetLineColor(ROOT.kBlue)
113
114# Draw all frames on a canvas
115ca = ROOT.TCanvas("rf702_efficiency_2D", "rf702_efficiency_2D", 1200, 400)
116ca.Divide(3)
117ca.cd(1)
118ROOT.gPad.SetLeftMargin(0.15)
119hh_data_all.GetZaxis().SetTitleOffset(1.8)
120hh_data_all.Draw("lego")
121ca.cd(2)
122ROOT.gPad.SetLeftMargin(0.15)
123hh_data_sel.GetZaxis().SetTitleOffset(1.8)
124hh_data_sel.Draw("lego")
125ca.cd(3)
126ROOT.gPad.SetLeftMargin(0.15)
127hh_eff.GetZaxis().SetTitleOffset(1.8)
128hh_eff.Draw("surf")
129
130ca.SaveAs("rf702_efficiency_2D.png")