Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rs701_BayesianCalculator.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook
4/// Bayesian calculator: basic example
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");
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[3]"); // 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");
42 data.add(RooArgSet(*(w->var("x"))), w->var("n")->getVal());
43
44 // to suppress messages when pdf goes to zero
46
47 RooArgSet *nuisPar = nullptr;
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();
78 gPad->SetLogy();
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}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
#define gPad
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:57
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:43
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:652
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
SimpleInterval is a concrete implementation of the ConfInterval interface.
virtual double UpperLimit()
return the interval upper limit
double ConfidenceLevel() const override
return the confidence interval
virtual double LowerLimit()
return the interval lower limit
Persistable container for RooFit projects.
The Canvas class.
Definition TCanvas.h:23
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
return c1
Definition legend1.C:41
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
Namespace for the RooStats classes.
Definition Asimov.h:19