Bayesian calculator: basic example
Author: Gregory Schott
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:18 PM.
%%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"
using namespace RooFit;
using namespace RooStats;
Arguments are defined.
bool useBkg = true;
double confLevel = 0.90;
RooWorkspace *w = new RooWorkspace("w");
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
create a data set with n observed events
RooDataSet data("data", "", {*w->var("x"), *w->var("n")}, RooFit::WeightVar("n"));
data.add({*(w->var("x"))}, w->var("n")->getVal());
input_line_58:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration RooDataSet data("data", "", {*w->var("x"), *w->var("n")}, RooFit::WeightVar("n")); ^
to suppress messages when pdf goes to zero
RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
RooArgSet *nuisPar = nullptr;
if (useBkg)
nuisPar = &nuisanceParameters;
if (!useBkg) ((RooRealVar *)w->var("b"))->setVal(0);
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_60:4:26: error: reference to 'data' is ambiguous BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar); ^ input_line_58:2:13: note: candidate found by name lookup is 'data' RooDataSet data("data", "", {*w->var("x"), *w->var("n")}, RooFit::WeightVar("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_60:4:26: error: unknown type name 'data' BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar); ^ input_line_60:4:33: error: a type specifier is required for all declarations BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar); ^ input_line_60:4:58: error: a type specifier is required for all declarations BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar); ^ input_line_60:4:68: error: unknown type name 'nuisPar' BayesianCalculator bcalc(data, *model, RooArgSet(*POI), *priorPOI, nuisPar); ^ input_line_60:18:27: error: reference to 'data' is ambiguous BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar); ^ input_line_58:2:13: note: candidate found by name lookup is 'data' RooDataSet data("data", "", {*w->var("x"), *w->var("n")}, RooFit::WeightVar("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_60:18:27: error: unknown type name 'data' BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar); ^ input_line_60:18:34: error: a type specifier is required for all declarations BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar); ^ input_line_60:18:59: error: a type specifier is required for all declarations BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar); ^ input_line_60:18:70: error: unknown type name 'nuisPar' BayesianCalculator bcalc2(data, *model, RooArgSet(*POI), *priorPOI2, nuisPar); ^ input_line_60: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
gROOT->GetListOfCanvases()->Draw()