Rs 6 0 3_ H L Factory Elaborate Example

High Level Factory: creating a complex combined model.

Author: Danilo Piparo
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Thursday, June 30, 2022 at 09:57 AM.

In [1]:
%%cpp -d
#include <fstream>
#include "TString.h"
#include "TFile.h"
#include "TROOT.h"
#include "RooGlobalFunc.h"
#include "RooWorkspace.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooStats/HLFactory.h"

Use this order for safety on library loading

In [2]:
%%cpp -d
// This is a workaround to make sure the namespace is used inside functions
using namespace RooFit;
using namespace RooStats;
In [3]:
using namespace std;

--- prepare the 2 needed datacards for this example ---

In [4]:
TString card_name("rs603_card_WsMaker.rs");
ofstream ofile(card_name);
ofile << "// The simplest card for combination\n\n";
ofile << "gauss1 = Gaussian(x[0,100],mean1[50,0,100],4);\n";
ofile << "flat1 = Polynomial(x,0);\n";
ofile << "sb_model1 = SUM(nsig1[120,0,300]*gauss1 , nbkg1[100,0,1000]*flat1);\n\n";
ofile << "echo In the middle!;\n\n";
ofile << "gauss2 = Gaussian(x,mean2[80,0,100],5);\n";
ofile << "flat2 = Polynomial(x,0);\n";
ofile << "sb_model2 = SUM(nsig2[90,0,400]*gauss2 , nbkg2[80,0,1000]*flat2);\n\n";
ofile << "echo At the end!;\n";
ofile.close();

TString card_name2("rs603_card.rs");
ofstream ofile2(card_name2);
ofile2 << "// The simplest card for combination\n\n";
ofile2 << "gauss1 = Gaussian(x[0,100],mean1[50,0,100],4);\n";
ofile2 << "flat1 = Polynomial(x,0);\n";
ofile2 << "sb_model1 = SUM(nsig1[120,0,300]*gauss1 , nbkg1[100,0,1000]*flat1);\n\n";
ofile2 << "echo In the middle!;\n\n";
ofile2 << "gauss2 = Gaussian(x,mean2[80,0,100],5);\n";
ofile2 << "flat2 = Polynomial(x,0);\n";
ofile2 << "sb_model2 = SUM(nsig2[90,0,400]*gauss2 , nbkg2[80,0,1000]*flat2);\n\n";
ofile2 << "#include rs603_included_card.rs;\n\n";
ofile2 << "echo At the end!;\n";
ofile2.close();

TString card_name3("rs603_included_card.rs");
ofstream ofile3(card_name3);
ofile3 << "echo Now reading the included file!;\n\n";
ofile3 << "echo Including datasets in a Workspace in a Root file...;\n";
ofile3 << "data1 = import(rs603_infile.root,\n";
ofile3 << "               rs603_ws,\n";
ofile3 << "               data1);\n\n";
ofile3 << "data2 = import(rs603_infile.root,\n";
ofile3 << "               rs603_ws,\n";
ofile3 << "               data2);\n";
ofile3.close();

--- produce the two separate datasets into a workspace ---

In [5]:
HLFactory hlf("HLFactoryComplexExample", "rs603_card_WsMaker.rs", false);

auto x = static_cast<RooRealVar *>(hlf.GetWs()->arg("x"));
auto pdf1 = hlf.GetWs()->pdf("sb_model1");
auto pdf2 = hlf.GetWs()->pdf("sb_model2");

RooWorkspace w("rs603_ws");

auto data1 = pdf1->generate(RooArgSet(*x), Extended());
data1->SetName("data1");
w.import(*data1);

auto data2 = pdf2->generate(RooArgSet(*x), Extended());
data2->SetName("data2");
w.import(*data2);
[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(HLFactoryComplexExample_ws) INFO: references to all objects in this workspace will be created in CINT in 'namespace HLFactoryComplexExample_ws'
[HLFactoryComplexExample] echo: In the middle!
[HLFactoryComplexExample] echo: At the end!
[#1] INFO:ObjectHandling -- RooWorkspace::import(rs603_ws) importing dataset data1
[#1] INFO:ObjectHandling -- RooWorkspace::import(rs603_ws) importing RooRealVar::x
[#1] INFO:ObjectHandling -- RooWorkspace::import(rs603_ws) importing dataset data2

--- write the workspace into a rootfile ---

In [6]:
TFile outfile("rs603_infile.root", "RECREATE");
w.Write();
outfile.Close();

cout << "-------------------------------------------------------------------\n"
     << " Rootfile and Workspace prepared \n"
     << "-------------------------------------------------------------------\n";

HLFactory hlf_2("HLFactoryElaborateExample", "rs603_card.rs", false);

x = hlf_2.GetWs()->var("x");
pdf1 = hlf_2.GetWs()->pdf("sb_model1");
pdf2 = hlf_2.GetWs()->pdf("sb_model2");

hlf_2.AddChannel("model1", "sb_model1", "flat1", "data1");
hlf_2.AddChannel("model2", "sb_model2", "flat2", "data2");

auto data = hlf_2.GetTotDataSet();
auto pdf = hlf_2.GetTotSigBkgPdf();
auto thecat = hlf_2.GetTotCategory();
-------------------------------------------------------------------
 Rootfile and Workspace prepared 
-------------------------------------------------------------------
[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(HLFactoryElaborateExample_ws) INFO: references to all objects in this workspace will be created in CINT in 'namespace HLFactoryElaborateExample_ws'
[HLFactoryElaborateExample] echo: In the middle!
[HLFactoryElaborateExample] echo: Now reading the included file!
[HLFactoryElaborateExample] echo: Including datasets in a Workspace in a Root file...
[#1] INFO:ObjectHandling -- RooWorkspace::import(HLFactoryElaborateExample_ws) importing dataset data1
[#1] INFO:ObjectHandling -- RooWorkspace::import(HLFactoryElaborateExample_ws) importing dataset data2
[HLFactoryElaborateExample] echo: At the end!
RooDataSet::data1[x] = 219 entries
RooCategory::HLFactoryElaborateExample_category = model2(idx = 1)

input_line_69:19:1: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
auto data = hlf_2.GetTotDataSet();
^

--- perform extended ml fit of composite pdf to toy data ---

In [7]:
pdf->fitTo(*data);
input_line_119:2:14: error: reference to 'data' is ambiguous
 pdf->fitTo(*data);
             ^
input_line_69:19:6: note: candidate found by name lookup is '__cling_N523::data'
auto data = hlf_2.GetTotDataSet();
     ^
/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
    ^

--- plot toy data and composite pdf overlaid ---

In [8]:
auto xframe = x->frame();

data->plotOn(xframe);
thecat->setIndex(0);
pdf->plotOn(xframe, Slice(*thecat), ProjWData(*thecat, *data));

thecat->setIndex(1);
pdf->plotOn(xframe, Slice(*thecat), ProjWData(*thecat, *data));

gROOT->SetStyle("Plain");

xframe->Draw();
input_line_120:4:1: error: reference to 'data' is ambiguous
data->plotOn(xframe);
^
input_line_69:19:6: note: candidate found by name lookup is '__cling_N523::data'
auto data = hlf_2.GetTotDataSet();
     ^
/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
    ^
input_line_120:6:57: error: reference to 'data' is ambiguous
pdf->plotOn(xframe, Slice(*thecat), ProjWData(*thecat, *data));
                                                        ^
input_line_69:19:6: note: candidate found by name lookup is '__cling_N523::data'
auto data = hlf_2.GetTotDataSet();
     ^
/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
    ^
input_line_120:9:57: error: reference to 'data' is ambiguous
pdf->plotOn(xframe, Slice(*thecat), ProjWData(*thecat, *data));
                                                        ^
input_line_69:19:6: note: candidate found by name lookup is '__cling_N523::data'
auto data = hlf_2.GetTotDataSet();
     ^
/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
    ^

Draw all canvases

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