Special pdf's: unbinned maximum likelihood fit of an efficiency eff(x) function to a dataset D(x,cut), cut is a category encoding a selection whose efficiency as function of x should be described by eff(x)
Author: Clemens Lange, Wouter Verkerke (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, April 17, 2024 at 11:19 AM.
import ROOT
flat = False
Declare variables x,mean, with associated name, title, value and allowed range
x = ROOT.RooRealVar("x", "x", -10, 10)
y = ROOT.RooRealVar("y", "y", -10, 10)
Efficiency function eff(x;a,b)
ax = ROOT.RooRealVar("ax", "ay", 0.6, 0, 1)
bx = ROOT.RooRealVar("bx", "by", 5)
cx = ROOT.RooRealVar("cx", "cy", -1, -10, 10)
ay = ROOT.RooRealVar("ay", "ay", 0.2, 0, 1)
by = ROOT.RooRealVar("by", "by", 5)
cy = ROOT.RooRealVar("cy", "cy", -1, -10, 10)
effFunc = ROOT.RooFormulaVar(
"effFunc", "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))", [ax, bx, cx, x, ay, by, cy, y]
)
Acceptance state cut (1 or 0)
cut = ROOT.RooCategory("cut", "cutr", {"accept": 1, "reject": 0})
Construct efficiency pdf eff(cut|x)
effPdf = ROOT.RooEfficiency("effPdf", "effPdf", effFunc, cut, "accept")
Construct global shape pdf shape(x) and product model(x,cut) = eff(cut|x)*shape(x) (These are only needed to generate some toy MC here to be used later)
shapePdfX = ROOT.RooPolynomial("shapePdfX", "shapePdfX", x, [0 if flat else -0.095])
shapePdfY = ROOT.RooPolynomial("shapePdfY", "shapePdfY", y, [0 if flat else +0.095])
shapePdf = ROOT.RooProdPdf("shapePdf", "shapePdf", [shapePdfX, shapePdfY])
model = ROOT.RooProdPdf("model", "model", {shapePdf}, Conditional=({effPdf}, {cut}))
Generate some toy data from model
data = model.generate({x, y, cut}, 10000)
[#0] WARNING:Generation -- RooAcceptReject::ctor(effPdf_Int[]_Norm[cut]) WARNING: performing accept/reject sampling on a p.d.f in 2 dimensions without prior knowledge on maximum value of p.d.f. Determining maximum value by taking 200000 trial samples. If p.d.f contains sharp peaks smaller than average distance between trial sampling points these may be missed and p.d.f. may be sampled incorrectly. [#0] WARNING:Generation -- RooAcceptReject::ctor(effPdf_Int[]_Norm[cut]): WARNING: 200000 trial samples requested by p.d.f for 2-dimensional accept/reject sampling, this may take some time
Fit conditional efficiency pdf to data
effPdf.fitTo(data, ConditionalObservables={x, y}, PrintLevel=-1)
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:Fitting -- RooAbsPdf::fitTo(effPdf_over_effPdf_Int[cut]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- using CPU computation library compiled with -mavx2 [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_effPdf_over_effPdf_Int[cut]_modelData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
Make 2D histograms of all data, data and efficiency function
hh_data_all = data.createHistogram("hh_data_all", x, Binning=8, YVar=dict(var=y, Binning=8))
hh_data_sel = data.createHistogram("hh_data_sel", x, Binning=8, YVar=dict(var=y, Binning=8), Cut="cut==cut::accept")
hh_eff = effFunc.createHistogram("hh_eff", x, Binning=50, YVar=dict(var=y, Binning=50))
Some adjustsment for good visualization
hh_data_all.SetMinimum(0)
hh_data_sel.SetMinimum(0)
hh_eff.SetMinimum(0)
hh_eff.SetLineColor(ROOT.kBlue)
Draw all frames on a canvas
ca = ROOT.TCanvas("rf702_efficiency_2D", "rf702_efficiency_2D", 1200, 400)
ca.Divide(3)
ca.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
hh_data_all.GetZaxis().SetTitleOffset(1.8)
hh_data_all.Draw("lego")
ca.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
hh_data_sel.GetZaxis().SetTitleOffset(1.8)
hh_data_sel.Draw("lego")
ca.cd(3)
ROOT.gPad.SetLeftMargin(0.15)
hh_eff.GetZaxis().SetTitleOffset(1.8)
hh_eff.Draw("surf")
ca.SaveAs("rf702_efficiency_2D.png")
Info in <TCanvas::Print>: png file rf702_efficiency_2D.png has been created
Draw all canvases
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()