diff --git a/Fitter/scripts/EFTPlotter.py b/Fitter/scripts/EFTPlotter.py index 8576a6a..cbffd2e 100644 --- a/Fitter/scripts/EFTPlotter.py +++ b/Fitter/scripts/EFTPlotter.py @@ -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 = [] @@ -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() @@ -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 + 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') @@ -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") @@ -677,13 +705,10 @@ 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) @@ -691,8 +716,8 @@ def LLPlot2DEFT(self, name='.test', wcs=[], ceiling=1, log=False, final=False): # 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 hname = '{}{}less{}'.format(wcs[0],wcs[1],ceiling) if log: hname += "_log" @@ -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: @@ -924,16 +952,13 @@ 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) @@ -941,8 +966,8 @@ def ContourPlotEFT(self, name='.test', wcs=[], final=False): # 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)) @@ -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) @@ -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!") diff --git a/Fitter/scripts/parse_nll.py b/Fitter/scripts/parse_nll.py index 9c20688..e74023c 100644 --- a/Fitter/scripts/parse_nll.py +++ b/Fitter/scripts/parse_nll.py @@ -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 @@ -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