import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)
%jsroot on
filename = "HistFactoryCombinedWorkspace.root"
file = ROOT.TFile.Open("./output/"+filename)
ws = file.Get("HistFactoryWorkspace")
data = ws.data("obsData")
mc = ws.obj("ModelConfig")
pdf = mc.GetPdf()
Here we take the standard Systematic workspace and look at the nominal fit (Maximum Likelihood with Minuit).
fit_result = pdf.fitTo(data,ROOT.RooFit.Save())
Since we will take many different fits this nominal fit with the parameters as they are should be our baseline
minNLL_hat = fit_result.minNll()
ws.saveSnapshot('Start', pdf.getParameters(data))
print minNLL_hat
Now we want to look at what happens when we fit each nuisance parameter to a fixed given value. Let's see which NPs exist in this workspace.
nuisance_params = mc.GetNuisanceParameters()
it = nuisance_params.createIterator()
var = it.Next()
while var:
print var.GetName()
var = it.Next()
The basic idea is to take each NP, this can be selected from the list of all parameters by it's name using the find command. In this first iteration we set it's value to 0. and then fix it there.
np_name = "alpha_background_shape"
nuis = nuisance_params.find(np_name)
ws.loadSnapshot('Start')
nuis.setVal(0.)
nuis.setConstant(1)
Now we look at the minimum of the log likelihood fit here.
result = pdf.fitTo(data,ROOT.RooFit.Save(),ROOT.RooFit.PrintLevel(-1))
result.minNll()
This NP is sort of unique though, it has an interesting range stated as 0-50 however the NLL is poorly behaved outside of 0-2.
nuis.Print("v")
We can however define a function to iterate over various possible values of this NP. Not forgetting to load the snapshot each time to make sure that the parameters are at the state that they were when they were created.
def get_min(ws,np_name,value):
ws.loadSnapshot('Start')
nuis = ws.obj("ModelConfig").GetNuisanceParameters().find(np_name)
nuis.setRange(-50,60)
nuis.setVal(value)
nuis.setConstant(1)
return ws.obj("ModelConfig").GetPdf().fitTo(ws.data("obsData"),ROOT.RooFit.Save(),ROOT.RooFit.PrintLevel(-1))
Let's test with the value that we know the result from.
print get_min(ws,"alpha_background_shape",0.).minNll()
Ok so now let's loop over a load of possible values for this parameter.
import numpy as np
np_values = [0.2*i for i in range(-5,5)] #range(0,2]
d_index = 0
nlls = np.empty(len(np_values))
for value in np_values:
nlls[d_index] = get_min(ws,"alpha_background_shape",value).minNll()
print "for value {0} we get a minNLL of {1}".format(value,nlls[d_index])
d_index+=1
c1 = ROOT.TCanvas()
graph = ROOT.TGraph(len(np_values),np.asarray(np_values),nlls)
graph.SetTitle("alpha_background_shape")
graph.Draw()
c1.Draw()
In this case we want to tranlate the y axis into 2$\Delta$[-log(L)] and draw on the 1 sigma lines.
To form the translation we compute the differenice with the nominal value which we have stored in the "Start" snapshot.
ws.loadSnapshot('Start')
minNLL_hat = fit_result.minNll()
nominal = minNLL_hat
print nominal
It's possible to do it with the graph but since we're only interested in certain values of 2$\Delta$[-log(L)] I want to be able to put this calculation into the function directly.
def get_min(ws,np_name,value,nominal):
ws.loadSnapshot('Start')
nuis = ws.obj("ModelConfig").GetNuisanceParameters().find(np_name)
nuis.setVal(value)
nuis.setConstant(1)
return 2*(ws.obj("ModelConfig").GetPdf().fitTo(ws.data("obsData"),ROOT.RooFit.Save(),ROOT.RooFit.PrintLevel(-1)).minNll() - nominal)
np_values = [0.2*i for i in range(-20,20)] #range(-4,4)
d_index = 0
nlls = np.empty(len(np_values))
for value in np_values:
nlls[d_index] = get_min(ws,"alpha_background_shape",value,nominal)
d_index+=1
c1.Clear()
graph = ROOT.TGraph(len(np_values),np.asarray(np_values),nlls)
graph.SetTitle("alpha_background_shape")
graph.GetYaxis().SetTitle("2#Delta[-log(L)]")
graph.Draw()
line = ROOT.TLine()
line.SetLineStyle(7)
line.SetLineWidth(2)
line.DrawLine(-4, 0, 4, 0)
line.DrawLine(-4, 10, 4, 10)
c1.Draw()
However we're only interested in the bits below 10
c1.Clear()
graph = ROOT.TGraph(len(np_values),np.asarray(np_values),nlls)
graph.SetTitle("alpha_background_shape")
graph.GetYaxis().SetTitle("2#Delta[-log(L)]")
graph.GetYaxis().SetRangeUser(-1,10)
graph.GetXaxis().SetRangeUser(-2,2)
graph.Draw()
line = ROOT.TLine()
line.SetLineStyle(7)
line.SetLineWidth(2)
line.DrawLine(-2, 0, 2, 0)
line.DrawLine(-2, 1, 2, 1)
c1.Draw()