import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)
%jsroot on
from matplotlib import pyplot
%matplotlib inline
First we read in the signal region model from the workspace building notebook.
infile = ROOT.TFile.Open("output/mH_BinModel.root")
w = infile.Get("w")
mc = w.obj("ModelConfig")
data = w.data("data")
mc.LoadSnapshot()
x = w.var("x")
The aim of this exercise is to compute the significance of observing a signal in the Gaussian plus Exponential model. It is better to re-generate the model using a lower value for the number of signal events (e.g. 500), in order not to have a too high significance, which will then be difficult to estimate if we use pseudo-experiments. To estimate the significance, we need to perform an hypothesis test. We want to disprove the null model, i.e the background only model against the alternate model, the background plus the signal. In RooStats, we do this by defining two two ModelConfig objects, one for the null model (the background only model in this case) and one for the alternate model (the signal plus background).
We have defined in the workspace saved in the ROOT file only the signal plus background model. We can create the background model from the signal plus background model by setting the value of nsig=0. We do this by cloning the S+B model and by defining a snapshot in the ModelConfig for nsig with a different value:
sbModel = w.obj("ModelConfig")
poi = sbModel.GetParametersOfInterest().first()
poi.setVal(500)
sbModel.SetSnapshot(ROOT.RooArgSet(poi))
bModel = sbModel.Clone()
bModel.SetName("B_only_model")
poi.setVal(0)
bModel.SetSnapshot(ROOT.RooArgSet(poi))
Using the given models plus the observed data set, we can create a RooStats hypothesis test calculator to compute the significance. We have the choice between
The first two require generate pseudo-experiment, while the Asymptotic calculator, requires only a fit to the two models. For the Frequentist and the Hybrid calculators we need to choose also the test statistics. We have various choices in RooStats. Since it is faster, we will use for the exercise the Asymptotic calculator. The Asymptotic calculator it is based on the Profile Likelihood test statistics, the one used at LHC (see slides for its definition). For testing the significance, we must use the one-side test statistic (we need to configure this explicitly in the class). Summarising, we will do:
infile = ROOT.TFile.Open("output/SRselection.root")
w = infile.Get("w")
data = w.data("data")
sbModel = w.obj("ModelConfig")
sbModel.SetName("S+B Model")
poi = sbModel.GetParametersOfInterest().first()
poi.setVal(500)
sbModel.SetSnapshot(ROOT.RooArgSet(poi))
bModel = sbModel.Clone()
bModel.SetName("B Model")
poi.setVal(0)
bModel.SetSnapshot(ROOT.RooArgSet(poi))
w.Print()
ac = ROOT.RooStats.AsymptoticCalculator(data, sbModel, bModel)
ac.SetOneSidedDiscovery(True)
asResult = ac.GetHypoTest()
asResult.Print()
And the frequentist appproach.
This needs test statistics.
fc = ROOT.RooStats.FrequentistCalculator(data, sbModel, bModel)
fc.SetToys(2000,500)
profll = ROOT.RooStats.ProfileLikelihoodTestStat(sbModel.GetPdf())
profll.SetOneSidedDiscovery(True)
toymcs = ROOT.RooStats.ToyMCSampler(fc.GetTestStatSampler())
toymcs.SetTestStatistic(profll)
if not sbModel.GetPdf().canBeExtended():
toymcs.SetNEventsPerToy(1)
print 'adjusting for non-extended formalism'
Run the test!
fqResult = fc.GetHypoTest()
fqResult.Print()
and plot!
c = ROOT.TCanvas()
plot = ROOT.RooStats.HypoTestPlot(fqResult)
plot.SetLogYaxis(True)
plot.Draw()
c.Draw()