Logo ROOT   6.12/07
Reference Guide
rs701_BayesianCalculator.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook
4 /// 'Bayesian Calculator' RooStats tutorial macro #701
5 ///
6 /// This tutorial shows an example of using the BayesianCalculator class
7 ///
8 /// \macro_image
9 /// \macro_output
10 /// \macro_code
11 ///
12 /// \author Gregory Schott
13 
14 #include "RooRealVar.h"
15 #include "RooWorkspace.h"
16 #include "RooDataSet.h"
17 #include "RooPlot.h"
18 #include "RooMsgService.h"
19 
22 #include "TCanvas.h"
23 
24 using namespace RooFit;
25 using namespace RooStats;
26 
27 void rs701_BayesianCalculator(bool useBkg = true, double confLevel = 0.90)
28 {
29 
30 
31  RooWorkspace* w = new RooWorkspace("w",true);
32  w->factory("SUM::pdf(s[0.001,15]*Uniform(x[0,1]),b[1,0,2]*Uniform(x))");
33  w->factory("Gaussian::prior_b(b,1,1)");
34  w->factory("PROD::model(pdf,prior_b)");
35  RooAbsPdf* model = w->pdf("model"); // pdf*priorNuisance
36  RooArgSet nuisanceParameters(*(w->var("b")));
37 
38 
39 
40  RooAbsRealLValue* POI = w->var("s");
41  RooAbsPdf* priorPOI = (RooAbsPdf *) w->factory("Uniform::priorPOI(s)");
42  RooAbsPdf* priorPOI2 = (RooAbsPdf *) w->factory("GenericPdf::priorPOI2('1/sqrt(@0)',s)");
43 
44  w->factory("n[3]"); // observed number of events
45  // create a data set with n observed events
46  RooDataSet data("data","",RooArgSet(*(w->var("x")),*(w->var("n"))),"n");
47  data.add(RooArgSet(*(w->var("x"))),w->var("n")->getVal());
48 
49  // to suppress messages when pdf goes to zero
51 
52  RooArgSet * nuisPar = 0;
53  if (useBkg) nuisPar = &nuisanceParameters;
54  //if (!useBkg) ((RooRealVar *)w->var("b"))->setVal(0);
55 
56  double size = 1.-confLevel;
57  std::cout << "\nBayesian Result using a Flat prior " << std::endl;
58  BayesianCalculator bcalc(data,*model,RooArgSet(*POI),*priorPOI, nuisPar);
59  bcalc.SetTestSize(size);
60  SimpleInterval* interval = bcalc.GetInterval();
61  double cl = bcalc.ConfidenceLevel();
62  std::cout << cl <<"% CL central interval: [ " << interval->LowerLimit() << " - " << interval->UpperLimit()
63  << " ] or "
64  << cl+(1.-cl)/2 << "% CL limits\n";
65  RooPlot * plot = bcalc.GetPosteriorPlot();
66  TCanvas * c1 = new TCanvas("c1","Bayesian Calculator Result");
67  c1->Divide(1,2);
68  c1->cd(1);
69  plot->Draw();
70  c1->Update();
71 
72  std::cout << "\nBayesian Result using a 1/sqrt(s) prior " << std::endl;
73  BayesianCalculator bcalc2(data,*model,RooArgSet(*POI),*priorPOI2,nuisPar);
74  bcalc2.SetTestSize(size);
75  SimpleInterval* interval2 = bcalc2.GetInterval();
76  cl = bcalc2.ConfidenceLevel();
77  std::cout << cl <<"% CL central interval: [ " << interval2->LowerLimit() << " - " << interval2->UpperLimit()
78  << " ] or "
79  << cl+(1.-cl)/2 << "% CL limits\n";
80 
81  RooPlot * plot2 = bcalc2.GetPosteriorPlot();
82  c1->cd(2);
83  plot2->Draw();
84  gPad->SetLogy();
85  c1->Update();
86 
87  // observe one event while expecting one background event -> the 95% CL upper limit on s is 4.10
88  // observe one event while expecting zero background event -> the 95% CL upper limit on s is 4.74
89 }
virtual void Draw(Option_t *option="")=0
Default Draw method for all objects.
return c1
Definition: legend1.C:41
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
static RooMsgService & instance()
Return reference to singleton instance.
virtual Double_t LowerLimit()
virtual Double_t ConfidenceLevel() const
return confidence level
void setGlobalKillBelow(RooFit::MsgLevel level)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:41
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 is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1153
#define gPad
Definition: TVirtualPad.h:285
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2248
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:42