ROOT   Reference Guide
Searching...
No Matches
rs701_BayesianCalculator.py File Reference

## Namespaces

namespace  rs701_BayesianCalculator

## Detailed Description

Bayesian calculator: basic example

Bayesian Result using a Flat prior
0.9% CL central interval: [ 0.506059550587157 - 6.893261312862926 ] or 0.95% CL limits
Bayesian Result using a 1/sqrt(s) prior
0.9% CL central interval: [ 0.07462937796837053 - 5.8542450825237315 ] or 0.95% CL limits
import ROOT
useBkg = True
confLevel = 0.90
w = ROOT.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)")
model = w.factory("PROD::model(pdf,prior_b)") # pdf*priorNuisance
nuisanceParameters = ROOT.RooArgSet(w["b"])
POI = w["s"]
priorPOI = w.factory("Uniform::priorPOI(s)")
priorPOI2 = w.factory("GenericPdf::priorPOI2('1/sqrt(@0)',s)")
w.factory("n[3]") # observed number of events
# create a data set with n observed events
data = ROOT.RooDataSet("data", "", {w["x"], w["n"]}, "n")
# to suppress messages when pdf goes to zero
ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.FATAL)
nuisPar = 0
if useBkg:
nuisPar = nuisanceParameters
else:
w["b"].setVal(0)
size = 1.0 - confLevel
print("\nBayesian Result using a Flat prior ")
bcalc = ROOT.RooStats.BayesianCalculator(data, model, {POI}, priorPOI, nuisPar)
bcalc.SetTestSize(size)
interval = bcalc.GetInterval()
cl = bcalc.ConfidenceLevel()
print(
"{}% CL central interval: [ {} - {} ] or {}% CL limits\n".format(
cl, interval.LowerLimit(), interval.UpperLimit(), cl + (1.0 - cl) / 2
)
)
plot = bcalc.GetPosteriorPlot()
c1 = ROOT.TCanvas("c1", "Bayesian Calculator Result")
c1.Divide(1, 2)
c1.cd(1)
plot.Draw()
c1.Update()
print("\nBayesian Result using a 1/sqrt(s) prior ")
bcalc2 = ROOT.RooStats.BayesianCalculator(data, model, {POI}, priorPOI2, nuisPar)
bcalc2.SetTestSize(size)
interval2 = bcalc2.GetInterval()
cl = bcalc2.ConfidenceLevel()
print(
"{}% CL central interval: [ {} - {} ] or {}% CL limits\n".format(
cl, interval2.LowerLimit(), interval2.UpperLimit(), cl + (1.0 - cl) / 2
)
)
plot2 = bcalc2.GetPosteriorPlot()
c1.cd(2)
plot2.Draw()
c1.Update()
# 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
c1.SaveAs("rs701_BayesianCalculator.png")
# TODO: The BayesianCalculator has to be destructed first. Otherwise, we can get
# segmentation faults depending on the destruction order, which is random in
# Python. Probably the issue is that some object has a non-owning pointer to
# another object, which it uses in its destructor. This should be fixed either
# in the design of RooStats in C++, or with phythonizations.
del bcalc, bcalc2
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t format
Date
July 2022

Definition in file rs701_BayesianCalculator.py.