HybridOriginalDemo.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roostats
/// \notebook -js
/// 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).
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \author Gregory Schott

#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"

void HybridOriginalDemo(int ntoys = 1000)
{
   using namespace RooFit;
   using namespace RooStats;

   // set RooFit random seed
   RooRandom::randomGenerator()->SetSeed(3007);

   /// build the models for background and signal+background
   RooRealVar x("x", "", -3, 3);
   RooArgList observables(x); // variables to be generated

   // gaussian signal
   RooGaussian sig_pdf("sig_pdf", "", x, RooConst(0.0), RooConst(0.8));
   RooRealVar sig_yield("sig_yield", "", 20, 0, 300);

   // flat background (extended PDF)
   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);
   sig_yield.setConstant(kTRUE);

   // total sig+bkg (extended PDF)
   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% )
   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
   RooDataSet *data = tot_pdf.generate(observables, RooFit::Extended());

   // run HybridCalculator on those inputs
   // use interface from HypoTest calculator by default

   HybridCalculatorOriginal myHybridCalc(*data, tot_pdf, bkg_ext_pdf, &nuisance_parameters, &bkg_yield_prior);

   // here I use the default test statistics: 2*lnQ (optional)
   myHybridCalc.SetTestStatistic(1);
   // myHybridCalc.SetTestStatistic(3); // profile likelihood ratio

   myHybridCalc.SetNumberOfToys(ntoys);
   myHybridCalc.UseNuisance(true);

   // for speed up generation (do binned data)
   myHybridCalc.SetGenerateBinned(false);

   // calculate by running ntoys for the S+B and B hypothesis and retrieve the result
   HybridResult *myHybridResult = myHybridCalc.GetHypoTest();

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

   /// nice plot of the results
   HybridPlot *myHybridPlot =
      myHybridResult->GetPlot("myHybridPlot", "Plot of results with HybridCalculatorOriginal", 100);
   myHybridPlot->Draw();

   /// recover and display the results
   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
   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;
}