I’m trying to convert RooStats tutorial “rs101_limitexample.C” to python
but getting
terminate called after throwing an instance of 'std::length_error'
what(): vector::reserve
*** Break *** abort
at the end of the execution, also not getting the proper graph.
Tried looking for a long time but could not able to locate the issue.
My python code for the above tutorial
rs101_limitexample.py
import ROOT
from ROOT import RooStats
# --------------------------------------
# An example of setting a limit in a number counting experiment with uncertainty on background and signal
# to time the macro
t = ROOT.TStopwatch()
t.Start()
# --------------------------------------
# The Model building stage
# --------------------------------------
wspace = ROOT.RooWorkspace()
wspace.factory("Poisson::countingModel(obs[150,0,300], "
"sum(s[50,0,120]*ratioSigEff[1.,0,3.],b[100]*ratioBkgEff[1.,0.,3.]))") # counting model
wspace.factory("Gaussian::sigConstraint(gSigEff[1,0,3],ratioSigEff,0.05)") # 5% signal efficiency uncertainty
wspace.factory("Gaussian::bkgConstraint(gSigBkg[1,0,3],ratioBkgEff,0.2)") # 10% background efficiency uncertainty
wspace.factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)") # product of terms
wspace.Print()
modelWithConstraints = wspace.pdf("modelWithConstraints") # get the model
obs = wspace.var("obs") # get the observable
s = wspace.var("s") # get the signal we care about
b = wspace.var("b") # get the background and set it to a constant. Uncertainty included in ratioBkgEff
b.setConstant()
ratioSigEff = wspace.var("ratioSigEff") # get uncertain parameter to constrain
ratioBkgEff = wspace.var("ratioBkgEff") # get uncertain parameter to constrain
constrainedParams = ROOT.RooArgSet(ratioSigEff,
ratioBkgEff) # need to constrain these in the fit (should change default behavior)
gSigEff = wspace.var("gSigEff") # global observables for signal efficiency
gSigBkg = wspace.var("gSigBkg") # global obs for background efficiency
gSigEff.setConstant()
gSigBkg.setConstant()
# Create an example dataset with 160 observed events
obs.setVal(160.)
data = ROOT.RooDataSet("exampleData", "exampleData", ROOT.RooArgSet(obs))
data.add(ROOT.RooArgSet(obs))
all = ROOT.RooArgSet(s, ratioBkgEff, ratioSigEff)
# not necessary
modelWithConstraints.fitTo(data, ROOT.RooFit.Constrain(ROOT.RooArgSet(ratioSigEff, ratioBkgEff)))
# Now let's make some confidence intervals for s, our parameter of interest
paramOfInterest = ROOT.RooArgSet(s)
modelConfig = RooStats.ModelConfig(wspace)
modelConfig.SetPdf(modelWithConstraints)
modelConfig.SetParametersOfInterest(paramOfInterest)
modelConfig.SetNuisanceParameters(constrainedParams)
modelConfig.SetObservables(obs)
modelConfig.SetGlobalObservables(ROOT.RooArgSet(gSigEff, gSigBkg))
modelConfig.SetName("ModelConfig")
wspace.Import(modelConfig)
wspace.Import(data)
wspace.SetName("w")
wspace.writeToFile("rs101_ws.root")
# First, let's use a Calculator based on the Profile Likelihood Ratio
# ProfileLikelihoodCalculator plc(*data, *modelWithConstraints, paramOfInterest);
plc = RooStats.ProfileLikelihoodCalculator(data, modelConfig)
plc.SetTestSize(.05)
lrinterval = plc.GetInterval() # that was easy.
# Let's make a plot
dataCanvas = ROOT.TCanvas("dataCanvas")
dataCanvas.Divide(2, 1)
dataCanvas.cd(1)
plotInt = RooStats.LikelihoodIntervalPlot(lrinterval)
plotInt.SetTitle("Profile Likelihood Ratio and Posterior for S")
plotInt.Draw()
# Second, use a Calculator based on the Feldman Cousins technique
fc = RooStats.FeldmanCousins(data, modelConfig)
fc.UseAdaptiveSampling(True)
fc.FluctuateNumDataEntries(False) # number counting analysis: dataset always has 1 entry with N events observed
fc.SetNBins(100) #number of points to test per parameter
fc.SetTestSize(.05)
# fc.SaveBeltToFile(true); # optional
# fcint = RooStats.ConfInterval()
fcint = fc.GetInterval() # that was easy.
fit = modelWithConstraints.fitTo(data, ROOT.RooFit.Save(True))
# Third, use a Calculator based on Markov Chain monte carlo
# Before configuring the calculator, let's make a ProposalFunction
# that will achieve a high acceptance rate
ph = RooStats.ProposalHelper()
ph.SetVariables(fit.floatParsFinal())
ph.SetCovMatrix(fit.covarianceMatrix())
ph.SetUpdateProposalParameters(True)
ph.SetCacheSize(100)
pdfProp = ph.GetProposalFunction() # that was easy
mc = RooStats.MCMCCalculator(data, modelConfig)
mc.SetNumIters(20000) # steps to propose in the chain
mc.SetTestSize(.05) # 95% CL
mc.SetNumBurnInSteps(40) # ignore first N steps in chain as "burn in"
mc.SetProposalFunction(pdfProp)
mc.SetLeftSideTailFraction(0.5) # find a "central" interval
mcInt = mc.GetInterval() # that was easy
# Get Lower and Upper limits from Profile Calculator
# print("Profile lower limit on s = {:.4f}" .format(RooStats.LikelihoodInterval(lrinterval).LowerLimit(s)))
# print("Profile upper limit on s = {:.4f}" .format(RooStats.LikelihoodInterval(lrinterval).UpperLimit(s)))
# lrinterval = RooStats.LikelihoodInterval(lrinterval)
print("Profile lower limit on s = {:.4f}" .format(lrinterval.LowerLimit(s)))
print("Profile lower limit on s = {:.4f}" .format(lrinterval.UpperLimit(s)))
# Get Lower and Upper limits from FeldmanCousins with profile construction
if (fcint != None) :
fcul = fcint.UpperLimit(s)
fcll = fcint.LowerLimit(s)
print("FC lower limit on s = {:.1f}".format(fcll))
print("FC upper limit on s = {:.1f}".format(fcul))
fcllLine = ROOT.TLine(fcll, 0, fcll, 1)
fculLine = ROOT.TLine(fcul, 0, fcul, 1)
fcllLine.SetLineColor(ROOT.kRed)
fculLine.SetLineColor(ROOT.kRed)
fcllLine.Draw("same")
fculLine.Draw("same")
dataCanvas.Update()
# Plot MCMC interval and print some statistics
mcPlot = RooStats.MCMCIntervalPlot(mcInt)
mcPlot.SetLineColor(ROOT.kMagenta)
mcPlot.SetLineWidth(2)
mcPlot.Draw("same")
mcul = mcInt.UpperLimit(s)
mcll = mcInt.LowerLimit(s)
print("MCMC lower limit on s = {:.4f}" .format(mcll))
print("MCMC upper limit on s = {:.4f}" .format(mcul))
print("MCMC Actual confidence level: {:.6f}" .format(mcInt.GetActualConfidenceLevel()))
# 3-d plot of the parameter points
dataCanvas.cd(2)
# also plot the points in the markov chain
# mcInt = mcInt = RooStats.MCMCInterval(mc.GetInterval())
chainData = mcInt.GetChainAsDataSet()
print("check 1")
assert chainData
print("plotting the chain data - nentries = {}" .format( chainData.numEntries()) )
chain = RooStats.GetAsTTree("chainTreeData", "chainTreeData", chainData)
assert chain
chain.SetMarkerStyle(6)
chain.SetMarkerColor(ROOT.kRed)
print("check 2")
chain.Draw("s:ratioSigEff:ratioBkgEff", "nll_MarkovChain_local_", "box") # 3-d box proportional to posterior
# the points used in the profile construction
parScanData = ROOT.RooDataSet(fc.GetPointsToScan())
assert parScanData
print("plotting the scanned points used in the frequentist construction - npoints = {}"
.format(parScanData.numEntries()) )
# getting the tree and drawing it -crashes (very strange....);
gr = ROOT.TGraph2D(parScanData.numEntries())
for ievt in range(parScanData.numEntries()):
evt = parScanData.get(ievt)
x = evt.getRealValue("ratioBkgEff")
y = evt.getRealValue("ratioSigEff")
z = evt.getRealValue("s")
gr.SetPoint(ievt, x, y, z)
# std::cout << ievt << " " << x << " " << y << " " << z << std::endl;
print("check 4")
gr.SetMarkerStyle(24)
gr.Draw("P SAME")
del wspace
del lrinterval
del mcInt
del fcint
del data
# print timing info
t.Stop()
t.Print()
Thank you for looking into this issue 
ROOT Version: 6.24/06
Platform: wsl2 using conda