import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)
%jsroot on

from matplotlib import pyplot
%matplotlib inline

Welcome to JupyROOT 6.07/07

Reading Workspaces and Recreating Models

First we read in the signal region model from the workspace building notebook.

infile = ROOT.TFile.Open("output/SRselection.root")
w = infile.Get("w")
mc = w.obj("ModelConfig")
data = w.data("data")
mc.LoadSnapshot()
x = w.var("x")

Bayesian Numerical Calculator

For comparison let's first ask the ProfileLikelihood calculator for the 68% interval

pl = ROOT.RooStats.ProfileLikelihoodCalculator(data,mc)
pl.SetConfidenceLevel(0.683)
PLinterval = pl.GetInterval()
firstPOI = mc.GetParametersOfInterest().first()
lowerLimit = PLinterval.LowerLimit(firstPOI)
upperLimit = PLinterval.UpperLimit(firstPOI)
print "68% interval on {0} is [{1},{2}]".format(firstPOI.GetName(),lowerLimit,upperLimit)

68% interval on n_sig is [300.405555929,589.564440644]

We can now try to use the numerical Bayesian calculator on the same model. As before we create the class and configure it setting for example the integration type or the number of iterations. Since the numerical integration can be tricky, it is better in this case to set before a reasonable range from the nuisance parameters, which will be integrated to obtain the marginalised posterior function. A sensible choice is to use for example an interval using 10-sigma around the best fit values for the nuisance parameters.

Note that we have not defined a prior distribution for the parameter of interest (the number of signal events). If a prior is not defined in the model, it is assumed that the prior distribution for the parameter of interest is the uniform distribution in the parameter range.

define reasanable range for nuisance parameters ( +/- 10 sigma around fitted values) to facilitate the integral

mc.LoadSnapshot()
mc.Print()

=== Using the following for ModelConfig ===
Observables:             RooArgSet:: = (x)
Parameters of Interest:  RooArgSet:: = (n_sig)
Nuisance Parameters:     RooArgSet:: = (n_bkg)
PDF:                     RooAddPdf::model[ n_sig * sigpdf + n_bkg * bkgpdf ] = 26.7114
Snapshot:                
  1) 0x8077c80 RooRealVar:: n_sig = 421.877  L(42.1877 - 4218.77)  "n_sig"

Unfortunately iterating over parameters is really ikky in python. Also pointers are harder to define.

nuisPar = mc.GetNuisanceParameters()
nuis_iter = nuisPar.createIterator()
var = nuis_iter.Next()
while var:
    var.setRange(var.getVal()-10*var.getError(), var.getVal()+10*var.getError())
    var = nuis_iter.Next()

As before, we define our calculator at 68.3% confidence level.

Again the interval can be defined but here central is default. SetLeftSideTailFraction = {0.5: central interval, 0.: upper limit, NULL:shortest}

the intergration can also be defined using integration type. Default is ADAPTIVE but also available are "VEGAS" , "MISER", "PLAIN" and "TOYMC" if using TOYMC you also need the following lines:

if "TOYMC" in integrationType: nuisPdf = ROOT.RooStats.MakeNuisancePdf(mc, "nuisance_pdf") bayesianCalc.ForceNuisancePdf(nuisPdf)

Now we compute interval by scanning the posterior function.

bayesianCalc = ROOT.RooStats.BayesianCalculator(data,mc)
bayesianCalc.SetConfidenceLevel(0.683)
bayesianCalc.SetLeftSideTailFraction(0.5)
bayesianCalc.SetNumIters(1000)
bayesianCalc.SetIntegrationType("") 
bayesianCalc.SetScanOfPosterior(50)

Now we compute interval by scanning the posterior function.

firstPOI = mc.GetParametersOfInterest().first()

interval = bayesianCalc.GetInterval()

lowerLimit = interval.LowerLimit()
upperLimit = interval.UpperLimit()

print "68% interval on {0} is : [{1},{2}]".format(firstPOI.GetName(),lowerLimit,upperLimit)

68% interval on n_sig is : [298.2307904,594.029283082]

c = ROOT.TCanvas()
plot = bayesianCalc.GetPosteriorPlot()
plot.Draw()
c.Draw()