diff --git a/fillHistosScouting.py b/fillHistosScouting.py index 92efe53..af39d9d 100644 --- a/fillHistosScouting.py +++ b/fillHistosScouting.py @@ -265,7 +265,6 @@ def getweight(era, ngen, frac=1.0, xsec=1000): fin.close() os.system('rm -f filein.txt') print("Found {} files matching criteria".format(len(files))) -print(files) index = int(args.splitIndex) pace = int(args.splitPace) @@ -1195,6 +1194,8 @@ def getweight(era, ngen, frac=1.0, xsec=1000): if not applyLzSelection(lz): continue mass = v.M() + uv = dmu_muvecdp[-1] + dmu_muvecdp[-2] + umass = uv.M() pt = v.Pt() nhitsbeforesvtotal = t.Muon_nhitsbeforesv[dmuidxs[int(vn*2)]] + t.Muon_nhitsbeforesv[dmuidxs[int(vn*2)+1]] # @@ -1277,6 +1278,15 @@ def getweight(era, ngen, frac=1.0, xsec=1000): dphisv2u = abs(dmu_muvecdp[int(vn*2)+1].Vect().DeltaPhi(svvec[vn])) detasvu = abs(vu.Vect().Eta()-svvec[vn].Eta()) detadphisvu = 1e6 + delta_phi = dphisv - dphisvu + if delta_phi>0: + addition_seagull = dmu_muvecdp[-1] + dmu_muvecdp[-2] + seagullumass = addition_seagull.M() + cowboyumass = 0 + elif delta_phi<0: + addition_cowboy = dmu_muvecdp[-1] + dmu_muvecdp[-2] + cowboyumass = addition_cowboy.M() + seagullumass = 0 if dphisvu>0.0: detadphisvu = detasvu/dphisvu a3dsvu = abs(vu.Vect().Angle(svvec[vn])) diff --git a/python/make_efficiencies.py b/python/make_efficiencies.py new file mode 100644 index 0000000..e9e7893 --- /dev/null +++ b/python/make_efficiencies.py @@ -0,0 +1,322 @@ +import ROOT +import os,sys,json +import argparse +from datetime import date +import numpy as np +import copy +import math +sys.path.append('utils') +import plotUtils +from IPython.display import Image +import argparse +import os +import mplhep as hep +import matplotlib.pyplot as plt +ROOT.gStyle.SetOptStat(0) +parser = argparse.ArgumentParser() +parser.add_argument('directory', type=str) +args = parser.parse_args() +files = [] +directory = args.directory + +latex = ROOT.TLatex() +latex.SetTextFont(42) +latex.SetTextAlign(31) +latex.SetTextSize(0.04) +latex.SetNDC(True) + +if not os.path.exists('cuts_plots/'): + os.mkdir('cuts_plots/') + +def getValues(histo): + values = [] + bins = [] + for n in range(1,histo.GetNbinsX()+1): + values.append(histo.GetBinContent(n)) + bins.append(histo.GetBinLowEdge(n)) + bins.append(histo.GetBinLowEdge(n) + histo.GetBinWidth(n)) + return np.array(values), np.array(bins) + +sum_100mm_allCuts = {} +sum_10mm_allCuts = {} +sum_1mm_allCuts = {} +sum_100mm_noCuts = {} +sum_10mm_noCuts = {} +sum_1mm_noCuts = {} +sum_100mm_Cut1 = {} +sum_10mm_Cut1 = {} +sum_1mm_Cut1 = {} +sum_100mm_Cut2 = {} +sum_10mm_Cut2 = {} +sum_1mm_Cut2 = {} +sum_100mm_Cut3 = {} +sum_10mm_Cut3 = {} +sum_1mm_Cut3 = {} +sum_100mm_Cut4 = {} +sum_10mm_Cut4 = {} +sum_1mm_Cut4 = {} +for subdir in os.listdir(directory): + if 'Jun-25' in subdir: + full_subdir_path = os.path.join(directory, subdir, '') + files = os.listdir(full_subdir_path) + for filename in files: + if '_all' in filename: + continue + if '-2p0_' in filename: + mass = 2 + elif '-7p0_' in filename: + mass = 7 + elif '-10p0_' in filename: + mass = 10 + elif '-20p0_' in filename: + mass = 20 + elif '-30p0_' in filename: + mass = 30 + else: + continue + file = ROOT.TFile.Open(full_subdir_path + filename) + if 'allCuts' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_allCuts = hdimuon_mass.Integral() + current_value_100mm_allCuts = sum_100mm_allCuts.get(mass, 0) + updated_value_100mm_allCuts = current_value_100mm_allCuts + n_hist_100mm_allCuts + sum_100mm_allCuts[mass] = updated_value_100mm_allCuts + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_allCuts = hdimuon_mass.Integral() + current_value_10mm_allCuts = sum_10mm_allCuts.get(mass, 0) + updated_value_10mm_allCuts = current_value_10mm_allCuts + n_hist_10mm_allCuts + sum_10mm_allCuts[mass] = updated_value_10mm_allCuts + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_allCuts = hdimuon_mass.Integral() + current_value_1mm_allCuts = sum_1mm_allCuts.get(mass, 0) + updated_value_1mm_allCuts = current_value_1mm_allCuts + n_hist_1mm_allCuts + sum_1mm_allCuts[mass] = updated_value_1mm_allCuts + + + elif 'noCuts' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_noCuts = hdimuon_mass.Integral() + current_value_100mm_noCuts = sum_100mm_noCuts.get(mass, 0) + updated_value_100mm_noCuts = current_value_100mm_noCuts + n_hist_100mm_noCuts + sum_100mm_noCuts[mass] = updated_value_100mm_noCuts + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_noCuts = hdimuon_mass.Integral() + current_value_10mm_noCuts = sum_10mm_noCuts.get(mass, 0) + updated_value_10mm_noCuts = current_value_10mm_noCuts + n_hist_10mm_noCuts + sum_10mm_noCuts[mass] = updated_value_10mm_noCuts + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_noCuts = hdimuon_mass.Integral() + current_value_1mm_noCuts = sum_1mm_noCuts.get(mass, 0) + updated_value_1mm_noCuts = current_value_1mm_noCuts + n_hist_1mm_noCuts + sum_1mm_noCuts[mass] = updated_value_1mm_noCuts + elif 'DimuonAngularSel' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut1 = hdimuon_mass.Integral() + current_value_100mm_Cut1 = sum_100mm_Cut1.get(mass, 0) + updated_value_100mm_Cut1 = current_value_100mm_Cut1 + n_hist_100mm_Cut1 + sum_100mm_Cut1[mass] = updated_value_100mm_Cut1 + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut1 = hdimuon_mass.Integral() + current_value_10mm_Cut1 = sum_10mm_Cut1.get(mass, 0) + updated_value_10mm_Cut1 = current_value_10mm_Cut1 + n_hist_10mm_Cut1 + sum_10mm_Cut1[mass] = updated_value_10mm_Cut1 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut1 = hdimuon_mass.Integral() + current_value_1mm_Cut1 = sum_1mm_Cut1.get(mass, 0) + updated_value_1mm_Cut1 = current_value_1mm_Cut1 + n_hist_1mm_Cut1 + sum_1mm_Cut1[mass] = updated_value_1mm_Cut1 + + elif 'MaterialVeto' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut2 = hdimuon_mass.Integral() + current_value_100mm_Cut2 = sum_100mm_Cut2.get(mass, 0) + updated_value_100mm_Cut2 = current_value_100mm_Cut2 + n_hist_100mm_Cut2 + sum_100mm_Cut2[mass] = updated_value_100mm_Cut2 + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut2 = hdimuon_mass.Integral() + current_value_10mm_Cut2 = sum_10mm_Cut2.get(mass, 0) + updated_value_10mm_Cut2 = current_value_10mm_Cut2 + n_hist_10mm_Cut2 + sum_10mm_Cut2[mass] = updated_value_10mm_Cut2 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut2 = hdimuon_mass.Integral() + current_value_1mm_Cut2 = sum_1mm_Cut2.get(mass, 0) + updated_value_1mm_Cut2 = current_value_1mm_Cut2 + n_hist_1mm_Cut2 + sum_1mm_Cut2[mass] = updated_value_1mm_Cut2 + + + elif 'MuonHitSel' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut3 = hdimuon_mass.Integral() + current_value_100mm_Cut3 = sum_100mm_Cut3.get(mass, 0) + updated_value_100mm_Cut3 = current_value_100mm_Cut3 + n_hist_100mm_Cut3 + sum_100mm_Cut3[mass] = updated_value_100mm_Cut3 + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut3 = hdimuon_mass.Integral() + current_value_10mm_Cut3 = sum_10mm_Cut3.get(mass, 0) + updated_value_10mm_Cut3 = current_value_10mm_Cut3 + n_hist_10mm_Cut3 + sum_10mm_Cut3[mass] = updated_value_10mm_Cut3 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut3 = hdimuon_mass.Integral() + current_value_1mm_Cut3 = sum_1mm_Cut3.get(mass, 0) + updated_value_1mm_Cut3 = current_value_1mm_Cut3 + n_hist_1mm_Cut3 + sum_1mm_Cut3[mass] = updated_value_1mm_Cut3 + + elif 'MuonIPSel' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut4 = hdimuon_mass.Integral() + current_value_100mm_Cut4 = sum_100mm_Cut4.get(mass, 0) + updated_value_100mm_Cut4 = current_value_100mm_Cut4 + n_hist_100mm_Cut4 + sum_100mm_Cut4[mass] = updated_value_100mm_Cut4 + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut4 = hdimuon_mass.Integral() + current_value_10mm_Cut4 = sum_10mm_Cut4.get(mass, 0) + updated_value_10mm_Cut4 = current_value_10mm_Cut4 + n_hist_10mm_Cut4 + sum_10mm_Cut4[mass] = updated_value_10mm_Cut4 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut4 = hdimuon_mass.Integral() + current_value_1mm_Cut4 = sum_1mm_Cut4.get(mass, 0) + updated_value_1mm_Cut4 = current_value_1mm_Cut4 + n_hist_1mm_Cut4 + sum_1mm_Cut4[mass] = updated_value_1mm_Cut4 + + +masses = [2,7,10,20,30] +ctaus =['1mm','10mm','100mm'] + +nbins = len(masses)*len(ctaus) +maxx = nbins +minx = 0 + +h_noCuts = ROOT.TH1F('h_noCuts', 'No cuts', nbins, minx, maxx) +h_allCuts = ROOT.TH1F('h_allCuts', 'All cuts', nbins, minx, maxx) +h_Cut1 = ROOT.TH1F('h_Cut1', 'Cut 1', nbins, minx, maxx) +h_Cut2 = ROOT.TH1F('h_Cut2', 'Cut 2', nbins, minx, maxx) +h_Cut3 = ROOT.TH1F('h_Cut3', 'Cut 3', nbins, minx, maxx) +h_Cut4 = ROOT.TH1F('h_Cut4', 'Cut 4', nbins, minx, maxx) +for massindex, mass in enumerate(masses): + h_noCuts.Fill(massindex,sum_1mm_noCuts[mass]/sum_1mm_noCuts[mass]) + h_Cut1.Fill(massindex,sum_1mm_Cut1[mass]/sum_1mm_noCuts[mass]) + h_Cut2.Fill(massindex,sum_1mm_Cut2[mass]/sum_1mm_noCuts[mass]) + h_Cut3.Fill(massindex,sum_1mm_Cut3[mass]/sum_1mm_noCuts[mass]) + h_Cut4.Fill(massindex,sum_1mm_Cut4[mass]/sum_1mm_noCuts[mass]) + h_allCuts.Fill(massindex,sum_1mm_allCuts[mass]/sum_1mm_noCuts[mass]) + + h_noCuts.Fill(massindex + len(masses),sum_10mm_noCuts[mass]/sum_10mm_noCuts[mass]) + h_Cut1.Fill(massindex + len(masses),sum_10mm_Cut1[mass]/sum_10mm_noCuts[mass]) + h_Cut2.Fill(massindex + len(masses),sum_10mm_Cut2[mass]/sum_10mm_noCuts[mass]) + h_Cut3.Fill(massindex + len(masses),sum_10mm_Cut3[mass]/sum_10mm_noCuts[mass]) + h_Cut4.Fill(massindex + len(masses),sum_10mm_Cut4[mass]/sum_10mm_noCuts[mass]) + h_allCuts.Fill(massindex + len(masses),sum_10mm_allCuts[mass]/sum_10mm_noCuts[mass]) + + h_noCuts.Fill(massindex + 2*len(masses),sum_100mm_noCuts[mass]/sum_100mm_noCuts[mass]) + h_Cut1.Fill(massindex + 2*len(masses),sum_100mm_Cut1[mass]/sum_100mm_noCuts[mass]) + h_Cut2.Fill(massindex + 2*len(masses),sum_100mm_Cut2[mass]/sum_100mm_noCuts[mass]) + h_Cut3.Fill(massindex + 2*len(masses),sum_100mm_Cut3[mass]/sum_100mm_noCuts[mass]) + h_Cut4.Fill(massindex + 2*len(masses),sum_100mm_Cut4[mass]/sum_100mm_noCuts[mass]) + h_allCuts.Fill(massindex + 2*len(masses),sum_100mm_allCuts[mass]/sum_100mm_noCuts[mass]) + +values_noCuts, bins_noCuts = getValues(h_noCuts) +values_allCuts, bins_allCuts = getValues(h_allCuts) +values_Cut1, bins_Cut1 = getValues(h_Cut1) +values_Cut2, bins_Cut2 = getValues(h_Cut2) +values_Cut3, bins_Cut3 = getValues(h_Cut3) +values_Cut4, bins_Cut4 = getValues(h_Cut4) + + +hist_array = [values_noCuts, values_Cut1, values_Cut2, values_Cut3, values_Cut4, values_allCuts] +label_array = ['No cuts','Dimuon Angular Selection','Material Veto','Muon Hit Selection','Muon IP Selection','All cuts'] +yearenergy="(2022, 13.6 TeV)" + +bin_labels = [] +j = -1 +masses_plot = masses*3 +''' +for i, mass in enumerate(masses_plot): + if i%len(masses) == 0: + j += 1 + bin_labels.append(f"m={mass},\n c$\\tau$={ctaus[j]}") +bin_labels.append('') +''' +hep.style.use("CMS") +fig, ax = plt.subplots(figsize=(10, 7.5)) +hep.histplot( + hist_array, + stack = False, + bins = bins_allCuts, + histtype = "step", + alpha = 1, + label = label_array, + ax = ax +) + + +legend_coords = (0.78,0.4) +ax.set_xlim(0, 15) +expoffset = 0.065 +hep.cms.label(rlabel="") + +ax.set_xticks(bins_allCuts) +j = -1 +for i in range(len(masses_plot)): + if i%len(masses) == 0: + j+=1 + x = (bins_allCuts[i] + bins_allCuts[i+1]) / 2 - 0.45 + y = -0.07 + mass_label = f"{masses_plot[i]} GeV" + ctau_label = f"{ctaus[j]}" + plt.text(x, y, mass_label + "\n" + ctau_label, fontsize=10) + +ax.set_xticks(bins_allCuts) +ax.set_xticklabels([''] * len(bins_allCuts)) +ax.legend(loc='center', bbox_to_anchor=legend_coords, fontsize=15) +ax.text(1.01, 1 + expoffset, yearenergy, transform=ax.transAxes, verticalalignment='top', horizontalalignment='right', fontsize=25) +plt.ylabel("Efficiency") +plt.savefig("cuts_plots/efficiencies_goodplots.png") +plt.close() + + diff --git a/python/make_histogramComparison.py b/python/make_histogramComparison.py new file mode 100644 index 0000000..ce868ab --- /dev/null +++ b/python/make_histogramComparison.py @@ -0,0 +1,257 @@ +import ROOT +import os,sys,json +import argparse +from datetime import date +import numpy as np +import copy +import math +sys.path.append('utils') +#import plotUtils +from IPython.display import Image +import argparse +import os + +ROOT.gStyle.SetOptStat(0) +ROOT.gROOT.SetBatch(1) + +user = os.environ.get("USER") +#today = date.today().strftime("%b-%d-%Y") + +parser = argparse.ArgumentParser() +parser.add_argument('directory', type=str) +args = parser.parse_args() +files = [] +directory = args.directory +indir = directory +hname = 'histograms' +year = '2022' +luminosity = 3.51 +# +latex = ROOT.TLatex() +latex.SetTextFont(42) +latex.SetTextAlign(31) +latex.SetTextSize(0.04) +latex.SetNDC(True) +# +latexCMS = ROOT.TLatex() +latexCMS.SetTextFont(61) +latexCMS.SetTextSize(0.055) +latexCMS.SetNDC(True) +# +latexCMSExtra = ROOT.TLatex() +latexCMSExtra.SetTextFont(52) +latexCMSExtra.SetTextSize(0.04) +latexCMSExtra.SetNDC(True) +# +latexExtra = ROOT.TLatex() +latexExtra.SetTextFont(42) +latexExtra.SetTextSize(0.034) +latexExtra.SetNDC(True) +# +latexExtraBold = ROOT.TLatex() +latexExtraBold.SetTextFont(62) +latexExtraBold.SetTextSize(0.034) +latexExtraBold.SetNDC(True) +# +cmsExtra = "Preliminary" +yearenergy = "%.2f fb^{-1} (%s, 13.6 TeV)"%(luminosity,year) +#directory = '/ceph/cms/store/user/garciaja/Run3ScoutingOutput/outputHistograms_Jun-06-2024_JPsi_lxy0p0to0p2/' +#directory = '/home/users/garciaja/run3Scouting/run3_scouting/outputHistograms_Jun-06-2024/' + +for subdir in os.listdir(directory): + if 'Jul-02' in subdir: + full_subdir_path = os.path.join(directory, subdir, '') + files = os.listdir(full_subdir_path) + data_files = [] + summed_histogram_mass = None + summed_histogram_umass = None + summed_histogram_cowboy = None + summed_histogram_seagull = None + + for filename in files: + if 'Data' in filename: + data_files.append(filename) + + for filename in data_files: + #for filename in files: + file = ROOT.TFile.Open(full_subdir_path + filename) + hdimuon_mass = file.Get("hdimuon_mass") + hdimuon_umass = file.Get("hdimuon_umass") + hdimuon_cowboy = file.Get("hdimuon_cowboy") + hdimuon_seagull = file.Get("hdimuon_seagull") + if summed_histogram_mass is None: + summed_histogram_mass = hdimuon_mass.Clone("summed_histogram_mass") + summed_histogram_mass.SetDirectory(0) + else: + summed_histogram_mass.Add(hdimuon_mass) + if summed_histogram_umass is None: + summed_histogram_umass = hdimuon_umass.Clone("summed_histogram_umass") + summed_histogram_umass.SetDirectory(0) + else: + summed_histogram_umass.Add(hdimuon_umass) + if summed_histogram_cowboy is None: + summed_histogram_cowboy = hdimuon_cowboy.Clone("summed_histogram_cowboy") + summed_histogram_cowboy.SetDirectory(0) + else: + summed_histogram_cowboy.Add(hdimuon_cowboy) + if summed_histogram_seagull is None: + summed_histogram_seagull = hdimuon_seagull.Clone("summed_histogram_seagull") + summed_histogram_seagull.SetDirectory(0) + else: + summed_histogram_seagull.Add(hdimuon_seagull) + + if '0p0to0p2' in subdir: + ymax = 220000 + title = 'lxy = 0.0 cm to 0.2 cm' + name = '0p0to0p2' + xCMS = 0.2 + if '0p2to1p0' in subdir: + ymax = 130000 + title = 'lxy = 0.2 cm to 1.0 cm' + name = '0p2to1p0' + xCMS = 0.2 + if '1p0to2p4' in subdir: + ymax = 10000 + title = 'lxy = 1.0 cm to 2.4 cm' + name = '1p0to2p4' + xCMS = 0.115 + if '2p4to3p1' in subdir: + ymax = 300 + title = 'lxy = 2.4 cm to 3.1 cm' + name = '2p4to3p1' + xCMS = 0.115 + if '3p1to7p0' in subdir: + ymax = 130 + title = 'lxy = 3.1 cm to 7.0 cm' + name = '3p1to7p0' + xCMS = 0.115 + if '3p1to11p0' in subdir: + ymax = 130 + title = 'lxy = 3.1 cm to 11.0 cm' + name = '3p1to11p0' + xCMS = 0.115 + if '7p0to11p0' in subdir: + ymax = 7 + title = 'lxy = 7.0 cm to 11.0 cm' + name = '7p0to11p0' + xCMS = 0.115 + if '11p0to16p0' in subdir: + ymax = 5 + title = 'lxy = 11.0 cm to 16.0 cm' + name = '11p0to16p0' + xCMS = 0.115 + if '16p0to70p0' in subdir: + ymax = 20 + title = 'lxy = 16.0 cm to 70.0 cm' + name = '16p0to70p0' + xCMS = 0.115 + + + if not os.path.exists('JPsiPlots/'): + os.mkdir('JPsiPlots/') + + c1 = ROOT.TCanvas('c1', '', 600, 500) + c1.cd() + summed_histogram_mass.GetXaxis().SetRangeUser(2.6,3.6) + summed_histogram_mass.GetYaxis().SetRangeUser(0, ymax) + summed_histogram_mass.SetLineColor(ROOT.kBlack) + #summed_histogram_mass.SetTitle(title) + #c1.SetLogy(True) + summed_histogram_mass.Draw('HIST') + summed_histogram_umass.SetLineColor(ROOT.kRed) + summed_histogram_umass.Draw('HIST SAME') + summed_histogram_seagull.SetLineColor(ROOT.kCyan) + summed_histogram_seagull.Draw('HIST SAME') + summed_histogram_cowboy.SetLineColor(ROOT.kOrange) + summed_histogram_cowboy.Draw('HIST SAME') + #ROOT.gStyle.SetOptTitle(0) + #latex = ROOT.TLatex() + #latex.SetTextSize(0.04) + #latex.DrawLatexNDC(0.38, 0.92, title) + ''' + box = ROOT.TBox (3.27, ymax*0.78, 3.59, ymax*0.85) + box.SetFillStyle(0) + box.SetLineColor(1) + box.SetLineWidth(1) + box.Draw() + text = ROOT.TLatex(3.28,ymax*0.8, title) + text.SetTextSize(0.03) + text.SetTextFont(42) + text.Draw() + + ''' + legend = ROOT.TLegend(0.15, 0.6, 0.35, 0.8) + legend.SetTextSize(0.04) + legend.SetBorderSize(0) + legend.AddEntry(summed_histogram_mass, "Corrected mass", "l") + legend.AddEntry(summed_histogram_umass, "Uncorrected mass", "l") + legend.AddEntry(summed_histogram_seagull, "Seagull", "l") + legend.AddEntry(summed_histogram_cowboy, "Cowboy", "l") + legend.Draw() + latexCMS.DrawLatex(xCMS,0.91,"CMS"); + latexCMSExtra.DrawLatex(xCMS + 0.12,0.91, cmsExtra); + latex = ROOT.TLatex() + latex.SetTextSize(0.04) + latexExtra.DrawLatex(0.62, 0.91, yearenergy); + + + output_filename = f'JPsiPlots/lxy_{name}.png' + c1.SaveAs(output_filename) + Image(filename=output_filename) + #file.close() +''' + + if '-6p0_ctau-' in filename: + file = ROOT.TFile.Open(directory + filename) + #file = ROOT.TFile.Open(filename) + #if '100mm' in filename: + # file.ls() + + hdimuon_mass = file.Get("hdimuon_mass") + hdimuon_umass = file.Get("hdimuon_umass") + if not os.path.exists('comparisonPlots/'): + os.mkdir('comparisonPlots/') + + if '100mm' in filename: + ymax = 3 + if 'post' in filename: + title = 'Dimuon mass comparison. ctau = 100mm, postEE' + else: + title = 'Dimuon mass comparison. ctau = 100mm' + elif '10mm' in filename: + ymax = 10 + if 'post' in filename: + title = 'Dimuon mass comparison. ctau = 10mm, postEE' + else: + title = 'Dimuon mass comparison. ctau = 10mm' + elif '1mm' in filename: + ymax = 21 + if 'post' in filename: + title = 'Dimuon mass comparison. ctau = 1mm, postEE' + else: + title = 'Dimuon mass comparison. ctau = 1mm' + + c1 = ROOT.TCanvas('c1', '', 600, 500) + c1.cd() + hdimuon_mass.GetXaxis().SetRangeUser(4,8) + hdimuon_mass.GetYaxis().SetRangeUser(0, ymax) + hdimuon_mass.GetXaxis().SetTitle("Mass (GeV)") + hdimuon_mass.GetYaxis().SetTitle("Events/0.05") + hdimuon_mass.SetLineColor(ROOT.kRed) + hdimuon_mass.SetTitle(title) + #c1.SetLogy() + hdimuon_mass.Draw('HIST') + hdimuon_umass.SetLineColor(ROOT.kBlue) + hdimuon_umass.Draw('HIST, SAME') + legend = ROOT.TLegend(0.15, 0.7, 0.35, 0.8) + legend.AddEntry(hdimuon_mass, "Corrected mass", "l") + legend.AddEntry(hdimuon_umass, "Uncorrected mass", "l") + legend.Draw() + + filename = filename.replace('.root', '') + output_filename = f'comparisonPlots/hdimuon_mass_comparison_{filename}.png' + c1.SaveAs(output_filename) + Image(filename=output_filename) + + +''' diff --git a/python/make_nMinus1efficiencies.py b/python/make_nMinus1efficiencies.py new file mode 100644 index 0000000..21f2230 --- /dev/null +++ b/python/make_nMinus1efficiencies.py @@ -0,0 +1,307 @@ +import ROOT +import os,sys,json +import argparse +from datetime import date +import numpy as np +import copy +import math +sys.path.append('utils') +import plotUtils +from IPython.display import Image +import argparse +import os +import matplotlib.pyplot as plt +import mplhep as hep + +ROOT.gStyle.SetOptStat(0) +parser = argparse.ArgumentParser() +parser.add_argument('directory', type=str) +args = parser.parse_args() +files = [] +directory = args.directory + +def getValues(histo): + values = [] + bins = [] + for n in range(1,histo.GetNbinsX()+1): + values.append(histo.GetBinContent(n)) + bins.append(histo.GetBinLowEdge(n)) + bins.append(histo.GetBinLowEdge(n) + histo.GetBinWidth(n)) + return np.array(values), np.array(bins) + +if not os.path.exists('cuts_plots/'): + os.mkdir('cuts_plots/') + +sum_100mm_allCuts = {} +sum_10mm_allCuts = {} +sum_1mm_allCuts = {} +sum_100mm_noCuts = {} +sum_10mm_noCuts = {} +sum_1mm_noCuts = {} +sum_100mm_Cut1 = {} +sum_10mm_Cut1 = {} +sum_1mm_Cut1 = {} +sum_100mm_Cut2 = {} +sum_10mm_Cut2 = {} +sum_1mm_Cut2 = {} +sum_100mm_Cut3 = {} +sum_10mm_Cut3 = {} +sum_1mm_Cut3 = {} +sum_100mm_Cut4 = {} +sum_10mm_Cut4 = {} +sum_1mm_Cut4 = {} +for subdir in os.listdir(directory): + if 'Jun-14' in subdir: + full_subdir_path = os.path.join(directory, subdir, '') + files = os.listdir(full_subdir_path) + for filename in files: + if '_all' in filename: + continue + if '-2p0_' in filename: + mass = 2 + elif '-7p0_' in filename: + mass = 7 + elif '-10p0_' in filename: + mass = 10 + elif '-20p0_' in filename: + mass = 20 + elif '-30p0_' in filename: + mass = 30 + else: + continue + file = ROOT.TFile.Open(full_subdir_path + filename) + if 'allCuts' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_allCuts = hdimuon_mass.Integral() + current_value_100mm_allCuts = sum_100mm_allCuts.get(mass, 0) + updated_value_100mm_allCuts = current_value_100mm_allCuts + n_hist_100mm_allCuts + sum_100mm_allCuts[mass] = updated_value_100mm_allCuts + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_allCuts = hdimuon_mass.Integral() + current_value_10mm_allCuts = sum_10mm_allCuts.get(mass, 0) + updated_value_10mm_allCuts = current_value_10mm_allCuts + n_hist_10mm_allCuts + sum_10mm_allCuts[mass] = updated_value_10mm_allCuts + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_allCuts = hdimuon_mass.Integral() + current_value_1mm_allCuts = sum_1mm_allCuts.get(mass, 0) + updated_value_1mm_allCuts = current_value_1mm_allCuts + n_hist_1mm_allCuts + sum_1mm_allCuts[mass] = updated_value_1mm_allCuts + + + elif 'noCuts' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_noCuts = hdimuon_mass.Integral() + current_value_100mm_noCuts = sum_100mm_noCuts.get(mass, 0) + updated_value_100mm_noCuts = current_value_100mm_noCuts + n_hist_100mm_noCuts + sum_100mm_noCuts[mass] = updated_value_100mm_noCuts + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_noCuts = hdimuon_mass.Integral() + current_value_10mm_noCuts = sum_10mm_noCuts.get(mass, 0) + updated_value_10mm_noCuts = current_value_10mm_noCuts + n_hist_10mm_noCuts + sum_10mm_noCuts[mass] = updated_value_10mm_noCuts + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_noCuts = hdimuon_mass.Integral() + current_value_1mm_noCuts = sum_1mm_noCuts.get(mass, 0) + updated_value_1mm_noCuts = current_value_1mm_noCuts + n_hist_1mm_noCuts + sum_1mm_noCuts[mass] = updated_value_1mm_noCuts + elif 'noDiMuonAngularSel' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut1 = hdimuon_mass.Integral() + current_value_100mm_Cut1 = sum_100mm_Cut1.get(mass, 0) + updated_value_100mm_Cut1 = current_value_100mm_Cut1 + n_hist_100mm_Cut1 + sum_100mm_Cut1[mass] = updated_value_100mm_Cut1 + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut1 = hdimuon_mass.Integral() + current_value_10mm_Cut1 = sum_10mm_Cut1.get(mass, 0) + updated_value_10mm_Cut1 = current_value_10mm_Cut1 + n_hist_10mm_Cut1 + sum_10mm_Cut1[mass] = updated_value_10mm_Cut1 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut1 = hdimuon_mass.Integral() + current_value_1mm_Cut1 = sum_1mm_Cut1.get(mass, 0) + updated_value_1mm_Cut1 = current_value_1mm_Cut1 + n_hist_1mm_Cut1 + sum_1mm_Cut1[mass] = updated_value_1mm_Cut1 + + elif 'noMaterialVeto' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut2 = hdimuon_mass.Integral() + current_value_100mm_Cut2 = sum_100mm_Cut2.get(mass, 0) + updated_value_100mm_Cut2 = current_value_100mm_Cut2 + n_hist_100mm_Cut2 + sum_100mm_Cut2[mass] = updated_value_100mm_Cut2 + + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut2 = hdimuon_mass.Integral() + current_value_10mm_Cut2 = sum_10mm_Cut2.get(mass, 0) + updated_value_10mm_Cut2 = current_value_10mm_Cut2 + n_hist_10mm_Cut2 + sum_10mm_Cut2[mass] = updated_value_10mm_Cut2 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut2 = hdimuon_mass.Integral() + current_value_1mm_Cut2 = sum_1mm_Cut2.get(mass, 0) + updated_value_1mm_Cut2 = current_value_1mm_Cut2 + n_hist_1mm_Cut2 + sum_1mm_Cut2[mass] = updated_value_1mm_Cut2 + + + elif 'noMuonHitSel' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut3 = hdimuon_mass.Integral() + current_value_100mm_Cut3 = sum_100mm_Cut3.get(mass, 0) + updated_value_100mm_Cut3 = current_value_100mm_Cut3 + n_hist_100mm_Cut3 + sum_100mm_Cut3[mass] = updated_value_100mm_Cut3 + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut3 = hdimuon_mass.Integral() + current_value_10mm_Cut3 = sum_10mm_Cut3.get(mass, 0) + updated_value_10mm_Cut3 = current_value_10mm_Cut3 + n_hist_10mm_Cut3 + sum_10mm_Cut3[mass] = updated_value_10mm_Cut3 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut3 = hdimuon_mass.Integral() + current_value_1mm_Cut3 = sum_1mm_Cut3.get(mass, 0) + updated_value_1mm_Cut3 = current_value_1mm_Cut3 + n_hist_1mm_Cut3 + sum_1mm_Cut3[mass] = updated_value_1mm_Cut3 + + elif 'noMuonIPSel' in subdir: + if '100mm' in filename: + ctau = '100mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_100mm_Cut4 = hdimuon_mass.Integral() + current_value_100mm_Cut4 = sum_100mm_Cut4.get(mass, 0) + updated_value_100mm_Cut4 = current_value_100mm_Cut4 + n_hist_100mm_Cut4 + sum_100mm_Cut4[mass] = updated_value_100mm_Cut4 + elif '10mm' in filename: + ctau = '10mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_10mm_Cut4 = hdimuon_mass.Integral() + current_value_10mm_Cut4 = sum_10mm_Cut4.get(mass, 0) + updated_value_10mm_Cut4 = current_value_10mm_Cut4 + n_hist_10mm_Cut4 + sum_10mm_Cut4[mass] = updated_value_10mm_Cut4 + + elif '1mm' in filename: + ctau = '1mm' + hdimuon_mass = file.Get("hdimuon_mass") + n_hist_1mm_Cut4 = hdimuon_mass.Integral() + current_value_1mm_Cut4 = sum_1mm_Cut4.get(mass, 0) + updated_value_1mm_Cut4 = current_value_1mm_Cut4 + n_hist_1mm_Cut4 + sum_1mm_Cut4[mass] = updated_value_1mm_Cut4 + + +masses = [2,7,10,20,30] +ctaus =['1mm','10mm','100mm'] + +nbins = len(masses)*len(ctaus) +maxx = nbins +minx = 0 + +h_noCuts = ROOT.TH1F('h_noCuts', 'No cuts', nbins, minx, maxx) +h_allCuts = ROOT.TH1F('h_allCuts', 'All cuts', nbins, minx, maxx) +h_Cut1 = ROOT.TH1F('h_Cut1', 'Cut 1', nbins, minx, maxx) +h_Cut2 = ROOT.TH1F('h_Cut2', 'Cut 2', nbins, minx, maxx) +h_Cut3 = ROOT.TH1F('h_Cut3', 'Cut 3', nbins, minx, maxx) +h_Cut4 = ROOT.TH1F('h_Cut4', 'Cut 4', nbins, minx, maxx) +for massindex, mass in enumerate(masses): + h_noCuts.Fill(massindex,sum_1mm_noCuts[mass]/sum_1mm_noCuts[mass]) + h_Cut1.Fill(massindex,sum_1mm_Cut1[mass]/sum_1mm_noCuts[mass]) + h_Cut2.Fill(massindex,sum_1mm_Cut2[mass]/sum_1mm_noCuts[mass]) + h_Cut3.Fill(massindex,sum_1mm_Cut3[mass]/sum_1mm_noCuts[mass]) + h_Cut4.Fill(massindex,sum_1mm_Cut4[mass]/sum_1mm_noCuts[mass]) + h_allCuts.Fill(massindex,sum_1mm_allCuts[mass]/sum_1mm_noCuts[mass]) + + h_noCuts.Fill(massindex + len(masses),sum_10mm_noCuts[mass]/sum_10mm_noCuts[mass]) + h_Cut1.Fill(massindex + len(masses),sum_10mm_Cut1[mass]/sum_10mm_noCuts[mass]) + h_Cut2.Fill(massindex + len(masses),sum_10mm_Cut2[mass]/sum_10mm_noCuts[mass]) + h_Cut3.Fill(massindex + len(masses),sum_10mm_Cut3[mass]/sum_10mm_noCuts[mass]) + h_Cut4.Fill(massindex + len(masses),sum_10mm_Cut4[mass]/sum_10mm_noCuts[mass]) + h_allCuts.Fill(massindex + len(masses),sum_10mm_allCuts[mass]/sum_10mm_noCuts[mass]) + + h_noCuts.Fill(massindex + 2*len(masses),sum_100mm_noCuts[mass]/sum_100mm_noCuts[mass]) + h_Cut1.Fill(massindex + 2*len(masses),sum_100mm_Cut1[mass]/sum_100mm_noCuts[mass]) + h_Cut2.Fill(massindex + 2*len(masses),sum_100mm_Cut2[mass]/sum_100mm_noCuts[mass]) + h_Cut3.Fill(massindex + 2*len(masses),sum_100mm_Cut3[mass]/sum_100mm_noCuts[mass]) + h_Cut4.Fill(massindex + 2*len(masses),sum_100mm_Cut4[mass]/sum_100mm_noCuts[mass]) + h_allCuts.Fill(massindex + 2*len(masses),sum_100mm_allCuts[mass]/sum_100mm_noCuts[mass]) + + +values_noCuts, bins_noCuts = getValues(h_noCuts) +values_allCuts, bins_allCuts = getValues(h_allCuts) +values_Cut1, bins_Cut1 = getValues(h_Cut1) +values_Cut2, bins_Cut2 = getValues(h_Cut2) +values_Cut3, bins_Cut3 = getValues(h_Cut3) +values_Cut4, bins_Cut4 = getValues(h_Cut4) + +hist_array = [values_noCuts, values_Cut1, values_Cut2, values_Cut3, values_Cut4, values_allCuts] +label_array = ['No cuts','No Dimuon Angular Selection','No Material Veto','No Muon Hit Selection','No Muon IP Selection','All cuts'] +yearenergy="(2022, 13.6 TeV)" +bin_labels = [] +j = -1 +masses_plot = masses*3 + +hep.style.use("CMS") +fig, ax = plt.subplots(figsize=(10, 7.5)) +hep.histplot( + hist_array, + stack = False, + bins = bins_allCuts, + histtype = "step", + alpha = 1, + label = label_array, + ax = ax +) + +legend_coords = (0.78,0.4) +ax.set_xlim(0, 15) +expoffset = 0.065 +hep.cms.label(rlabel="") +ax.set_xticks(bins_allCuts) +j = -1 +for i in range(len(masses_plot)): + if i%len(masses) == 0: + j+=1 + x = (bins_allCuts[i] + bins_allCuts[i+1]) / 2 - 0.45 + y = -0.07 + mass_label = f"{masses_plot[i]} GeV" + ctau_label = f"{ctaus[j]}" + plt.text(x, y, mass_label + "\n" + ctau_label, fontsize=10) + +ax.set_xticks(bins_allCuts) +ax.set_xticklabels([''] * len(bins_allCuts)) +ax.legend(loc='center', bbox_to_anchor=legend_coords, fontsize=15) +ax.text(1.01, 1 + expoffset, yearenergy, transform=ax.transAxes, verticalalignment='top', horizontalalignment='right', fontsize=25) +plt.ylabel("N - 1 Efficiency") +plt.savefig("cuts_plots/nMinus1efficiencies_goodplots.png") +plt.close() + diff --git a/utils/histDefinition.py b/utils/histDefinition.py index 778628a..fe35b40 100644 --- a/utils/histDefinition.py +++ b/utils/histDefinition.py @@ -927,6 +927,26 @@ def hist1dDefinition(nbins, low, high, xtitle, ytitle, labels, variable): variable["dimu3danglesv"] = "a3dsv" # + nbins ["dimuumass"] = 15000 + low ["dimuumass"] = 0 + high ["dimuumass"] = 150 + xtitle ["dimuumass"] = "Uncorrected m_{#mu#mu} [GeV]" + ytitle ["dimuumass"] = "Events / 0.01" + variable["dimuumass"] = "umass" + + nbins ["seagullumass"] = 15000 + low ["seagullumass"] = 0 + high ["seagullumass"] = 150 + xtitle ["seagullumass"] = "Uncorrected m_{#mu#mu} [GeV]" + ytitle ["seagullumass"] = "Events / 0.01" + variable["seagullumass"] = "seagullumass" + + nbins ["cowboyumass"] = 15000 + low ["cowboyumass"] = 0 + high ["cowboyumass"] = 150 + xtitle ["cowboyumass"] = "Uncorrected m_{#mu#mu} [GeV]" + ytitle ["cowboyumass"] = "Events / 0.01" + variable["cowboyumass"] = "cowboyumass" nbins ["dimuupdr"] = 100 low ["dimuupdr"] = 0 @@ -3128,6 +3148,14 @@ def histBooking(presel=True, dimuon=True, fourmuon=True, fourmuonosv=True): histname2d["dimuon"].append("hdimuon_genjpsi_lxycomp") histtype2d[histname2d["dimuon"][-1]]="lxyvslxygen" # + histname["dimuon"].append("hdimuon_umass") + histtype[histname["dimuon"][-1]]="dimuumass" + # + histname["dimuon"].append("hdimuon_cowboy") + histtype[histname["dimuon"][-1]]="cowboyumass" + # + histname["dimuon"].append("hdimuon_seagull") + histtype[histname["dimuon"][-1]]="seagullumass" ## # Di-muon from overlapping SV #