Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
rf501_simultaneouspdf.py
Go to the documentation of this file.
1 ## \file
2 ## \ingroup tutorial_roofit
3 ## \notebook
4 ##
5 ## Organization and simultaneous fits: using simultaneous p.d.f.s to describe simultaneous fits to multiple datasets
6 ##
7 ## \macro_code
8 ##
9 ## \date February 2018
10 ## \author Clemens Lange, Wouter Verkerke (C++ version)
11 
12 import ROOT
13 
14 
15 # Create model for physics sample
16 # -------------------------------------------------------------
17 
18 # Create observables
19 x = ROOT.RooRealVar("x", "x", -8, 8)
20 
21 # Construct signal pdf
22 mean = ROOT.RooRealVar("mean", "mean", 0, -8, 8)
23 sigma = ROOT.RooRealVar("sigma", "sigma", 0.3, 0.1, 10)
24 gx = ROOT.RooGaussian("gx", "gx", x, mean, sigma)
25 
26 # Construct background pdf
27 a0 = ROOT.RooRealVar("a0", "a0", -0.1, -1, 1)
28 a1 = ROOT.RooRealVar("a1", "a1", 0.004, -1, 1)
29 px = ROOT.RooChebychev("px", "px", x, ROOT.RooArgList(a0, a1))
30 
31 # Construct composite pdf
32 f = ROOT.RooRealVar("f", "f", 0.2, 0., 1.)
33 model = 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
41 mean_ctl = ROOT.RooRealVar("mean_ctl", "mean_ctl", -3, -8, 8)
42 gx_ctl = ROOT.RooGaussian("gx_ctl", "gx_ctl", x, mean_ctl, sigma)
43 
44 # Construct the background pdf
45 a0_ctl = ROOT.RooRealVar("a0_ctl", "a0_ctl", -0.1, -1, 1)
46 a1_ctl = ROOT.RooRealVar("a1_ctl", "a1_ctl", 0.5, -0.1, 1)
47 px_ctl = ROOT.RooChebychev(
48  "px_ctl", "px_ctl", x, ROOT.RooArgList(a0_ctl, a1_ctl))
49 
50 # Construct the composite model
51 f_ctl = ROOT.RooRealVar("f_ctl", "f_ctl", 0.5, 0., 1.)
52 model_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
64 data = model.generate(ROOT.RooArgSet(x), 100)
65 data_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
71 sample = ROOT.RooCategory("sample", "sample")
72 sample.defineType("physics")
73 sample.defineType("control")
74 
75 # Construct combined dataset in (x,sample)
76 combData = 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
92 simPdf = ROOT.RooSimultaneous("simPdf", "simultaneous pdf", sample)
93 
94 # Associate model with the physics state and model_ctl with the control
95 # state
96 simPdf.addPdf(model, "physics")
97 simPdf.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
103 simPdf.fitTo(combData)
104 
105 # Plot model slices on data slices
106 # ----------------------------------------------------------------
107 
108 # Make a frame for the physics sample
109 frame1 = x.frame(ROOT.RooFit.Bins(30), ROOT.RooFit.Title("Physics sample"))
110 
111 # Plot all data tagged as physics sample
112 combData.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.
120 sampleSet = ROOT.RooArgSet(sample)
121 simPdf.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
125 frame2 = x.frame(ROOT.RooFit.Bins(30), ROOT.RooFit.Title("Control sample"))
126 combData.plotOn(frame2, ROOT.RooFit.Cut("sample==sample::control"))
127 simPdf.plotOn(frame2, ROOT.RooFit.Slice(sample, "control"),
128  ROOT.RooFit.ProjWData(sampleSet, combData))
129 simPdf.plotOn(frame2, ROOT.RooFit.Slice(sample, "control"), ROOT.RooFit.Components(
130  "px_ctl"), ROOT.RooFit.ProjWData(sampleSet, combData), ROOT.RooFit.LineStyle(ROOT.kDashed))
131 
132 c = ROOT.TCanvas("rf501_simultaneouspdf",
133  "rf501_simultaneouspdf", 800, 400)
134 c.Divide(2)
135 c.cd(1)
136 ROOT.gPad.SetLeftMargin(0.15)
137 frame1.GetYaxis().SetTitleOffset(1.4)
138 frame1.Draw()
139 c.cd(2)
140 ROOT.gPad.SetLeftMargin(0.15)
141 frame2.GetYaxis().SetTitleOffset(1.4)
142 frame2.Draw()
143 
144 c.SaveAs("rf501_simultaneouspdf.png")