Logo ROOT   6.07/09
Reference Guide
rs701_BayesianCalculator.C File Reference

Detailed Description

View in nbviewer Open in SWAN 'Bayesian Calculator' RooStats tutorial macro #701

This tutorial shows an example of using the BayesianCalculator class

pict1_rs701_BayesianCalculator.C.png
Processing /mnt/vdb/lsf/workspace/root-makedoc-v608/rootspi/rdoc/src/v6-08-00-patches/tutorials/roostats/rs701_BayesianCalculator.C...
RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(w) INFO: references to all objects in this workspace will be created in CINT in 'namespace w'
Bayesian Result using a Flat prior
0.9% CL central interval: [ 0.50606 - 6.89326 ] or 0.95% CL limits
Bayesian Result using a 1/sqrt(s) prior
0.9% CL central interval: [ 0.0746294 - 5.85425 ] or 0.95% CL limits
#include "RooRealVar.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooMsgService.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 messages when pdf goes to zero
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
}
Author
Gregory Schott

Definition in file rs701_BayesianCalculator.C.