#include "RooRandom.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooArgSet.h"
#include "RooAddPdf.h"
#include "RooDataSet.h"
#include "RooExtendPdf.h"
#ifndef __CINT__ // problem including this file with CINT
#include "RooGlobalFunc.h"
#endif
#include "RooStats/HybridCalculator.h"
#include "RooStats/HybridResult.h"
#include "RooStats/HybridPlot.h"
void rs201_hybridcalculator(int ntoys = 3000)
{
using namespace RooFit;
using namespace RooStats;
RooRandom::randomGenerator()->SetSeed(3007);
RooRealVar x("x","",-3,3);
RooArgList observables(x);
RooRealVar sig_mean("sig_mean","",0);
RooRealVar sig_sigma("sig_sigma","",0.8);
RooGaussian sig_pdf("sig_pdf","",x,sig_mean,sig_sigma);
RooRealVar sig_yield("sig_yield","",20,0,300);
RooRealVar bkg_slope("bkg_slope","",0);
RooPolynomial bkg_pdf("bkg_pdf","",x,bkg_slope);
RooRealVar bkg_yield("bkg_yield","",40,0,300);
RooExtendPdf bkg_ext_pdf("bkg_ext_pdf","",bkg_pdf,bkg_yield);
RooAddPdf tot_pdf("tot_pdf","",RooArgList(sig_pdf,bkg_pdf),RooArgList(sig_yield,bkg_yield));
RooGaussian bkg_yield_prior("bkg_yield_prior","",bkg_yield,RooConst(bkg_yield.getVal()),RooConst(10.));
RooArgSet nuisance_parameters(bkg_yield);
RooDataSet* data = tot_pdf.generate(observables,RooFit::Extended());
#ifndef USE_OLD_API
HybridCalculator myHybridCalc("myHybridCalc","HybridCalculator example",
*data, tot_pdf , bkg_ext_pdf ,
&nuisance_parameters, &bkg_yield_prior);
myHybridCalc.SetNumberOfToys(ntoys);
HybridResult* myHybridResult = myHybridCalc.GetHypoTest();
#else // use old api
HybridCalculator myHybridCalc("myHybridCalc","HybridCalculator example",tot_pdf,bkg_ext_pdf,observables,nuisance_parameters,bkg_yield_prior);
myHybridCalc.SetTestStatistics(1);
HybridResult* myHybridResult = myHybridCalc.Calculate(*data,ntoys,true);
#endif
if (! myHybridResult) {
std::cerr << "\nError returned from Hypothesis test" << std::endl;
return;
}
HybridPlot* myHybridPlot = myHybridResult->GetPlot("myHybridPlot","Plot of results with HybridCalculator",100);
myHybridPlot->Draw();
double clsb_data = myHybridResult->CLsplusb();
double clb_data = myHybridResult->CLb();
double cls_data = myHybridResult->CLs();
double data_significance = myHybridResult->Significance();
double mean_sb_toys_test_stat = myHybridPlot->GetSBmean();
myHybridResult->SetDataTestStatistics(mean_sb_toys_test_stat);
double toys_significance = myHybridResult->Significance();
std::cout << "Completed HybridCalculator example:\n";
std::cout << " - -2lnQ = " << myHybridResult->GetTestStat_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;
}
|
|