Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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, ROOT.RooArgList(a0, a1))
30
31# Construct composite pdf
32f = ROOT.RooRealVar("f", "f", 0.2, 0., 1.)
33model = ROOT.RooAddPdf(
34 "model", "model", ROOT.RooArgList(gx, px), ROOT.RooArgList(f))
35
36# Create model for control sample
37# --------------------------------------------------------------
38
39# Construct signal pdf.
40# NOTE that sigma is shared with the signal sample model
41mean_ctl = ROOT.RooRealVar("mean_ctl", "mean_ctl", -3, -8, 8)
42gx_ctl = ROOT.RooGaussian("gx_ctl", "gx_ctl", x, mean_ctl, sigma)
43
44# Construct the background pdf
45a0_ctl = ROOT.RooRealVar("a0_ctl", "a0_ctl", -0.1, -1, 1)
46a1_ctl = ROOT.RooRealVar("a1_ctl", "a1_ctl", 0.5, -0.1, 1)
47px_ctl = ROOT.RooChebychev(
48 "px_ctl", "px_ctl", x, ROOT.RooArgList(a0_ctl, a1_ctl))
49
50# Construct the composite model
51f_ctl = ROOT.RooRealVar("f_ctl", "f_ctl", 0.5, 0., 1.)
52model_ctl = ROOT.RooAddPdf(
53 "model_ctl",
54 "model_ctl",
55 ROOT.RooArgList(
56 gx_ctl,
57 px_ctl),
58 ROOT.RooArgList(f_ctl))
59
60# Generate events for both samples
61# ---------------------------------------------------------------
62
63# Generate 1000 events in x and y from model
64data = model.generate(ROOT.RooArgSet(x), 100)
65data_ctl = model_ctl.generate(ROOT.RooArgSet(x), 2000)
66
67# Create index category and join samples
68# ---------------------------------------------------------------------------
69
70# Define category to distinguish physics and control samples events
71sample = ROOT.RooCategory("sample", "sample")
72sample.defineType("physics")
73sample.defineType("control")
74
75# Construct combined dataset in (x,sample)
76combData = ROOT.RooDataSet(
77 "combData",
78 "combined data",
79 ROOT.RooArgSet(x),
80 ROOT.RooFit.Index(sample),
81 ROOT.RooFit.Import(
82 "physics",
83 data),
84 ROOT.RooFit.Import(
85 "control",
86 data_ctl))
87
88# Construct a simultaneous pdf in (x, sample)
89# -----------------------------------------------------------------------------------
90
91# Construct a simultaneous pdf using category sample as index
92simPdf = ROOT.RooSimultaneous("simPdf", "simultaneous pdf", sample)
93
94# Associate model with the physics state and model_ctl with the control
95# state
96simPdf.addPdf(model, "physics")
97simPdf.addPdf(model_ctl, "control")
98
99# Perform a simultaneous fit
100# ---------------------------------------------------
101
102# Perform simultaneous fit of model to data and model_ctl to data_ctl
103simPdf.fitTo(combData)
104
105# Plot model slices on data slices
106# ----------------------------------------------------------------
107
108# Make a frame for the physics sample
109frame1 = x.frame(ROOT.RooFit.Bins(30), ROOT.RooFit.Title("Physics sample"))
110
111# Plot all data tagged as physics sample
112combData.plotOn(frame1, ROOT.RooFit.Cut("sample==sample::physics"))
113
114# Plot "physics" slice of simultaneous pdf.
115# NB: You *must* project the sample index category with data using ProjWData
116# as a RooSimultaneous makes no prediction on the shape in the index category
117# and can thus not be integrated
118# NB2: The sampleSet *must* be named. It will not work to pass this as a temporary
119# because python will delete it. The same holds for fitTo() and plotOn() below.
120sampleSet = ROOT.RooArgSet(sample)
121simPdf.plotOn(frame1, ROOT.RooFit.Slice(sample, "physics"), ROOT.RooFit.Components(
122 "px"), ROOT.RooFit.ProjWData(sampleSet, combData), ROOT.RooFit.LineStyle(ROOT.kDashed))
123
124# The same plot for the control sample slice
125frame2 = x.frame(ROOT.RooFit.Bins(30), ROOT.RooFit.Title("Control sample"))
126combData.plotOn(frame2, ROOT.RooFit.Cut("sample==sample::control"))
127simPdf.plotOn(frame2, ROOT.RooFit.Slice(sample, "control"),
128 ROOT.RooFit.ProjWData(sampleSet, combData))
129simPdf.plotOn(frame2, ROOT.RooFit.Slice(sample, "control"), ROOT.RooFit.Components(
130 "px_ctl"), ROOT.RooFit.ProjWData(sampleSet, combData), ROOT.RooFit.LineStyle(ROOT.kDashed))
131
132c = ROOT.TCanvas("rf501_simultaneouspdf",
133 "rf501_simultaneouspdf", 800, 400)
134c.Divide(2)
135c.cd(1)
136ROOT.gPad.SetLeftMargin(0.15)
137frame1.GetYaxis().SetTitleOffset(1.4)
138frame1.Draw()
139c.cd(2)
140ROOT.gPad.SetLeftMargin(0.15)
141frame2.GetYaxis().SetTitleOffset(1.4)
142frame2.Draw()
143
144c.SaveAs("rf501_simultaneouspdf.png")