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 
14 import ROOT
15 
16 
17 flat = 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
23 x = ROOT.RooRealVar("x", "x", -10, 10)
24 y = ROOT.RooRealVar("y", "y", -10, 10)
25 
26 # Efficiency function eff(x;a,b)
27 ax = ROOT.RooRealVar("ax", "ay", 0.6, 0, 1)
28 bx = ROOT.RooRealVar("bx", "by", 5)
29 cx = ROOT.RooRealVar("cx", "cy", -1, -10, 10)
30 
31 ay = ROOT.RooRealVar("ay", "ay", 0.2, 0, 1)
32 by = ROOT.RooRealVar("by", "by", 5)
33 cy = ROOT.RooRealVar("cy", "cy", -1, -10, 10)
34 
35 effFunc = 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)
49 cut = ROOT.RooCategory("cut", "cutr")
50 cut.defineType("accept", 1)
51 cut.defineType("reject", 0)
52 
53 # Construct conditional efficiency pdf E(cut|x,y)
54 # ---------------------------------------------------------------------------------------------
55 
56 # Construct efficiency p.d.f eff(cut|x)
57 effPdf = 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)
64 shapePdfX = ROOT.RooPolynomial(
65  "shapePdfX", "shapePdfX", x, ROOT.RooArgList(
66  ROOT.RooFit.RooConst(
67  0 if flat else -0.095)))
68 shapePdfY = ROOT.RooPolynomial(
69  "shapePdfY", "shapePdfY", y, ROOT.RooArgList(
70  ROOT.RooFit.RooConst(
71  0 if flat else +0.095)))
72 shapePdf = ROOT.RooProdPdf(
73  "shapePdf",
74  "shapePdf",
75  ROOT.RooArgList(
76  shapePdfX,
77  shapePdfY))
78 model = 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
87 data = 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
93 effPdf.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
99 hh_data_all = ROOT.RooAbsData.createHistogram(
100  data, "hh_data_all", x, ROOT.RooFit.Binning(8), ROOT.RooFit.YVar(
101  y, ROOT.RooFit.Binning(8)))
102 hh_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"))
105 hh_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
110 hh_data_all.SetMinimum(0)
111 hh_data_sel.SetMinimum(0)
112 hh_eff.SetMinimum(0)
113 hh_eff.SetLineColor(ROOT.kBlue)
114 
115 # Draw all frames on a canvas
116 ca = ROOT.TCanvas("rf702_efficiency_2D", "rf702_efficiency_2D", 1200, 400)
117 ca.Divide(3)
118 ca.cd(1)
119 ROOT.gPad.SetLeftMargin(0.15)
120 hh_data_all.GetZaxis().SetTitleOffset(1.8)
121 hh_data_all.Draw("lego")
122 ca.cd(2)
123 ROOT.gPad.SetLeftMargin(0.15)
124 hh_data_sel.GetZaxis().SetTitleOffset(1.8)
125 hh_data_sel.Draw("lego")
126 ca.cd(3)
127 ROOT.gPad.SetLeftMargin(0.15)
128 hh_eff.GetZaxis().SetTitleOffset(1.8)
129 hh_eff.Draw("surf")
130 
131 ca.SaveAs("rf702_efficiency_2D.png")