Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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_image
9## \macro_code
10## \macro_output
11##
12## \date February 2018
13## \authors Clemens Lange, Wouter Verkerke (C++ version)
14
15import ROOT
16
17
18flat = False
19# Construct efficiency function e(x,y)
20# -----------------------------------------------------------------------
21
22# Declare variables x,mean, with associated name, title, value and allowed
23# range
24x = ROOT.RooRealVar("x", "x", -10, 10)
25y = ROOT.RooRealVar("y", "y", -10, 10)
26
27# Efficiency function eff(x;a,b)
28ax = ROOT.RooRealVar("ax", "ay", 0.6, 0, 1)
29bx = ROOT.RooRealVar("bx", "by", 5)
30cx = ROOT.RooRealVar("cx", "cy", -1, -10, 10)
31
32ay = ROOT.RooRealVar("ay", "ay", 0.2, 0, 1)
33by = ROOT.RooRealVar("by", "by", 5)
34cy = ROOT.RooRealVar("cy", "cy", -1, -10, 10)
35
36effFunc = ROOT.RooFormulaVar(
37 "effFunc", "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))", [ax, bx, cx, x, ay, by, cy, y]
38)
39
40# Acceptance state cut (1 or 0)
41cut = ROOT.RooCategory("cut", "cutr", {"accept": 1, "reject": 0})
42
43# Construct conditional efficiency pdf E(cut|x,y)
44# ---------------------------------------------------------------------------------------------
45
46# Construct efficiency pdf eff(cut|x)
47effPdf = ROOT.RooEfficiency("effPdf", "effPdf", effFunc, cut, "accept")
48
49# Generate data(x,y,cut) from a toy model
50# -------------------------------------------------------------------------------
51
52# Construct global shape pdf shape(x) and product model(x,cut) = eff(cut|x)*shape(x)
53# (These are _only_ needed to generate some toy MC here to be used later)
54shapePdfX = ROOT.RooPolynomial("shapePdfX", "shapePdfX", x, [0 if flat else -0.095])
55shapePdfY = ROOT.RooPolynomial("shapePdfY", "shapePdfY", y, [0 if flat else +0.095])
56shapePdf = ROOT.RooProdPdf("shapePdf", "shapePdf", [shapePdfX, shapePdfY])
57model = ROOT.RooProdPdf("model", "model", {shapePdf}, Conditional=({effPdf}, {cut}))
58
59# Generate some toy data from model
60data = model.generate({x, y, cut}, 10000)
61
62# Fit conditional efficiency pdf to data
63# --------------------------------------------------------------------------
64
65# Fit conditional efficiency pdf to data
66effPdf.fitTo(data, ConditionalObservables={x, y}, PrintLevel=-1)
67
68# Plot fitted, data efficiency
69# --------------------------------------------------------
70
71# Make 2D histograms of all data, data and efficiency function
72hh_data_all = data.createHistogram("hh_data_all", x, Binning=8, YVar=dict(var=y, Binning=8))
73hh_data_sel = data.createHistogram("hh_data_sel", x, Binning=8, YVar=dict(var=y, Binning=8), Cut="cut==cut::accept")
74hh_eff = effFunc.createHistogram("hh_eff", x, Binning=50, YVar=dict(var=y, Binning=50))
75
76# Some adjustsment for good visualization
77hh_data_all.SetMinimum(0)
78hh_data_sel.SetMinimum(0)
79hh_eff.SetMinimum(0)
80hh_eff.SetLineColor(ROOT.kBlue)
81
82# Draw all frames on a canvas
83ca = ROOT.TCanvas("rf702_efficiency_2D", "rf702_efficiency_2D", 1200, 400)
84ca.Divide(3)
85ca.cd(1)
86ROOT.gPad.SetLeftMargin(0.15)
87hh_data_all.GetZaxis().SetTitleOffset(1.8)
88hh_data_all.Draw("lego")
89ca.cd(2)
90ROOT.gPad.SetLeftMargin(0.15)
91hh_data_sel.GetZaxis().SetTitleOffset(1.8)
92hh_data_sel.Draw("lego")
93ca.cd(3)
94ROOT.gPad.SetLeftMargin(0.15)
95hh_eff.GetZaxis().SetTitleOffset(1.8)
96hh_eff.Draw("surf")
97
98ca.SaveAs("rf702_efficiency_2D.png")