 ROOT   Reference Guide rs701_BayesianCalculator.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook
4/// Bayesian calculator: basic exmple
5///
6/// \macro_image
7/// \macro_output
8/// \macro_code
9///
10/// \author Gregory Schott
11
12#include "RooRealVar.h"
13#include "RooWorkspace.h"
14#include "RooDataSet.h"
15#include "RooPlot.h"
16#include "RooMsgService.h"
17
20#include "TCanvas.h"
21
22using namespace RooFit;
23using namespace RooStats;
24
25void rs701_BayesianCalculator(bool useBkg = true, double confLevel = 0.90)
26{
27
28 RooWorkspace *w = new RooWorkspace("w", true);
29 w->factory("SUM::pdf(s[0.001,15]*Uniform(x[0,1]),b[1,0,2]*Uniform(x))");
30 w->factory("Gaussian::prior_b(b,1,1)");
31 w->factory("PROD::model(pdf,prior_b)");
32 RooAbsPdf *model = w->pdf("model"); // pdf*priorNuisance
33 RooArgSet nuisanceParameters(*(w->var("b")));
34
35 RooAbsRealLValue *POI = w->var("s");
36 RooAbsPdf *priorPOI = (RooAbsPdf *)w->factory("Uniform::priorPOI(s)");
37 RooAbsPdf *priorPOI2 = (RooAbsPdf *)w->factory("GenericPdf::priorPOI2('1/sqrt(@0)',s)");
38
39 w->factory("n"); // observed number of events
40 // create a data set with n observed events
41 RooDataSet data("data", "", RooArgSet(*(w->var("x")), *(w->var("n"))), "n");
43
44 // to suppress messages when pdf goes to zero
46
47 RooArgSet *nuisPar = 0;
48 if (useBkg)
49 nuisPar = &nuisanceParameters;
50 // if (!useBkg) ((RooRealVar *)w->var("b"))->setVal(0);
51
52 double size = 1. - confLevel;
53 std::cout << "\nBayesian Result using a Flat prior " << std::endl;
54 BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar);
55 bcalc.SetTestSize(size);
56 SimpleInterval *interval = bcalc.GetInterval();
57 double cl = bcalc.ConfidenceLevel();
58 std::cout << cl << "% CL central interval: [ " << interval->LowerLimit() << " - " << interval->UpperLimit()
59 << " ] or " << cl + (1. - cl) / 2 << "% CL limits\n";
60 RooPlot *plot = bcalc.GetPosteriorPlot();
61 TCanvas *c1 = new TCanvas("c1", "Bayesian Calculator Result");
62 c1->Divide(1, 2);
63 c1->cd(1);
64 plot->Draw();
65 c1->Update();
66
67 std::cout << "\nBayesian Result using a 1/sqrt(s) prior " << std::endl;
68 BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar);
69 bcalc2.SetTestSize(size);
70 SimpleInterval *interval2 = bcalc2.GetInterval();
71 cl = bcalc2.ConfidenceLevel();
72 std::cout << cl << "% CL central interval: [ " << interval2->LowerLimit() << " - " << interval2->UpperLimit()
73 << " ] or " << cl + (1. - cl) / 2 << "% CL limits\n";
74
75 RooPlot *plot2 = bcalc2.GetPosteriorPlot();
76 c1->cd(2);
77 plot2->Draw();
79 c1->Update();
80
81 // observe one event while expecting one background event -> the 95% CL upper limit on s is 4.10
82 // observe one event while expecting zero background event -> the 95% CL upper limit on s is 4.74
83}
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
SimpleInterval is a concrete implementation of the ConfInterval interface.
virtual Double_t ConfidenceLevel() const
return confidence level
virtual Double_t UpperLimit()
virtual Double_t LowerLimit()
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
RooFactoryWSTool & factory()
Return instance to factory tool.
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
The Canvas class.
Definition: TCanvas.h:27
return c1
Definition: legend1.C:41
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition: Asimov.h:19