Logo ROOT  
Reference Guide
rf501_simultaneouspdf.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Organization and simultaneous fits: using simultaneous pdfs to describe simultaneous
5## fits to multiple datasets
6##
7## \macro_code
8##
9## \date February 2018
10## \authors Clemens Lange, Wouter Verkerke (C++ version)
11
12import ROOT
13
14
15# Create model for physics sample
16# -------------------------------------------------------------
17
18# Create observables
19x = ROOT.RooRealVar("x", "x", -8, 8)
20
21# Construct signal pdf
22mean = ROOT.RooRealVar("mean", "mean", 0, -8, 8)
23sigma = ROOT.RooRealVar("sigma", "sigma", 0.3, 0.1, 10)
24gx = ROOT.RooGaussian("gx", "gx", x, mean, sigma)
25
26# Construct background pdf
27a0 = ROOT.RooRealVar("a0", "a0", -0.1, -1, 1)
28a1 = ROOT.RooRealVar("a1", "a1", 0.004, -1, 1)
29px = ROOT.RooChebychev("px", "px", x, [a0, a1])
30
31# Construct composite pdf
32f = ROOT.RooRealVar("f", "f", 0.2, 0.0, 1.0)
33model = ROOT.RooAddPdf("model", "model", [gx, px], [f])
34
35# Create model for control sample
36# --------------------------------------------------------------
37
38# Construct signal pdf.
39# NOTE that sigma is shared with the signal sample model
40mean_ctl = ROOT.RooRealVar("mean_ctl", "mean_ctl", -3, -8, 8)
41gx_ctl = ROOT.RooGaussian("gx_ctl", "gx_ctl", x, mean_ctl, sigma)
42
43# Construct the background pdf
44a0_ctl = ROOT.RooRealVar("a0_ctl", "a0_ctl", -0.1, -1, 1)
45a1_ctl = ROOT.RooRealVar("a1_ctl", "a1_ctl", 0.5, -0.1, 1)
46px_ctl = ROOT.RooChebychev("px_ctl", "px_ctl", x, [a0_ctl, a1_ctl])
47
48# Construct the composite model
49f_ctl = ROOT.RooRealVar("f_ctl", "f_ctl", 0.5, 0.0, 1.0)
50model_ctl = ROOT.RooAddPdf("model_ctl", "model_ctl", [gx_ctl, px_ctl], [f_ctl])
51
52# Generate events for both samples
53# ---------------------------------------------------------------
54
55# Generate 1000 events in x and y from model
56data = model.generate({x}, 100)
57data_ctl = model_ctl.generate({x}, 2000)
58
59# Create index category and join samples
60# ---------------------------------------------------------------------------
61
62# Define category to distinguish physics and control samples events
63sample = ROOT.RooCategory("sample", "sample")
64sample.defineType("physics")
65sample.defineType("control")
66
67# Construct combined dataset in (x,sample)
68combData = ROOT.RooDataSet(
69 "combData",
70 "combined data",
71 {x},
72 Index=sample,
73 Import={"physics": data, "control": data_ctl},
74)
75
76# Construct a simultaneous pdf in (x, sample)
77# -----------------------------------------------------------------------------------
78
79# Construct a simultaneous pdf using category sample as index: associate model
80# with the physics state and model_ctl with the control state
81simPdf = ROOT.RooSimultaneous("simPdf", "simultaneous pdf", {"physics": model, "control": model_ctl}, sample)
82
83# Perform a simultaneous fit
84# ---------------------------------------------------
85
86# Perform simultaneous fit of model to data and model_ctl to data_ctl
87fitResult = simPdf.fitTo(combData, PrintLevel=-1, Save=True)
88fitResult.Print()
89
90# Plot model slices on data slices
91# ----------------------------------------------------------------
92
93# Make a frame for the physics sample
94frame1 = x.frame(Bins=30, Title="Physics sample")
95
96# Plot all data tagged as physics sample
97combData.plotOn(frame1, Cut="sample==sample::physics")
98
99# Plot "physics" slice of simultaneous pdf.
100# NB: You *must* project the sample index category with data using ProjWData as
101# a RooSimultaneous makes no prediction on the shape in the index category and
102# can thus not be integrated. In other words: Since the PDF doesn't know the
103# number of events in the different category states, it doesn't know how much
104# of each component it has to project out. This info is read from the data.
105simPdf.plotOn(frame1, Slice=(sample, "physics"), ProjWData=(sample, combData))
106simPdf.plotOn(frame1, Slice=(sample, "physics"), Components="px", ProjWData=(sample, combData), LineStyle="--")
107
108# The same plot for the control sample slice
109frame2 = x.frame(Bins=30, Title="Control sample")
110combData.plotOn(frame2, Cut="sample==sample::control")
111simPdf.plotOn(frame2, Slice=(sample, "control"), ProjWData=(sample, combData))
112simPdf.plotOn(frame2, Slice=(sample, "control"), Components="px_ctl", ProjWData=(sample, combData), LineStyle="--")
113
114c = ROOT.TCanvas("rf501_simultaneouspdf", "rf501_simultaneouspdf", 800, 400)
115c.Divide(2)
116c.cd(1)
117ROOT.gPad.SetLeftMargin(0.15)
118frame1.GetYaxis().SetTitleOffset(1.4)
119frame1.Draw()
120c.cd(2)
121ROOT.gPad.SetLeftMargin(0.15)
122frame2.GetYaxis().SetTitleOffset(1.4)
123frame2.Draw()
124
125c.SaveAs("rf501_simultaneouspdf.png")