rf501_simultaneouspdf

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

Author: Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, November 30, 2022 at 11:23 AM.

In [1]:
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
using namespace RooFit;

Create model for physics sample

Create observables

In [2]:
RooRealVar x("x", "x", -8, 8);

Construct signal pdf

In [3]:
RooRealVar mean("mean", "mean", 0, -8, 8);
RooRealVar sigma("sigma", "sigma", 0.3, 0.1, 10);
RooGaussian gx("gx", "gx", x, mean, sigma);

Construct background pdf

In [4]:
RooRealVar a0("a0", "a0", -0.1, -1, 1);
RooRealVar a1("a1", "a1", 0.004, -1, 1);
RooChebychev px("px", "px", x, RooArgSet(a0, a1));

Construct composite pdf

In [5]:
RooRealVar f("f", "f", 0.2, 0., 1.);
RooAddPdf model("model", "model", RooArgList(gx, px), f);

Create model for control sample

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

In [6]:
RooRealVar mean_ctl("mean_ctl", "mean_ctl", -3, -8, 8);
RooGaussian gx_ctl("gx_ctl", "gx_ctl", x, mean_ctl, sigma);

Construct the background pdf

In [7]:
RooRealVar a0_ctl("a0_ctl", "a0_ctl", -0.1, -1, 1);
RooRealVar a1_ctl("a1_ctl", "a1_ctl", 0.5, -0.1, 1);
RooChebychev px_ctl("px_ctl", "px_ctl", x, RooArgSet(a0_ctl, a1_ctl));

Construct the composite model

In [8]:
RooRealVar f_ctl("f_ctl", "f_ctl", 0.5, 0., 1.);
RooAddPdf model_ctl("model_ctl", "model_ctl", RooArgList(gx_ctl, px_ctl), f_ctl);

Generate events for both samples

Generate 1000 events in x and y from model

In [9]:
RooDataSet *data = model.generate(RooArgSet(x), 100);
RooDataSet *data_ctl = model_ctl.generate(RooArgSet(x), 2000);
input_line_57:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
 RooDataSet *data = model.generate(RooArgSet(x), 100);
 ^

Create index category and join samples

Define category to distinguish physics and control samples events

In [10]:
RooCategory sample("sample", "sample");
sample.defineType("physics");
sample.defineType("control");
input_line_58:2:2: warning: 'sample' shadows a declaration with the same name in the 'std' namespace; use '::sample' to reference this declaration
 RooCategory sample("sample", "sample");
 ^

Construct combined dataset in (x,sample)

In [11]:
RooDataSet combData("combData", "combined data", x, Index(sample),
                    Import({{"physics", data}, {"control", data_ctl}}));
input_line_59:2:60: error: reference to 'sample' is ambiguous
 RooDataSet combData("combData", "combined data", x, Index(sample),
                                                           ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_59:3:41: error: reference to 'data' is ambiguous
                    Import({{"physics", data}, {"control", data_ctl}}));
                                        ^
input_line_57:2:14: note: candidate found by name lookup is '__cling_N528::data'
 RooDataSet *data = model.generate(RooArgSet(x), 100);
             ^
/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^

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]:
RooSimultaneous simPdf("simPdf", "simultaneous pdf", {{"physics", &model}, {"control", &model_ctl}}, sample);
input_line_60:2:103: error: reference to 'sample' is ambiguous
 RooSimultaneous simPdf("simPdf", "simultaneous pdf", {{"physics", &model}, {"control", &model_ctl}}, sample);
                                                                                                      ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^

Perform a simultaneous fit

Perform simultaneous fit of model to data and model_ctl to data_ctl

In [13]:
std::unique_ptr<RooFitResult> fitResult{simPdf.fitTo(combData, PrintLevel(-1), Save())};
fitResult->Print();
input_line_61:2:65: error: cannot initialize an array element of type 'void *' with an rvalue of type 'RooCmdArg (*)(Int_t)' (aka 'RooCmdArg (*)(int)')
 std::unique_ptr<RooFitResult> fitResult{simPdf.fitTo(combData, PrintLevel(-1), Save())};
                                                                ^~~~~~~~~~
input_line_61:2:81: error: cannot initialize an array element of type 'void *' with an rvalue of type 'RooCmdArg (*)(bool)'
 std::unique_ptr<RooFitResult> fitResult{simPdf.fitTo(combData, PrintLevel(-1), Save())};
                                                                                ^~~~

Plot model slices on data slices

Make a frame for the physics sample

In [14]:
RooPlot *frame1 = x.frame(Bins(30), Title("Physics sample"));

Plot all data tagged as physics sample

In [15]:
combData.plotOn(frame1, Cut("sample==sample::physics"));
input_line_63:2:26: error: cannot initialize an array element of type 'void *' with an rvalue of type 'RooCmdArg (*)(const char *)'
 combData.plotOn(frame1, Cut("sample==sample::physics"));
                         ^~~

Plot "physics" slice of simultaneous pdf. NBL 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 information 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(kDashed));
input_line_64:2:30: error: reference to 'sample' is ambiguous
 simPdf.plotOn(frame1, Slice(sample, "physics"), ProjWData(sample, combData));
                             ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_64:2:60: error: reference to 'sample' is ambiguous
 simPdf.plotOn(frame1, Slice(sample, "physics"), ProjWData(sample, combData));
                                                           ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_64:2:68: error: use of undeclared identifier 'combData'
 simPdf.plotOn(frame1, Slice(sample, "physics"), ProjWData(sample, combData));
                                                                   ^
input_line_64:3:1: error: use of undeclared identifier 'simPdf'
simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData), LineStyle(kDashed));
^
input_line_64:3:29: error: reference to 'sample' is ambiguous
simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData), LineStyle(kDashed));
                            ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_64:3:77: error: reference to 'sample' is ambiguous
simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData), LineStyle(kDashed));
                                                                            ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_64:3:85: error: use of undeclared identifier 'combData'
simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData), LineStyle(kDashed));
                                                                                    ^

The same plot for the control sample slice

In [17]:
RooPlot *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(kDashed));

TCanvas *c = new TCanvas("rf501_simultaneouspdf", "rf403_simultaneouspdf", 800, 400);
c->Divide(2);
c->cd(1);
gPad->SetLeftMargin(0.15);
frame1->GetYaxis()->SetTitleOffset(1.4);
frame1->Draw();
c->cd(2);
gPad->SetLeftMargin(0.15);
frame2->GetYaxis()->SetTitleOffset(1.4);
frame2->Draw();
input_line_65:4:29: error: reference to 'sample' is ambiguous
simPdf.plotOn(frame2, Slice(sample, "control"), ProjWData(sample, combData));
                            ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_65:4:59: error: reference to 'sample' is ambiguous
simPdf.plotOn(frame2, Slice(sample, "control"), ProjWData(sample, combData));
                                                          ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_65:4:67: error: use of undeclared identifier 'combData'
simPdf.plotOn(frame2, Slice(sample, "control"), ProjWData(sample, combData));
                                                                  ^
input_line_65:5:1: error: use of undeclared identifier 'simPdf'
simPdf.plotOn(frame2, Slice(sample, "control"), Components("px_ctl"), ProjWData(sample, combData),
^
input_line_65:5:29: error: reference to 'sample' is ambiguous
simPdf.plotOn(frame2, Slice(sample, "control"), Components("px_ctl"), ProjWData(sample, combData),
                            ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_65:5:81: error: reference to 'sample' is ambiguous
simPdf.plotOn(frame2, Slice(sample, "control"), Components("px_ctl"), ProjWData(sample, combData),
                                                                                ^
input_line_58:2:14: note: candidate found by name lookup is '__cling_N529::sample'
 RooCategory sample("sample", "sample");
             ^
/usr/include/c++/9/bits/stl_algo.h:5855:5: note: candidate found by name lookup is 'std::sample'
    sample(_PopulationIterator __first, _PopulationIterator __last,
    ^
input_line_65:5:89: error: use of undeclared identifier 'combData'
simPdf.plotOn(frame2, Slice(sample, "control"), Components("px_ctl"), ProjWData(sample, combData),
                                                                                        ^
input_line_65:3:25: error: cannot initialize an array element of type 'void *' with an rvalue of type 'RooCmdArg (*)(const char *)'
combData.plotOn(frame2, Cut("sample==sample::control"));
                        ^~~

Draw all canvases

In [18]:
%jsroot on
gROOT->GetListOfCanvases()->Draw()