Logo ROOT   6.21/01
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 
22 using namespace RooFit;
23 using namespace RooStats;
24 
25 void 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[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 = 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();
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 }
return c1
Definition: legend1.C:41
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
static RooMsgService & instance()
Return reference to singleton instance.
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual Double_t LowerLimit()
virtual Double_t ConfidenceLevel() const
return confidence level
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:89
void setGlobalKillBelow(RooFit::MsgLevel level)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
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.
virtual Double_t UpperLimit()
SimpleInterval is a concrete implementation of the ConfInterval interface.
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
#define gPad
Definition: TVirtualPad.h:286
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...