Logo ROOT  
Reference Guide
rf603_multicpu.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Likelihood and minimization: setting up a multi-core parallelized unbinned maximum likelihood fit
5##
6## \macro_code
7##
8## \date February 2018
9## \authors Clemens Lange, Wouter Verkerke (C++ version)
10
11import ROOT
12
13
14# Create 3D pdf and data
15# -------------------------------------------
16
17# Create observables
18x = ROOT.RooRealVar("x", "x", -5, 5)
19y = ROOT.RooRealVar("y", "y", -5, 5)
20z = ROOT.RooRealVar("z", "z", -5, 5)
21
22# Create signal pdf gauss(x)*gauss(y)*gauss(z)
23gx = ROOT.RooGaussian(
24 "gx", "gx", x, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
25gy = ROOT.RooGaussian(
26 "gy", "gy", y, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
27gz = ROOT.RooGaussian(
28 "gz", "gz", z, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
29sig = ROOT.RooProdPdf("sig", "sig", ROOT.RooArgList(gx, gy, gz))
30
31# Create background pdf poly(x)*poly(y)*poly(z)
32px = ROOT.RooPolynomial("px", "px", x, ROOT.RooArgList(
33 ROOT.RooFit.RooConst(-0.1), ROOT.RooFit.RooConst(0.004)))
34py = ROOT.RooPolynomial("py", "py", y, ROOT.RooArgList(
35 ROOT.RooFit.RooConst(0.1), ROOT.RooFit.RooConst(-0.004)))
36pz = ROOT.RooPolynomial("pz", "pz", z)
37bkg = ROOT.RooProdPdf("bkg", "bkg", ROOT.RooArgList(px, py, pz))
38
39# Create composite pdf sig+bkg
40fsig = ROOT.RooRealVar("fsig", "signal fraction", 0.1, 0., 1.)
41model = ROOT.RooAddPdf("model", "model", ROOT.RooArgList(
42 sig, bkg), ROOT.RooArgList(fsig))
43
44# Generate large dataset
45data = model.generate(ROOT.RooArgSet(x, y, z), 200000)
46
47# Parallel fitting
48# -------------------------------
49
50# In parallel mode the likelihood calculation is split in N pieces,
51# that are calculated in parallel and added a posteriori before passing
52# it back to MINUIT.
53
54# Use four processes and time results both in wall time and CPU time
55model.fitTo(data, ROOT.RooFit.NumCPU(4), ROOT.RooFit.Timer(ROOT.kTRUE))
56
57# Parallel MC projections
58# ----------------------------------------------
59
60# Construct signal, likelihood projection on (y,z) observables and
61# likelihood ratio
62sigyz = sig.createProjection(ROOT.RooArgSet(x))
63totyz = model.createProjection(ROOT.RooArgSet(x))
64llratio_func = ROOT.RooFormulaVar(
65 "llratio", "log10(@0)-log10(@1)", ROOT.RooArgList(sigyz, totyz))
66
67# Calculate likelihood ratio for each event, subset of events with high
68# signal likelihood
69data.addColumn(llratio_func)
70dataSel = data.reduce(ROOT.RooFit.Cut("llratio>0.7"))
71
72# Make plot frame and plot data
73frame = x.frame(ROOT.RooFit.Title(
74 "Projection on X with LLratio(y,z)>0.7"), ROOT.RooFit.Bins(40))
75dataSel.plotOn(frame)
76
77# Perform parallel projection using MC integration of pdf using given input dataSet.
78# In self mode the data-weighted average of the pdf is calculated by splitting the
79# input dataset in N equal pieces and calculating in parallel the weighted average
80# one each subset. The N results of those calculations are then weighted into the
81# final result
82
83# Use four processes
84model.plotOn(frame, ROOT.RooFit.ProjWData(dataSel), ROOT.RooFit.NumCPU(4))
85
86c = ROOT.TCanvas("rf603_multicpu", "rf603_multicpu", 600, 600)
87ROOT.gPad.SetLeftMargin(0.15)
88frame.GetYaxis().SetTitleOffset(1.6)
89frame.Draw()
90
91c.SaveAs("rf603_multicpu.png")