-
Notifications
You must be signed in to change notification settings - Fork 18
/
rf603_multicpu.py
99 lines (77 loc) · 3.47 KB
/
rf603_multicpu.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#####################################
#
# 'LIKELIHOOD AND MINIMIZATION' ROOT.RooFit tutorial macro #603
#
# Setting up a multi-core parallelized unbinned maximum likelihood fit
#
#
#
# 07/2008 - Wouter Verkerke
#
# /
import ROOT
def rf603_multicpu():
# C r e a t e 3 D p d f a n d d a t a
# -------------------------------------------
# Create observables
x = ROOT.RooRealVar("x", "x", -5, 5)
y = ROOT.RooRealVar("y", "y", -5, 5)
z = ROOT.RooRealVar("z", "z", -5, 5)
# Create signal pdf gauss(x)*gauss(y)*gauss(z)
gx = ROOT.RooGaussian(
"gx", "gx", x, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
gy = ROOT.RooGaussian(
"gy", "gy", y, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
gz = ROOT.RooGaussian(
"gz", "gz", z, ROOT.RooFit.RooConst(0), ROOT.RooFit.RooConst(1))
sig = ROOT.RooProdPdf("sig", "sig", ROOT.RooArgList(gx, gy, gz))
# Create background pdf poly(x)*poly(y)*poly(z)
px = ROOT.RooPolynomial("px", "px", x, ROOT.RooArgList(
ROOT.RooFit.RooConst(-0.1), ROOT.RooFit.RooConst(0.004)))
py = ROOT.RooPolynomial("py", "py", y, ROOT.RooArgList(
ROOT.RooFit.RooConst(0.1), ROOT.RooFit.RooConst(-0.004)))
pz = ROOT.RooPolynomial("pz", "pz", z)
bkg = ROOT.RooProdPdf("bkg", "bkg", ROOT.RooArgList(px, py, pz))
# Create composite pdf sig+bkg
fsig = ROOT.RooRealVar("fsig", "signal fraction", 0.1, 0., 1.)
model = ROOT.RooAddPdf("model", "model", ROOT.RooArgList(
sig, bkg), ROOT.RooArgList(fsig))
# Generate large dataset
data = model.generate(ROOT.RooArgSet(x, y, z), 200000)
# P a r a l l e l f i t t i n g
# -------------------------------
# In parallel mode the likelihood calculation is split in N pieces,
# that are calculated in parallel and added a posteriori before passing
# it back to MINUIT.
# Use four processes and time results both in wall time and CPU time
model.fitTo(data, ROOT.RooFit.NumCPU(4), ROOT.RooFit.Timer(ROOT.kTRUE))
# P a r a l l e l M C p r o j e c t i o n s
# ----------------------------------------------
# Construct signal, likelihood projection on (y,z) observables and
# likelihood ratio
sigyz = sig.createProjection(ROOT.RooArgSet(x))
totyz = model.createProjection(ROOT.RooArgSet(x))
llratio_func = ROOT.RooFormulaVar(
"llratio", "log10(@0)-log10(@1)", ROOT.RooArgList(sigyz, totyz))
# Calculate likelihood ratio for each event, subset of events with high
# signal likelihood
data.addColumn(llratio_func)
dataSel = data.reduce(ROOT.RooFit.Cut("llratio>0.7"))
# Make plot frame and plot data
frame = x.frame(ROOT.RooFit.Title(
"Projection on X with LLratio(y,z)>0.7"), ROOT.RooFit.Bins(40))
dataSel.plotOn(frame)
# Perform parallel projection using MC integration of pdf using given input dataSet.
# In self mode the data-weighted average of the pdf is calculated by splitting the
# input dataset in N equal pieces and calculating in parallel the weighted average
# one each subset. ROOT.The N results of those calculations are then weighted into the
# final result
# Use four processes
model.plotOn(frame, ROOT.RooFit.ProjWData(dataSel), ROOT.RooFit.NumCPU(4))
c = ROOT.TCanvas("rf603_multicpu", "rf603_multicpu", 600, 600)
ROOT.gPad.SetLeftMargin(0.15)
frame.GetYaxis().SetTitleOffset(1.6)
frame.Draw()
c.SaveAs("rf603_multicpu.png")
if __name__ == "__main__":
rf603_multicpu()