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