Rs 6 0 3_ H L Factory Elaborate Example

Author: Danilo Piparo
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Friday, September 18, 2020 at 09:53 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);
RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt

[#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)

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

In [7]:
pdf->fitTo(*data);
[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
RooAbsTestStatistic::initSimMode: creating slave calculator #0 for state model1 (219 dataset entries)
RooAbsTestStatistic::initSimMode: creating slave calculator #1 for state model2 (164 dataset entries)
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 2 slave calculators.
[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (flat1)
[#1] INFO:Minization --  The following expressions will be evaluated in cache-and-track mode: (gauss1)
[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (flat2)
[#1] INFO:Minization --  The following expressions will be evaluated in cache-and-track mode: (gauss2)
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 mean1        5.00000e+01  1.00000e+01    0.00000e+00  1.00000e+02
     2 mean2        8.00000e+01  1.00000e+01    0.00000e+00  1.00000e+02
     3 nbkg1        1.00000e+02  5.00000e+01    0.00000e+00  1.00000e+03
     4 nbkg2        8.00000e+01  4.00000e+01    0.00000e+00  1.00000e+03
     5 nsig1        1.20000e+02  3.00000e+01    0.00000e+00  3.00000e+02
     6 nsig2        9.00000e+01  4.00000e+01    0.00000e+00  4.00000e+02
 **********
 **    3 **SET ERR         0.5
 **********
 **********
 **    4 **SET PRINT           1
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD        3000           1
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
 START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03
 FCN=168.565 FROM MIGRAD    STATUS=INITIATE       24 CALLS          25 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  mean1        5.00000e+01   1.00000e+01   2.01358e-01  -1.82836e+02
   2  mean2        8.00000e+01   1.00000e+01   2.57889e-01   1.04014e+01
   3  nbkg1        1.00000e+02   5.00000e+01   1.72186e-01   2.51699e+01
   4  nbkg2        8.00000e+01   4.00000e+01   1.52384e-01   2.98132e+01
   5  nsig1        1.20000e+02   3.00000e+01   2.05758e-01  -9.05084e+00
   6  nsig2        9.00000e+01   4.00000e+01   2.45245e-01  -5.18071e+00
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=165.577 FROM MIGRAD    STATUS=CONVERGED     105 CALLS         106 TOTAL
                     EDM=1.82025e-06    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  mean1        5.06657e+01   4.11801e-01   7.32781e-05  -9.51127e-02
   2  mean2        7.98756e+01   6.29188e-01   1.39674e-04   1.05318e-02
   3  nbkg1        8.74181e+01   1.04372e+01   3.24289e-04   1.44260e-02
   4  nbkg2        6.75739e+01   9.49233e+00   3.28828e-04   2.19919e-02
   5  nsig1        1.31590e+02   1.23613e+01   7.28660e-04   3.75153e-03
   6  nsig2        9.64299e+01   1.09011e+01   5.54180e-04  -3.77030e-03
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  6    ERR DEF=0.5
  1.696e-01  0.000e+00 -9.815e-02  0.000e+00  9.813e-02  0.000e+00 
  0.000e+00  3.959e-01  0.000e+00  1.140e-01  0.000e+00 -1.141e-01 
 -9.815e-02  0.000e+00  1.090e+02  0.000e+00 -2.155e+01  0.000e+00 
  0.000e+00  1.140e-01  0.000e+00  9.015e+01  0.000e+00 -2.256e+01 
  9.813e-02  0.000e+00 -2.155e+01  0.000e+00  1.532e+02  0.000e+00 
  0.000e+00 -1.141e-01  0.000e+00 -2.256e+01  0.000e+00  1.190e+02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4      5      6
        1  0.02769   1.000  0.000 -0.023  0.000  0.019  0.000
        2  0.02296   0.000  1.000  0.000  0.019  0.000 -0.017
        3  0.16798  -0.023  0.000  1.000  0.000 -0.167  0.000
        4  0.21836   0.000  0.019  0.000  1.000  0.000 -0.218
        5  0.16755   0.019  0.000 -0.167  0.000  1.000  0.000
        6  0.21817   0.000 -0.017  0.000 -0.218  0.000  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        3000
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=165.577 FROM HESSE     STATUS=OK             40 CALLS         146 TOTAL
                     EDM=1.8204e-06    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  mean1        5.06657e+01   4.11802e-01   1.46556e-05   1.33141e-02
   2  mean2        7.98756e+01   6.29186e-01   2.79348e-05   6.40394e-01
   3  nbkg1        8.74181e+01   1.04374e+01   6.48579e-05  -9.70492e-01
   4  nbkg2        6.75739e+01   9.49253e+00   6.57656e-05  -1.04486e+00
   5  nsig1        1.31590e+02   1.23615e+01   1.45732e-04  -1.23045e-01
   6  nsig2        9.64299e+01   1.09013e+01   1.10836e-04  -5.44336e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  6    ERR DEF=0.5
  1.696e-01  0.000e+00 -9.850e-02  0.000e+00  9.849e-02  0.000e+00 
  0.000e+00  3.959e-01  0.000e+00  1.133e-01  0.000e+00 -1.133e-01 
 -9.850e-02  0.000e+00  1.090e+02  0.000e+00 -2.157e+01  0.000e+00 
  0.000e+00  1.133e-01  0.000e+00  9.015e+01  0.000e+00 -2.257e+01 
  9.849e-02  0.000e+00 -2.157e+01  0.000e+00  1.532e+02  0.000e+00 
  0.000e+00 -1.133e-01  0.000e+00 -2.257e+01  0.000e+00  1.190e+02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4      5      6
        1  0.02779   1.000  0.000 -0.023  0.000  0.019  0.000
        2  0.02281   0.000  1.000  0.000  0.019  0.000 -0.017
        3  0.16807  -0.023  0.000  1.000  0.000 -0.167  0.000
        4  0.21846   0.000  0.019  0.000  1.000  0.000 -0.218
        5  0.16763   0.019  0.000 -0.167  0.000  1.000  0.000
        6  0.21827   0.000 -0.017  0.000 -0.218  0.000  1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization

--- 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();
[#1] INFO:Plotting -- RooSimultaneous::plotOn(HLFactoryElaborateExample_sigbkg) plot on x represents a slice in the index category (HLFactoryElaborateExample_category)
[#1] INFO:Plotting -- RooAbsReal::plotOn(sb_model1) slice variable HLFactoryElaborateExample_category was not projected anyway
[#1] INFO:Plotting -- RooSimultaneous::plotOn(HLFactoryElaborateExample_sigbkg) plot on x represents a slice in the index category (HLFactoryElaborateExample_category)
[#1] INFO:Plotting -- RooAbsReal::plotOn(sb_model2) slice variable HLFactoryElaborateExample_category was not projected anyway
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1

Draw all canvases

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