rf501_simultaneouspdf¶

Organization and simultaneous fits: using simultaneous pdfs to describe simultaneous fits to multiple datasets

Author: Clemens Lange, Wouter Verkerke (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Monday, February 06, 2023 at 11:18 AM.

In [1]:

import ROOT

Welcome to JupyROOT 6.29/01

Create model for physics sample¶

Create observables

In [2]:

x = ROOT.RooRealVar("x", "x", -8, 8)

Construct signal pdf

In [3]:

mean = ROOT.RooRealVar("mean", "mean", 0, -8, 8)
sigma = ROOT.RooRealVar("sigma", "sigma", 0.3, 0.1, 10)
gx = ROOT.RooGaussian("gx", "gx", x, mean, sigma)

Construct background pdf

In [4]:

a0 = ROOT.RooRealVar("a0", "a0", -0.1, -1, 1)
a1 = ROOT.RooRealVar("a1", "a1", 0.004, -1, 1)
px = ROOT.RooChebychev("px", "px", x, [a0, a1])

Construct composite pdf

In [5]:

f = ROOT.RooRealVar("f", "f", 0.2, 0.0, 1.0)
model = ROOT.RooAddPdf("model", "model", [gx, px], [f])

Create model for control sample¶

Construct signal pdf. NOTE that sigma is shared with the signal sample model

In [6]:

mean_ctl = ROOT.RooRealVar("mean_ctl", "mean_ctl", -3, -8, 8)
gx_ctl = ROOT.RooGaussian("gx_ctl", "gx_ctl", x, mean_ctl, sigma)

Construct the background pdf

In [7]:

a0_ctl = ROOT.RooRealVar("a0_ctl", "a0_ctl", -0.1, -1, 1)
a1_ctl = ROOT.RooRealVar("a1_ctl", "a1_ctl", 0.5, -0.1, 1)
px_ctl = ROOT.RooChebychev("px_ctl", "px_ctl", x, [a0_ctl, a1_ctl])

Construct the composite model

In [8]:

f_ctl = ROOT.RooRealVar("f_ctl", "f_ctl", 0.5, 0.0, 1.0)
model_ctl = ROOT.RooAddPdf("model_ctl", "model_ctl", [gx_ctl, px_ctl], [f_ctl])

Generate events for both samples¶

Generate 1000 events in x and y from model

In [9]:

data = model.generate({x}, 100)
data_ctl = model_ctl.generate({x}, 2000)

Create index category and join samples¶

Define category to distinguish physics and control samples events

In [10]:

sample = ROOT.RooCategory("sample", "sample")
sample.defineType("physics")
sample.defineType("control")

Out[10]:

False

Construct combined dataset in (x,sample)

In [11]:

combData = ROOT.RooDataSet(
    "combData",
    "combined data",
    {x},
    Index=sample,
    Import={"physics": data, "control": data_ctl},
)

Construct a simultaneous pdf in (x, sample)¶

Construct a simultaneous pdf using category sample as index: associate model with the physics state and model_ctl with the control state

In [12]:

simPdf = ROOT.RooSimultaneous("simPdf", "simultaneous pdf", {"physics": model, "control": model_ctl}, sample)

Perform a simultaneous fit¶

Perform simultaneous fit of model to data and model_ctl to data_ctl

In [13]:

fitResult = simPdf.fitTo(combData, PrintLevel=-1, Save=True)
fitResult.Print()

RooAbsTestStatistic::initSimMode: creating slave calculator #0 for state control (2000 dataset entries)
RooAbsTestStatistic::initSimMode: creating slave calculator #1 for state physics (100 dataset entries)
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 2 slave calculators.
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (gx_ctl,px_ctl)
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (gx,px)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

  RooFitResult: minimized FCN value: 5773.95, estimated distance to minimum: 9.09933e-05
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 

    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                    a0    2.4369e-01 +/-  1.75e-01
                a0_ctl   -4.3957e-03 +/-  5.29e-02
                    a1    5.2692e-02 +/-  1.78e-01
                a1_ctl    5.4474e-01 +/-  3.70e-02
                     f    6.8391e-02 +/-  3.86e-02
                 f_ctl    5.0279e-01 +/-  1.24e-02
                  mean   -4.6480e-01 +/-  2.34e-01
              mean_ctl   -3.0263e+00 +/-  1.08e-02
                 sigma    3.0786e-01 +/-  8.77e-03

Plot model slices on data slices¶

Make a frame for the physics sample

In [14]:

frame1 = x.frame(Bins=30, Title="Physics sample")

Plot all data tagged as physics sample

In [15]:

combData.plotOn(frame1, Cut="sample==sample::physics")

Out[15]:

<cppyy.gbl.RooPlot object at 0x11281630>

[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 100 events out of 2100 total events

Plot "physics" slice of simultaneous pdf. NB: You must project the sample index category with data using ProjWData as a RooSimultaneous makes no prediction on the shape in the index category and can thus not be integrated. In other words: Since the PDF doesn't know the number of events in the different category states, it doesn't know how much of each component it has to project out. This info is read from the data.

In [16]:

simPdf.plotOn(frame1, Slice=(sample, "physics"), ProjWData=(sample, combData))
simPdf.plotOn(frame1, Slice=(sample, "physics"), Components="px", ProjWData=(sample, combData), LineStyle="--")

Out[16]:

<cppyy.gbl.RooPlot object at 0x11281630>

[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) slice variable sample was not projected anyway
[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (px)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) slice variable sample was not projected anyway

The same plot for the control sample slice

In [17]:

frame2 = x.frame(Bins=30, Title="Control sample")
combData.plotOn(frame2, Cut="sample==sample::control")
simPdf.plotOn(frame2, Slice=(sample, "control"), ProjWData=(sample, combData))
simPdf.plotOn(frame2, Slice=(sample, "control"), Components="px_ctl", ProjWData=(sample, combData), LineStyle="--")

c = ROOT.TCanvas("rf501_simultaneouspdf", "rf501_simultaneouspdf", 800, 400)
c.Divide(2)
c.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
frame1.GetYaxis().SetTitleOffset(1.4)
frame1.Draw()
c.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
frame2.GetYaxis().SetTitleOffset(1.4)
frame2.Draw()

c.SaveAs("rf501_simultaneouspdf.png")

[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 2000 events out of 2100 total events
[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsReal::plotOn(model_ctl) slice variable sample was not projected anyway
[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model_ctl) directly selected PDF components: (px_ctl)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model_ctl) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsReal::plotOn(model_ctl) slice variable sample was not projected anyway

Info in <TCanvas::Print>: png file rf501_simultaneouspdf.png has been created

Draw all canvases

In [18]:

from ROOT import gROOT 
gROOT.GetListOfCanvases().Draw()