# 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 Thursday, June 30, 2022 at 09:45 AM.

In [ ]:
%%cpp -d
#include "RooRandom.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooArgSet.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 [ ]:
int ntoys = 1000;

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


/ build the models for background and signal+background

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


Gaussian signal

In [ ]:
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 [ ]:
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 [ ]:
sig_yield.setConstant(kTRUE);


Total sig+bkg (extended pdf)

In [ ]:
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 [ ]:
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 [ ]:
RooDataSet *data = tot_pdf.generate(observables, RooFit::Extended());


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

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


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

In [ ]:
myHybridCalc.SetTestStatistic(1);


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

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


For speed up generation (do binned data)

In [ ]:
myHybridCalc.SetGenerateBinned(false);


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

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

if (!myHybridResult) {
std::cerr << "\nError returned from Hypothesis test" << std::endl;
return;
}


/ nice plot of the results

In [ ]:
HybridPlot *myHybridPlot =
myHybridResult->GetPlot("myHybridPlot", "Plot of results with HybridCalculatorOriginal", 100);
myHybridPlot->Draw();


/ recover and display the results

In [ ]:
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 [ ]:
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;


Draw all canvases

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