rs102_hypotestwithshapes.C: rs102_hypotestwithshapes for RooStats project | Roostats tutorials | rs301_splot.C: SPlot tutorial |
#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 |
|