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

#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)
{
  //***********************************************************************//
  // This macro show an example on how to use RooStats/HybridCalculator    //
  //***********************************************************************//
  //
  // With this example, you should get: CL_sb = 0.130 and CL_b = 0.946
  // (if data had -2lnQ = -3.0742). You can compare to the expected plot:
  // http://www-ekp.physik.uni-karlsruhe.de/~schott/roostats/hybridplot_example.png

  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
  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);

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

  // 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
#ifndef USE_OLD_API

    HybridCalculator myHybridCalc("myHybridCalc","HybridCalculator example",
                                *data, tot_pdf , bkg_ext_pdf ,
                                &nuisance_parameters, &bkg_yield_prior); 

  myHybridCalc.SetNumberOfToys(ntoys); 
  //myHybridCalc.UseNuisance(false);                            

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

#else // use old api 
  HybridCalculator myHybridCalc("myHybridCalc","HybridCalculator example",tot_pdf,bkg_ext_pdf,observables,nuisance_parameters,bkg_yield_prior);

  // here I use the default test statistics: 2*lnQ (optional)
  myHybridCalc.SetTestStatistics(1);

  /// run ntoys toys with gaussian prior on the background yield
  HybridResult* myHybridResult = myHybridCalc.Calculate(*data,ntoys,true);
#endif

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

  /// run 1000 toys without gaussian prior on the background yield
  //HybridResult* myHybridResult = myHybridCalc.Calculate(*data,1000,false);

  /// save the toy-MC results to file, this way splitting into sub-batch jobs is possible
  //TFile fileOut("some_hybridresult.root","RECREATE");
  //fileOut.cd();
  //myHybridResult.Write();
  //fileOut.Close();

  /// read the results from a file
  //TFile fileIn("some_hybridresult.root");
  //HybridResult* myOtherHybridResult = (HybridResult*) fileIn.Get("myHybridCalc");

  /// example on how to merge with toy-MC results obtained in another job
  //HybridResult* mergedHybridResult = new HybridResult("mergedHybridResult","this object holds merged results");
  //mergedHybridResult->Add(myHybridResult);
  //mergedHybridResult->Add(myOtherHybridResult);
  /// or
  //myHybridResult->Add(myOtherHybridResult);

  /// nice plot of the results
  HybridPlot* myHybridPlot = myHybridResult->GetPlot("myHybridPlot","Plot of results with HybridCalculator",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();

  /// 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 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;
}


 rs201_hybridcalculator.C:1
 rs201_hybridcalculator.C:2
 rs201_hybridcalculator.C:3
 rs201_hybridcalculator.C:4
 rs201_hybridcalculator.C:5
 rs201_hybridcalculator.C:6
 rs201_hybridcalculator.C:7
 rs201_hybridcalculator.C:8
 rs201_hybridcalculator.C:9
 rs201_hybridcalculator.C:10
 rs201_hybridcalculator.C:11
 rs201_hybridcalculator.C:12
 rs201_hybridcalculator.C:13
 rs201_hybridcalculator.C:14
 rs201_hybridcalculator.C:15
 rs201_hybridcalculator.C:16
 rs201_hybridcalculator.C:17
 rs201_hybridcalculator.C:18
 rs201_hybridcalculator.C:19
 rs201_hybridcalculator.C:20
 rs201_hybridcalculator.C:21
 rs201_hybridcalculator.C:22
 rs201_hybridcalculator.C:23
 rs201_hybridcalculator.C:24
 rs201_hybridcalculator.C:25
 rs201_hybridcalculator.C:26
 rs201_hybridcalculator.C:27
 rs201_hybridcalculator.C:28
 rs201_hybridcalculator.C:29
 rs201_hybridcalculator.C:30
 rs201_hybridcalculator.C:31
 rs201_hybridcalculator.C:32
 rs201_hybridcalculator.C:33
 rs201_hybridcalculator.C:34
 rs201_hybridcalculator.C:35
 rs201_hybridcalculator.C:36
 rs201_hybridcalculator.C:37
 rs201_hybridcalculator.C:38
 rs201_hybridcalculator.C:39
 rs201_hybridcalculator.C:40
 rs201_hybridcalculator.C:41
 rs201_hybridcalculator.C:42
 rs201_hybridcalculator.C:43
 rs201_hybridcalculator.C:44
 rs201_hybridcalculator.C:45
 rs201_hybridcalculator.C:46
 rs201_hybridcalculator.C:47
 rs201_hybridcalculator.C:48
 rs201_hybridcalculator.C:49
 rs201_hybridcalculator.C:50
 rs201_hybridcalculator.C:51
 rs201_hybridcalculator.C:52
 rs201_hybridcalculator.C:53
 rs201_hybridcalculator.C:54
 rs201_hybridcalculator.C:55
 rs201_hybridcalculator.C:56
 rs201_hybridcalculator.C:57
 rs201_hybridcalculator.C:58
 rs201_hybridcalculator.C:59
 rs201_hybridcalculator.C:60
 rs201_hybridcalculator.C:61
 rs201_hybridcalculator.C:62
 rs201_hybridcalculator.C:63
 rs201_hybridcalculator.C:64
 rs201_hybridcalculator.C:65
 rs201_hybridcalculator.C:66
 rs201_hybridcalculator.C:67
 rs201_hybridcalculator.C:68
 rs201_hybridcalculator.C:69
 rs201_hybridcalculator.C:70
 rs201_hybridcalculator.C:71
 rs201_hybridcalculator.C:72
 rs201_hybridcalculator.C:73
 rs201_hybridcalculator.C:74
 rs201_hybridcalculator.C:75
 rs201_hybridcalculator.C:76
 rs201_hybridcalculator.C:77
 rs201_hybridcalculator.C:78
 rs201_hybridcalculator.C:79
 rs201_hybridcalculator.C:80
 rs201_hybridcalculator.C:81
 rs201_hybridcalculator.C:82
 rs201_hybridcalculator.C:83
 rs201_hybridcalculator.C:84
 rs201_hybridcalculator.C:85
 rs201_hybridcalculator.C:86
 rs201_hybridcalculator.C:87
 rs201_hybridcalculator.C:88
 rs201_hybridcalculator.C:89
 rs201_hybridcalculator.C:90
 rs201_hybridcalculator.C:91
 rs201_hybridcalculator.C:92
 rs201_hybridcalculator.C:93
 rs201_hybridcalculator.C:94
 rs201_hybridcalculator.C:95
 rs201_hybridcalculator.C:96
 rs201_hybridcalculator.C:97
 rs201_hybridcalculator.C:98
 rs201_hybridcalculator.C:99
 rs201_hybridcalculator.C:100
 rs201_hybridcalculator.C:101
 rs201_hybridcalculator.C:102
 rs201_hybridcalculator.C:103
 rs201_hybridcalculator.C:104
 rs201_hybridcalculator.C:105
 rs201_hybridcalculator.C:106
 rs201_hybridcalculator.C:107
 rs201_hybridcalculator.C:108
 rs201_hybridcalculator.C:109
 rs201_hybridcalculator.C:110
 rs201_hybridcalculator.C:111
 rs201_hybridcalculator.C:112
 rs201_hybridcalculator.C:113
 rs201_hybridcalculator.C:114
 rs201_hybridcalculator.C:115
 rs201_hybridcalculator.C:116
 rs201_hybridcalculator.C:117
 rs201_hybridcalculator.C:118
 rs201_hybridcalculator.C:119
 rs201_hybridcalculator.C:120
 rs201_hybridcalculator.C:121
 rs201_hybridcalculator.C:122
 rs201_hybridcalculator.C:123
 rs201_hybridcalculator.C:124
 rs201_hybridcalculator.C:125
 rs201_hybridcalculator.C:126
 rs201_hybridcalculator.C:127
 rs201_hybridcalculator.C:128
 rs201_hybridcalculator.C:129
 rs201_hybridcalculator.C:130
 rs201_hybridcalculator.C:131
 rs201_hybridcalculator.C:132
 rs201_hybridcalculator.C:133
 rs201_hybridcalculator.C:134
 rs201_hybridcalculator.C:135
 rs201_hybridcalculator.C:136
 rs201_hybridcalculator.C:137
 rs201_hybridcalculator.C:138
 rs201_hybridcalculator.C:139
 rs201_hybridcalculator.C:140
thumb