Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 48 additions & 18 deletions Fitter/scripts/EFTPlotter.py
Original file line number Diff line number Diff line change
Expand Up @@ -152,9 +152,9 @@ def AppendStrToItemsInLst(self,in_lst,in_str):


# Takes as input the name of a root file (assumed to be in ../fit_files)
# Retruns [wc vals in the scan, delta nll vals at each point]
# Returns [wc vals in the scan, delta nll vals at each point]
# Optionally removes duplicate wc points (choosing min nll)
def GetWCsNLLFromRoot(self,base_name_lst,wc,unique=False):
def GetWCsNLLFromRoot(self,base_name_lst,wc_lst,unique=False):

graphwcs = []
graphnlls = []
Expand All @@ -170,7 +170,13 @@ def GetWCsNLLFromRoot(self,base_name_lst,wc,unique=False):
# Get coordinates for TGraph
for entry in range(limitTree.GetEntries()):
limitTree.GetEntry(entry)
graphwcs.append(limitTree.GetLeaf(wc).GetValue(0))
if len(wc_lst) == 1:
graphwcs.append(limitTree.GetLeaf(wc_lst[0]).GetValue(0))
elif len(wc_lst) == 2:
graphwcs.append([limitTree.GetLeaf(wc_lst[0]).GetValue(0),limitTree.GetLeaf(wc_lst[1]).GetValue(0)])
else:
logging.error("Cannot handle the case with more than 2 WCs!")

graphnlls.append(2*limitTree.GetLeaf('deltaNLL').GetValue(0))

rootFile.Close()
Expand All @@ -185,6 +191,28 @@ def GetWCsNLLFromRoot(self,base_name_lst,wc,unique=False):
return [graphwcs,graphnlls]


def CreateNewLimitTreefor2DScan(self,name_lst,wc_lst):
graphwcs, graphnlls = self.GetWCsNLLFromRoot(name_lst,wc_lst,unique=True)
newRootFile = ROOT.TFile("tmp.root","RECREATE") #temporary root file with new limit tree
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if it would be more efficient to pass around the tree instead of making a temporary file. Don't worry about doing that in this PR, but maybe just something to think about.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a good idea. I think my logic back then was that this function would allow us to dump out a root file (if need be), and we could investigate the TTree inside this file if we wanted to. Since this was 2 years ago, I definitely didn't have as good of an understanding of the code as I have now, so I would probably have opted for a different way.

limit = ROOT.TTree("limit","limit")
wc0 = numpy.array([0],dtype='f')
wc1 = numpy.array([0],dtype='f')
deltanll = numpy.array([0],dtype='f')
#Create new branches with two WCs being scanned and deltaNLL
limit.Branch("{}".format(wc_lst[0]),wc0,"{}".format(wc_lst[0]))
limit.Branch("{}".format(wc_lst[1]),wc1,"{}".format(wc_lst[1]))
limit.Branch("deltaNLL",deltanll,"deltaNLL")

for idx in range(len(graphnlls)):
wc0[0] = graphwcs[idx][0]
wc1[0] = graphwcs[idx][1]
deltanll[0] = graphnlls[idx]
limit.Fill()
limit.Write()

return newRootFile



def ResetHistoFile(self, name=''):
ROOT.TFile('Histos{}.root'.format(name),'RECREATE')
Expand All @@ -200,7 +228,7 @@ def LLPlot1DEFT(self, name_lst=['.test'], frozen=False, wc='', log=False):
ROOT.gROOT.SetBatch(True)
canvas = ROOT.TCanvas()

graphwcs, graphnlls = self.GetWCsNLLFromRoot(name_lst,wc,unique=True)
graphwcs, graphnlls = self.GetWCsNLLFromRoot(name_lst,[wc],unique=True)
if graphwcs == [] or graphnlls == []:
# Something went wrong
print("Error, probably could not find root file")
Expand Down Expand Up @@ -677,22 +705,19 @@ def BatchOverlayZoomLLPlot1DEFT(self, basename1='.EFT.SM.Float', basename2='.EFT
for wc in wcs:
self.OverlayZoomLLPlot1DEFT(basename1+'.'+wc, basename2+'.'+wc, wc, log)

def LLPlot2DEFT(self, name='.test', wcs=[], ceiling=1, log=False, final=False):
def LLPlot2DEFT(self, name_lst=['.test'], wcs=[], ceiling=1, log=False, final=False):
if len(wcs)!=2:
logging.error("Function 'LLPlot2D' requires exactly two wcs!")
return
if not os.path.exists('../fit_files/higgsCombine{}.MultiDimFit.root'.format(name)):
logging.error("File higgsCombine{}.MultiDimFit.root does not exist!".format(name))
return

ROOT.gROOT.SetBatch(True)

canvas = ROOT.TCanvas('c','c',800,800)

# Open file and draw 2D histogram
# wcs[0] is y-axis variable, wcs[1] is x-axis variable
rootFile = ROOT.TFile.Open('../fit_files/higgsCombine{}.MultiDimFit.root'.format(name))
limitTree = rootFile.Get('limit')
rootFile = self.CreateNewLimitTreefor2DScan(name_lst,wcs)
limitTree = rootFile.limit
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line actually works? I would have expected it to be

limitTree = rootFile.Get('limit')

hname = '{}{}less{}'.format(wcs[0],wcs[1],ceiling)
if log:
hname += "_log"
Expand Down Expand Up @@ -752,7 +777,10 @@ def LLPlot2DEFT(self, name='.test', wcs=[], ceiling=1, log=False, final=False):
hist.Write()
outfile.Close()
hist.SetDirectory(0)
#return hist

#close the tmp root file
rootFile.Close()
os.remove("tmp.root")

'''
Example:
Expand Down Expand Up @@ -924,25 +952,22 @@ def LLPlot2DVarEFT(self, names=['.test'], wcs=[], ceiling=1, log=False, final=Fa
rootFile.Close()
os.system('rm ../fit_files/higgsCombine.var{}.MultiDimFit.root'.format(names[0]))

def ContourPlotEFT(self, name='.test', wcs=[], final=False):
def ContourPlotEFT(self, name_lst=['.test'], wcs=[], final=False):
#hist2d = self.LLPlot2DEFT(name,wcs,9,False)
#hist = ROOT.TH2F('hist', hname, points, self.wc_ranges[wcs[1]][0], self.wc_ranges[wcs[1]][1], points, self.wc_ranges[wcs[0]][0], self.wc_ranges[wcs[0]][1])
#limitTree.Draw('2*(deltaNLL-{}):{}:{}>>hist({},{},{},{},{},{})'.format(minZ,wcs[1],wcs[0],points,self.wc_ranges[wcs[0]][0],self.wc_ranges[wcs[0]][1],points,self.wc_ranges[wcs[1]][0],self.wc_ranges[wcs[1]][1]), '2*(deltaNLL-{})<{}'.format(minZ,ceiling), 'prof colz')
if len(wcs)!=2:
logging.error("Function 'ContourPlot' requires exactly two wcs!")
return
if not os.path.exists('../fit_files/higgsCombine{}.MultiDimFit.root'.format(name)):
logging.error("File higgsCombine{}.MultiDimFit.root does not exist!".format(name))
return

best2DeltaNLL = 1000000
ROOT.gROOT.SetBatch(True)
canvas = ROOT.TCanvas('c','c',800,800)

# Get Grid scan and copy to h_contour
# wcs[0] is y-axis variable, wcs[1] is x-axis variable
gridFile = ROOT.TFile.Open('../fit_files/higgsCombine{}.MultiDimFit.root'.format(name))
gridTree = gridFile.Get('limit')
gridFile = self.CreateNewLimitTreefor2DScan(name_lst,wcs)
gridTree = gridFile.limit
minZ = gridTree.GetMinimum('deltaNLL')
points = 100
#gridTree.Draw('2*(deltaNLL-{}):{}:{}>>+hist'.format(minZ,wcs[0],wcs[1]), '2*(deltaNLL-{})<{}'.format(minZ,9))
Expand Down Expand Up @@ -1042,7 +1067,8 @@ def ContourPlotEFT(self, name='.test', wcs=[], final=False):
self.CMS_extra.SetTextFont(52)
if not final: self.CMS_extra.Draw('same')
scan_name = 'Other WCs profiled'
if 'Froz' in name or 'Freeze' in name or 'frozen' in name:

if 'Froz' in [name for name in name_lst] or 'Freeze' in [name for name in name_lst] or 'frozen' in [name for name in name_lst]:
scan_name = 'Other WCs fixed to SM'
self.scan_type = ROOT.TLatex(0.15, 0.885, scan_name)
self.scan_type.SetNDC(1)
Expand Down Expand Up @@ -1168,6 +1194,10 @@ def ContourPlotEFT(self, name='.test', wcs=[], final=False):

ROOT.gStyle.SetPalette(57)

#close the tmp root file
gridFile.Close()
os.remove("tmp.root")

def LLPlot1DSM(self, name='.test', param='', log=False):
if not param:
logging.error("No parameter specified!")
Expand Down
4 changes: 3 additions & 1 deletion Fitter/scripts/parse_nll.py
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,8 @@ def get_unique_points(in_dict,scan_var,minimize_var):
scan_var_val_lst_unique = {}
for idx in range(ref_len):
scan_var_val = in_dict[scan_var][idx]
if any(isinstance(i, list) for i in in_dict[scan_var]): #relevant only for 2D scans
scan_var_val = tuple(scan_var_val)
minimize_var_val = in_dict[minimize_var][idx]
if scan_var_val not in scan_var_val_lst_unique:
# This x value is not yet in our unique list
Expand All @@ -137,7 +139,7 @@ def get_unique_points(in_dict,scan_var,minimize_var):
idx_to_keep = np.array(idx_to_keep)
for var_name in in_dict.keys():
var_val_arr = np.array(in_dict[var_name])
var_val_arr_unique = np.take(var_val_arr,idx_to_keep)
var_val_arr_unique = var_val_arr[idx_to_keep]
out_dict[var_name] = list(var_val_arr_unique)

return out_dict
Expand Down