Rs 7 0 1_ Bayesian Calculator

Bayesian calculator: basic exmple

Author: Gregory Schott
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Thursday, June 30, 2022 at 09:58 AM.

In [1]:
%%cpp -d
#include "RooRealVar.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooMsgService.h"

#include "RooStats/BayesianCalculator.h"
#include "RooStats/SimpleInterval.h"
#include "TCanvas.h"
In [2]:
%%cpp -d
// This is a workaround to make sure the namespace is used inside functions
using namespace RooFit;
using namespace RooStats;

Arguments are defined.

In [3]:
bool useBkg = true;
double confLevel = 0.90;
In [4]:
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
[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(w) INFO: references to all objects in this workspace will be created in CINT in 'namespace w'
input_line_53:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooRealVar& s = *(RooRealVar *)0x7fe178aa9610 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_54:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooRealVar& x = *(RooRealVar *)0x7fe178ab2540 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_55:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooUniform& pdf_1 = *(RooUniform *)0x7fe178bbf470 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_56:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooRealVar& b = *(RooRealVar *)0x7fe178bdf7a0 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_57:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooUniform& pdf_2 = *(RooUniform *)0x7fe178bbeef0 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_58:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooAddPdf& pdf = *(RooAddPdf *)0x7fe178bc46e0 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_60:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooGaussian& prior_b = *(RooGaussian *)0x7fe178aadc20 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_61:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooProdPdf& model = *(RooProdPdf *)0x7fe178d6d990 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_63:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooUniform& priorPOI = *(RooUniform *)0x7fe178d6b4b0 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_68:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooGenericPdf& priorPOI2 = *(RooGenericPdf *)0x7fe178d65b70 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^
input_line_69:1:11: error: redefinition of 'w' as different kind of symbol
namespace w { RooRealVar& n = *(RooRealVar *)0x7fe178b94de0 ; }
          ^
input_line_50:2:16: note: previous definition is here
 RooWorkspace *w = new RooWorkspace("w", true);
               ^

Create a data set with n observed events

In [5]:
RooDataSet data("data", "", RooArgSet(*(w->var("x")), *(w->var("n"))), "n");
data.add(RooArgSet(*(w->var("x"))), w->var("n")->getVal());
input_line_70:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
 RooDataSet data("data", "", RooArgSet(*(w->var("x")), *(w->var("n"))), "n");
 ^

To suppress messages when pdf goes to zero

In [6]:
RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);

RooArgSet *nuisPar = 0;
if (useBkg)
   nuisPar = &nuisanceParameters;

If (!usebkg) ((roorealvar *)w->var("b"))->setval(0);

In [7]:
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();
input_line_72:4:26: error: reference to 'data' is ambiguous
BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar);
                         ^
input_line_70:2:13: note: candidate found by name lookup is '__cling_N523::data'
 RooDataSet data("data", "", RooArgSet(*(w->var("x")), *(w->var("n"))), "n");
            ^
/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_72:4:26: error: unknown type name 'data'
BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar);
                         ^
input_line_72:4:33: error: C++ requires a type specifier for all declarations
BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar);
                                ^
input_line_72:4:58: error: C++ requires a type specifier for all declarations
BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar);
                                                         ^
input_line_72:4:68: error: unknown type name 'nuisPar'
BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar);
                                                                   ^
input_line_72:18:27: error: reference to 'data' is ambiguous
BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar);
                          ^
input_line_70:2:13: note: candidate found by name lookup is '__cling_N523::data'
 RooDataSet data("data", "", RooArgSet(*(w->var("x")), *(w->var("n"))), "n");
            ^
/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_72:18:27: error: unknown type name 'data'
BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar);
                          ^
input_line_72:18:34: error: C++ requires a type specifier for all declarations
BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar);
                                 ^
input_line_72:18:59: error: C++ requires a type specifier for all declarations
BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar);
                                                          ^
input_line_72:18:70: error: unknown type name 'nuisPar'
BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar);
                                                                     ^
input_line_72:2:2: warning: 'size' shadows a declaration with the same name in the 'std' namespace; use '::size' to reference this declaration
 double size = 1. - confLevel;
 ^

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

Draw all canvases

In [8]:
gROOT->GetListOfCanvases()->Draw()