///////////////////////////////////////////////////////////////////////// // // 'Bayesian Calculator' RooStats tutorial macro #701 // author: Gregory Schott // date Sep 2009 // // This tutorial shows an example of using the BayesianCalculator class // ///////////////////////////////////////////////////////////////////////// #include "RooRealVar.h" #include "RooWorkspace.h" #include "RooDataSet.h" #include "RooPlot.h" #include "RooMsgService.h" #include "RooStats/BayesianCalculator.h" #include "RooStats/SimpleInterval.h" #include "TCanvas.h" using namespace RooFit; using namespace RooStats; void rs701_BayesianCalculator(bool useBkg = true, double confLevel = 0.90) { RooWorkspace* w = new RooWorkspace("w",true); w->factory("SUM::pdf(s[0.001,15]*Uniform(x[0,1]),b[1,0,2]*Uniform(x))"); w->factory("Gaussian::prior_b(b,1,1)"); w->factory("PROD::model(pdf,prior_b)"); RooAbsPdf* model = w->pdf("model"); // pdf*priorNuisance RooArgSet nuisanceParameters(*(w->var("b"))); RooAbsRealLValue* POI = w->var("s"); RooAbsPdf* priorPOI = (RooAbsPdf *) w->factory("Uniform::priorPOI(s)"); RooAbsPdf* priorPOI2 = (RooAbsPdf *) w->factory("GenericPdf::priorPOI2('1/sqrt(@0)',s)"); w->factory("n[3]"); // observed number of events // create a data set with n observed events RooDataSet data("data","",RooArgSet(*(w->var("x")),*(w->var("n"))),"n"); data.add(RooArgSet(*(w->var("x"))),w->var("n")->getVal()); // to suppress messgaes when pdf goes to zero RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL) ; RooArgSet * nuisPar = 0; if (useBkg) nuisPar = &nuisanceParameters; //if (!useBkg) ((RooRealVar *)w->var("b"))->setVal(0); double size = 1.-confLevel; std::cout << "\nBayesian Result using a Flat prior " << std::endl; BayesianCalculator bcalc(data,*model,RooArgSet(*POI),*priorPOI, nuisPar); bcalc.SetTestSize(size); SimpleInterval* interval = bcalc.GetInterval(); double cl = bcalc.ConfidenceLevel(); std::cout << cl <<"% CL central interval: [ " << interval->LowerLimit() << " - " << interval->UpperLimit() << " ] or " << cl+(1.-cl)/2 << "% CL limits\n"; RooPlot * plot = bcalc.GetPosteriorPlot(); TCanvas * c1 = new TCanvas("c1","Bayesian Calculator Result"); c1->Divide(1,2); c1->cd(1); plot->Draw(); c1->Update(); std::cout << "\nBayesian Result using a 1/sqrt(s) prior " << std::endl; BayesianCalculator bcalc2(data,*model,RooArgSet(*POI),*priorPOI2,nuisPar); bcalc2.SetTestSize(size); SimpleInterval* interval2 = bcalc2.GetInterval(); cl = bcalc2.ConfidenceLevel(); std::cout << cl <<"% CL central interval: [ " << interval2->LowerLimit() << " - " << interval2->UpperLimit() << " ] or " << cl+(1.-cl)/2 << "% CL limits\n"; RooPlot * plot2 = bcalc2.GetPosteriorPlot(); c1->cd(2); plot2->Draw(); gPad->SetLogy(); c1->Update(); // observe one event while expecting one background event -> the 95% CL upper limit on s is 4.10 // observe one event while expecting zero background event -> the 95% CL upper limit on s is 4.74 } rs701_BayesianCalculator.C:1 rs701_BayesianCalculator.C:2 rs701_BayesianCalculator.C:3 rs701_BayesianCalculator.C:4 rs701_BayesianCalculator.C:5 rs701_BayesianCalculator.C:6 rs701_BayesianCalculator.C:7 rs701_BayesianCalculator.C:8 rs701_BayesianCalculator.C:9 rs701_BayesianCalculator.C:10 rs701_BayesianCalculator.C:11 rs701_BayesianCalculator.C:12 rs701_BayesianCalculator.C:13 rs701_BayesianCalculator.C:14 rs701_BayesianCalculator.C:15 rs701_BayesianCalculator.C:16 rs701_BayesianCalculator.C:17 rs701_BayesianCalculator.C:18 rs701_BayesianCalculator.C:19 rs701_BayesianCalculator.C:20 rs701_BayesianCalculator.C:21 rs701_BayesianCalculator.C:22 rs701_BayesianCalculator.C:23 rs701_BayesianCalculator.C:24 rs701_BayesianCalculator.C:25 rs701_BayesianCalculator.C:26 rs701_BayesianCalculator.C:27 rs701_BayesianCalculator.C:28 rs701_BayesianCalculator.C:29 rs701_BayesianCalculator.C:30 rs701_BayesianCalculator.C:31 rs701_BayesianCalculator.C:32 rs701_BayesianCalculator.C:33 rs701_BayesianCalculator.C:34 rs701_BayesianCalculator.C:35 rs701_BayesianCalculator.C:36 rs701_BayesianCalculator.C:37 rs701_BayesianCalculator.C:38 rs701_BayesianCalculator.C:39 rs701_BayesianCalculator.C:40 rs701_BayesianCalculator.C:41 rs701_BayesianCalculator.C:42 rs701_BayesianCalculator.C:43 rs701_BayesianCalculator.C:44 rs701_BayesianCalculator.C:45 rs701_BayesianCalculator.C:46 rs701_BayesianCalculator.C:47 rs701_BayesianCalculator.C:48 rs701_BayesianCalculator.C:49 rs701_BayesianCalculator.C:50 rs701_BayesianCalculator.C:51 rs701_BayesianCalculator.C:52 rs701_BayesianCalculator.C:53 rs701_BayesianCalculator.C:54 rs701_BayesianCalculator.C:55 rs701_BayesianCalculator.C:56 rs701_BayesianCalculator.C:57 rs701_BayesianCalculator.C:58 rs701_BayesianCalculator.C:59 rs701_BayesianCalculator.C:60 rs701_BayesianCalculator.C:61 rs701_BayesianCalculator.C:62 rs701_BayesianCalculator.C:63 rs701_BayesianCalculator.C:64 rs701_BayesianCalculator.C:65 rs701_BayesianCalculator.C:66 rs701_BayesianCalculator.C:67 rs701_BayesianCalculator.C:68 rs701_BayesianCalculator.C:69 rs701_BayesianCalculator.C:70 rs701_BayesianCalculator.C:71 rs701_BayesianCalculator.C:72 rs701_BayesianCalculator.C:73 rs701_BayesianCalculator.C:74 rs701_BayesianCalculator.C:75 rs701_BayesianCalculator.C:76 rs701_BayesianCalculator.C:77 rs701_BayesianCalculator.C:78 rs701_BayesianCalculator.C:79 rs701_BayesianCalculator.C:80 rs701_BayesianCalculator.C:81 rs701_BayesianCalculator.C:82 rs701_BayesianCalculator.C:83 rs701_BayesianCalculator.C:84 rs701_BayesianCalculator.C:85 rs701_BayesianCalculator.C:86 rs701_BayesianCalculator.C:87 rs701_BayesianCalculator.C:88 |
|