Hybrid Original Demo

Example on how to use the HybridCalculatorOriginal class

With this example, you should get: CL_sb = 0.130 and CL_b = 0.946 (if data had -2lnQ = -3.0742).

Author: Gregory Schott
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Monday, February 17, 2020 at 03:13 AM.

In [1]:
%%cpp -d
#include "RooRandom.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooArgSet.h"
#include "RooAddPdf.h"
#include "RooDataSet.h"
#include "RooExtendPdf.h"
#include "RooConstVar.h"
#include "RooStats/HybridCalculatorOriginal.h"
#include "RooStats/HybridResult.h"
#include "RooStats/HybridPlot.h"

Arguments are defined.

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

Set roofit random seed

In [4]:
RooRandom::randomGenerator()->SetSeed(3007);
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

/ build the models for background and signal+background

In [5]:
RooRealVar x("x", "", -3, 3);
RooArgList observables(x); // variables to be generated

Gaussian signal

In [6]:
RooGaussian sig_pdf("sig_pdf", "", x, RooConst(0.0), RooConst(0.8));
RooRealVar sig_yield("sig_yield", "", 20, 0, 300);

Flat background (extended pdf)

In [7]:
RooPolynomial bkg_pdf("bkg_pdf", "", x, RooConst(0));
RooRealVar bkg_yield("bkg_yield", "", 40, 0, 300);
RooExtendPdf bkg_ext_pdf("bkg_ext_pdf", "", bkg_pdf, bkg_yield);

Bkg_yield.setconstant(ktrue);

In [8]:
sig_yield.setConstant(kTRUE);

Total sig+bkg (extended pdf)

In [9]:
RooAddPdf tot_pdf("tot_pdf", "", RooArgList(sig_pdf, bkg_pdf), RooArgList(sig_yield, bkg_yield));

Build the prior pdf on the parameters to be integrated gaussian contraint on the background yield ( N_B = 40 +/- 10 ie. 25% )

In [10]:
RooGaussian bkg_yield_prior("bkg_yield_prior", "", bkg_yield, RooConst(bkg_yield.getVal()), RooConst(10.));

RooArgSet nuisance_parameters(bkg_yield); // variables to be integrated

/ generate a data sample

In [11]:
RooDataSet *data = tot_pdf.generate(observables, RooFit::Extended());

Run hybridcalculator on those inputs use interface from HypoTest calculator by default

In [12]:
HybridCalculatorOriginal myHybridCalc(*data, tot_pdf, bkg_ext_pdf, &nuisance_parameters, &bkg_yield_prior);

Here i use the default test statistics: 2*lnq (optional)

In [13]:
myHybridCalc.SetTestStatistic(1);

Myhybridcalc.setteststatistic(3); // profile likelihood ratio

In [14]:
myHybridCalc.SetNumberOfToys(ntoys);
myHybridCalc.UseNuisance(true);

For speed up generation (do binned data)

In [15]:
myHybridCalc.SetGenerateBinned(false);

Calculate by running ntoys for the s+b and b hypothesis and retrieve the result

In [16]:
HybridResult *myHybridResult = myHybridCalc.GetHypoTest();

if (!myHybridResult) {
   std::cerr << "\nError returned from Hypothesis test" << std::endl;
   return;
}
Test statistics has been evaluated for data
HybridCalculatorOriginal: run 1000 toy-MC experiments
with test statistics index: 1
marginalize nuisance parameters 
....... toy number 0 / 1000
....... toy number 500 / 1000

/ nice plot of the results

In [17]:
HybridPlot *myHybridPlot =
   myHybridResult->GetPlot("myHybridPlot", "Plot of results with HybridCalculatorOriginal", 100);
myHybridPlot->Draw();
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1

/ recover and display the results

In [18]:
double clsb_data = myHybridResult->CLsplusb();
double clb_data = myHybridResult->CLb();
double cls_data = myHybridResult->CLs();
double data_significance = myHybridResult->Significance();
double min2lnQ_data = myHybridResult->GetTestStat_data();

/ compute the mean expected significance from toys

In [19]:
double mean_sb_toys_test_stat = myHybridPlot->GetSBmean();
myHybridResult->SetDataTestStatistics(mean_sb_toys_test_stat);
double toys_significance = myHybridResult->Significance();

std::cout << "Completed HybridCalculatorOriginal example:\n";
std::cout << " - -2lnQ = " << min2lnQ_data << endl;
std::cout << " - CL_sb = " << clsb_data << std::endl;
std::cout << " - CL_b  = " << clb_data << std::endl;
std::cout << " - CL_s  = " << cls_data << std::endl;
std::cout << " - significance of data  = " << data_significance << std::endl;
std::cout << " - mean significance of toys  = " << toys_significance << std::endl;
Completed HybridCalculatorOriginal example:
 - -2lnQ = -5.97789
 - CL_sb = 0.194
 - CL_b  = 0.04
 - CL_s  = 4.85
 - significance of data  = 1.75069
 - mean significance of toys  = 2.65207

Draw all canvases

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