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