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