From c951e8c65aca18fd5ddcca9b87f321ba688d5599 Mon Sep 17 00:00:00 2001 From: skeshri Date: Thu, 11 Apr 2024 17:27:58 +0200 Subject: [PATCH 1/5] add protection against no-root files to prevent the crash --- helpers.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/helpers.py b/helpers.py index 4880aaa..4a214cc 100644 --- a/helpers.py +++ b/helpers.py @@ -29,7 +29,9 @@ def files_from_path(pathname): ''' FileNameList = [] for file in os.listdir(pathname): - FileNameList.append(os.path.join(pathname, file)) + ext=os.path.splitext(file)[1] + if ext == ".root": + FileNameList.append(os.path.join(pathname, file)) return FileNameList def hadd_anatuple(output_tmp_folder, output_root_file): @@ -70,4 +72,4 @@ def delta_r2(v1, v2): def delta_r(v1, v2): ''' Return Delta r between two objects ''' - return np.sqrt(delta_r2(v1, v2)) \ No newline at end of file + return np.sqrt(delta_r2(v1, v2)) From aedd63ea3c99ed17e15d28f3c122dabb103dc2ee Mon Sep 17 00:00:00 2001 From: skeshri Date: Tue, 16 Apr 2024 08:27:54 +0200 Subject: [PATCH 2/5] added ETau and MuTau class --- ComputeEfficiency/compute_eff.py | 14 ++++ ComputeEfficiency/compute_eff_algo.py | 26 ++++++- ComputeEfficiency/tasks.py | 32 +++++++- ComputeRate/compute_rate.py | 2 +- ComputeRate/compute_rate_quick.py | 13 +++- ComputeRate/tasks.py | 20 ++++- HLTClass/SingleTauDataset.py | 1 + HLTClass/dataset.py | 108 +++++++++++++++++++++++++- config.ini | 38 +++++---- 9 files changed, 229 insertions(+), 25 deletions(-) diff --git a/ComputeEfficiency/compute_eff.py b/ComputeEfficiency/compute_eff.py index 337e317..9668cbe 100644 --- a/ComputeEfficiency/compute_eff.py +++ b/ComputeEfficiency/compute_eff.py @@ -45,6 +45,20 @@ else: print(f"Total number of GenTaus matching with Tau/L1Tau passing {HLT_name} requirements: {N_num}") +if HLT_name == 'HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1': + print(f"Total number of hadronic GenTaus with with vis. pt <= 20 & eta< 2.1 and 1 Muon from Tau: {N_den}") + if PNetMode: + print(f"Total number of GenTaus matching with Jet/L1Tau passing MuTau path with Pnet param {PNetparam}: {N_num}") + else: + print(f"Total number of GenTaus matching with Tau/L1Tau passing {HLT_name} requirements: {N_num}") + +if HLT_name == 'HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1': + print(f"Total number of hadronic GenTaus with with vis. pt <= 20 & eta< 2.1 and 1 Electron from Tau: {N_den}") + if PNetMode: + print(f"Total number of GenTaus matching with Jet/L1Tau passing ETau path with Pnet param {PNetparam}: {N_num}") + else: + print(f"Total number of GenTaus matching with Tau/L1Tau passing {HLT_name} requirements: {N_num}") + print('Computed Eff: ') print(f"Eff : {eff}") print(f"Eff_up : {err_low}") diff --git a/ComputeEfficiency/compute_eff_algo.py b/ComputeEfficiency/compute_eff_algo.py index b26b228..eac1988 100644 --- a/ComputeEfficiency/compute_eff_algo.py +++ b/ComputeEfficiency/compute_eff_algo.py @@ -16,7 +16,7 @@ # only one file otherwise it's too long -FileNameList_eff = f"/afs/cern.ch/work/p/pdebryas/PNetAtHLT/data/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" +FileNameList_eff = f"/afs/cern.ch/work/s/skeshri/TauHLT/Braden/TauHLTOptimzation/PnetAtHLT/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" #FileNameList_eff = files_from_path(f"/afs/cern.ch/work/p/pdebryas/PNetAtHLT/data/EfficiencyDen/{HLT_name}/") if HLT_name == 'HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1': @@ -41,7 +41,29 @@ print(f'Quick EffAlgo computation for {HLT_name}:') N_den, N_num = dataset_eff.ComputeEffAlgo_HLTLooseDeepTauPFTauHPS180_L2NN_eta2p1_v3() +if HLT_name == 'HLT_Ele24_eta2p1_WPTight_Gsf_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1': + from HLTClass.ETauDataset import ETauDataset + dataset_eff = ETauDataset(FileNameList_eff) + + if PNetMode: + print(f'Quick EffAlgo computation for SingleTau path with PNet param {PNetparam}:') + N_den, N_num = dataset_eff.ComputeEffAlgo_ETauPNet(PNetparam) + else: + print(f'Quick EffAlgo computation for {HLT_name}:') + N_den, N_num = dataset_eff.ComputeEffAlgo_ETauDeepNet() + +if HLT_name == 'HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1': + from HLTClass.MuTauDataset import MuTauDataset + dataset_eff = MuTauDataset(FileNameList_eff) + + if PNetMode: + print(f'Quick EffAlgo computation for SingleTau path with PNet param {PNetparam}:') + N_den, N_num = dataset_eff.ComputeEffAlgo_MuTauPNet(PNetparam) + else: + print(f'Quick EffAlgo computation for {HLT_name}:') + N_den, N_num = dataset_eff.ComputeEffAlgo_MuTauDeepNet() + EffAlgo, EffAlgo_low, EffAlgo_up = compute_ratio_witherr(N_num, N_den) print(f"Eff : {EffAlgo}") print(f"Eff_up : {EffAlgo_up}") -print(f"Eff_down : {EffAlgo_low}") \ No newline at end of file +print(f"Eff_down : {EffAlgo_low}") diff --git a/ComputeEfficiency/tasks.py b/ComputeEfficiency/tasks.py index f3a9386..822499b 100644 --- a/ComputeEfficiency/tasks.py +++ b/ComputeEfficiency/tasks.py @@ -38,7 +38,7 @@ def run(self): shutil.rmtree(output_tmp_folder) os.makedirs(output_tmp_folder) - HLT_config = ['HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1', 'HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3', 'HLT_DoubleTauOrSingleTau'] + HLT_config = ['HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1', 'HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3', 'HLT_DoubleTauOrSingleTau','HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1','HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1'] if self.HLT_name not in HLT_config: print(f'HLT name {self.HLT_name} not implemented in the code') raise @@ -64,6 +64,16 @@ def run(self): MC_dataset = DoubleORSingleTauDataset(FileName) MC_dataset.save_Event_Nden_eff_DoubleORSingleTau(output_tmp_file) + if self.HLT_name == 'HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1': + from HLTClass.ETauDataset import ETauDataset + MC_dataset = ETauDataset(FileName) + MC_dataset.save_Event_Nden_eff_ETau(output_tmp_file) + + if self.HLT_name == 'HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1': + from HLTClass.MuTauDataset import MuTauDataset + MC_dataset = MuTauDataset(FileName) + MC_dataset.save_Event_Nden_eff_MuTau(output_tmp_file) + # Hadd the tmp files to a single root file hadd_anatuple(output_tmp_folder, output_root_file) shutil.rmtree(output_tmp_folder) @@ -113,7 +123,7 @@ def run(self): if not os.path.exists(input_root_file): raise('Input root file does not exist') - HLT_config = ['HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1', 'HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3', 'HLT_DoubleTauOrSingleTau'] + HLT_config = ['HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1', 'HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3', 'HLT_DoubleTauOrSingleTau','HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1','HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1'] if self.HLT_name not in HLT_config: print(f'HLT name {self.HLT_name} not implemented in the code') raise @@ -144,3 +154,21 @@ def run(self): MC_dataset.produceRoot_DoubleORSinglePNet(output_root_file, self.PNetparam) else: MC_dataset.produceRoot_DoubleORSingleDeepTau(output_root_file) + + if self.HLT_name == 'HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1': + from HLTClass.ETauDataset import ETauDataset + + MC_dataset = ETauDataset(input_root_file) + if self.PNetMode: + MC_dataset.produceRoot_ETauPNet(output_root_file, self.PNetparam) + else: + MC_dataset.produceRoot_ETauDeepNet(output_root_file) + + if self.HLT_name == 'HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1': + from HLTClass.MuTauDataset import MuTauDataset + + MC_dataset = MuTauDataset(input_root_file) + if self.PNetMode: + MC_dataset.produceRoot_MuTauPNet(output_root_file, self.PNetparam) + else: + MC_dataset.produceRoot_MuTauDeepNet(output_root_file) diff --git a/ComputeRate/compute_rate.py b/ComputeRate/compute_rate.py index f700c32..6629b7d 100644 --- a/ComputeRate/compute_rate.py +++ b/ComputeRate/compute_rate.py @@ -46,4 +46,4 @@ print('Computed rate: ') print(f"rate : {eff*L1A_physics}") print(f"rate_up : {err_low*L1A_physics}") -print(f"rate_down : {err_up*L1A_physics}") \ No newline at end of file +print(f"rate_down : {err_up*L1A_physics}") diff --git a/ComputeRate/compute_rate_quick.py b/ComputeRate/compute_rate_quick.py index 493c281..6461524 100644 --- a/ComputeRate/compute_rate_quick.py +++ b/ComputeRate/compute_rate_quick.py @@ -53,8 +53,19 @@ else: N_den, N_num = dataset_rate.get_Nnum_Nden_HLT_DoubleORSingleDeepTau() +if HLT_name == 'HLT_Ele24_eta2p1_WPTight_Gsf_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1': + from HLTClass.ETauDataset import ETauDataset + + dataset_rate = ETauDataset(FileNameList_rate) + + if PNetMode: + N_den, N_num = dataset_rate.get_Nnum_Nden_ETauPNet(PNetparam) + else: + N_den, N_num = dataset_rate.get_Nnum_Nden_ETauDeepNet() + + rate, rate_low, rate_up = compute_ratio_witherr(N_num, N_den) print(f"Rate : {rate*L1A_physics}") print(f"Rate_up : {rate_up*L1A_physics}") -print(f"Rate_down : {rate_low*L1A_physics}") \ No newline at end of file +print(f"Rate_down : {rate_low*L1A_physics}") diff --git a/ComputeRate/tasks.py b/ComputeRate/tasks.py index 7f16021..d57b892 100644 --- a/ComputeRate/tasks.py +++ b/ComputeRate/tasks.py @@ -107,7 +107,7 @@ def run(self): event_counter[intput_root_file] = {} print(f"For {os.path.basename(intput_root_file)}:") - HLT_config = ['HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1', 'HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3', 'HLT_DoubleTauOrSingleTau'] + HLT_config = ['HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1', 'HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3', 'HLT_DoubleTauOrSingleTau','HLT_Ele24_eta2p1_WPTight_Gsf_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1','HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1'] if self.HLT_name not in HLT_config: print(f'HLT name {self.HLT_name} not implemented in the code') raise @@ -139,6 +139,24 @@ def run(self): else: N_den_i, N_num_i = eph_dataset.get_Nnum_Nden_HLT_DoubleORSingleDeepTau() + if self.HLT_name == 'HLT_Ele24_eta2p1_WPTight_Gsf_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1': + from HLTClass.ETauDataset import ETauDataset + + eph_dataset = ETauDataset(intput_root_file) + if self.PNetMode: + N_den_i, N_num_i = eph_dataset.get_Nnum_Nden_ETauPNet(self.PNetparam) + else: + N_den_i, N_num_i = eph_dataset.get_Nnum_Nden_HLT_ETauDeepNet() + + if self.HLT_name == 'HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1': + from HLTClass.MuTauDataset import MuTauDataset + + eph_dataset = MuTauDataset(intput_root_file) + if self.PNetMode: + N_den_i, N_num_i = eph_dataset.get_Nnum_Nden_MuTauPNet(self.PNetparam) + else: + N_den_i, N_num_i = eph_dataset.get_Nnum_Nden_HLT_MuTauDeepNet() + event_counter[intput_root_file]['N_den'] = N_den_i event_counter[intput_root_file]['N_num'] = N_num_i diff --git a/HLTClass/SingleTauDataset.py b/HLTClass/SingleTauDataset.py index c9b27d3..5434f5e 100644 --- a/HLTClass/SingleTauDataset.py +++ b/HLTClass/SingleTauDataset.py @@ -95,6 +95,7 @@ def evt_sel_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3(events, is_gen = False): SingleTau_mask = Tau_selection_SingleTau(events) # at least 1 L1tau/ recoTau should pass SingleTau_evt_mask = (ak.sum(L1Tau_Tau120er2p1L2NN_mask, axis=-1) >= 1) & (ak.sum(SingleTau_mask, axis=-1) >= 1) + print("SingleTau_evt_mask",SingleTau_evt_mask) if is_gen: # if MC data, at least 1 GenTau should also pass GenTau_mask = hGenTau_selection(events) diff --git a/HLTClass/dataset.py b/HLTClass/dataset.py index dce5991..9988a74 100644 --- a/HLTClass/dataset.py +++ b/HLTClass/dataset.py @@ -5,11 +5,24 @@ import awkward as ak from helpers import delta_r +### Get L1Taus def get_L1Taus(events): L1taus_dict = {"pt": events["L1Tau_pt"].compute(), "eta": events["L1Tau_eta"].compute(), "phi": events["L1Tau_phi"].compute(), "Iso": events["L1Tau_hwIso"].compute()} L1Taus = ak.zip(L1taus_dict) return L1Taus +### Get L1Egamma +def get_L1Egamma(events): + L1electrons_dict = {"pt": events["L1Egamma_pt"].compute(), "eta": events["L1Egamma_eta"].compute(), "phi": events["L1Egamma_phi"].compute(), "Iso": events["L1Egamma_hwIso"].compute()} + L1Electrons = ak.zip(L1electrons_dict) + return L1Electrons + +### Get L1Muons +def get_L1Muons(events): + L1muons_dict = {"pt": events["L1Muon_pt"].compute(), "eta": events["L1Muon_eta"].compute(), "phi": events["L1Muon_phi"].compute()} + L1Muons = ak.zip(L1muons_dict) + return L1Muons + def get_Taus(events): taus_dict = {"pt": events["Tau_pt"].compute(), "eta": events["Tau_eta"].compute(), "phi": events["Tau_phi"].compute(), "deepTauVSjet": events["Tau_deepTauVSjet"].compute()} Taus = ak.zip(taus_dict) @@ -20,6 +33,17 @@ def get_Jets(events): Jets = ak.zip(jets_dict) return Jets +def get_Muons(events): + muons_dict = {"pt": events["Muon_pt"].compute(), "eta": events["Muon_eta"].compute(), "phi": events["Muon_phi"].compute(), "muon_passSingleMuon": events['Muon_passSingleMuon'].compute(), "muon_passMuTau": events['Muon_passMuTau'].compute()} + Muons = ak.zip(muons_dict) + return Muons + +def get_Electrons(events): + electrons_dict = {"pt": events["Electron_pt"].compute(), "eta": events["Electron_eta"].compute(), "phi": events["Electron_phi"].compute(), "electron_passSingleElectron": events['Electron_passSingleElectron'].compute(), "electron_passETau": events['Electron_passETau'].compute()} + Electrons = ak.zip(electrons_dict) + return Electrons + + def get_GenTaus(events): gentaus_dict = {"pt": events["GenLepton_pt"].compute(), "eta": events["GenLepton_eta"].compute(), "phi": events["GenLepton_phi"].compute(), "nChargedHad": events["GenLepton_nChargedHad"].compute(), "nNeutralHad": events["GenLepton_nNeutralHad"].compute(), "DecayMode": events["GenLepton_DecayMode"].compute(), "charge": events["GenLepton_charge"].compute()} GenTaus = ak.zip(gentaus_dict) @@ -30,6 +54,28 @@ def hGenTau_selection(events): hGenTau_mask = (events['GenLepton_pt'].compute() >= 20) & (np.abs(events['GenLepton_eta'].compute()) <= 2.3) & (events['GenLepton_kind'].compute() == 5) return hGenTau_mask +def get_GenMuons(events): + genmuons_dict = {"pt": events["GenLepton_pt"].compute(), "eta": events["GenLepton_eta"].compute(), "phi": events["GenLepton_phi"].compute(), "nChargedHad": events["GenLepton_nChargedHad"].compute(), "nNeutralHad": events["GenLepton_nNeutralHad"].compute(), "DecayMode": events["GenLepton_DecayMode"].compute(), "charge": events["GenLepton_charge"].compute()} + GenMuons = ak.zip(genmuons_dict) + return GenMuons + +def get_GenElectrons(events): + genelectrons_dict = {"pt": events["GenLepton_pt"].compute(), "eta": events["GenLepton_eta"].compute(), "phi": events["GenLepton_phi"].compute(), "nChargedHad": events["GenLepton_nChargedHad"].compute(), "nNeutralHad": events["GenLepton_nNeutralHad"].compute(), "DecayMode": events["GenLepton_DecayMode"].compute(), "charge": events["GenLepton_charge"].compute()} + GenElectrons = ak.zip(genelectrons_dict) + return GenElectrons + + +def GenElectron_selection(events): + # return mask for GenLepton for Tau decay into electron passing minimal selection + GenElectron_mask = (events['GenLepton_pt'].compute() >= 22 ) & (np.abs(events['GenLepton_eta'].compute()) <= 2.1) & (events['GenLepton_kind'].compute() == 3) + return GenElectron_mask + +def GenMuon_selection(events): + # return mask for GenLepton for Tau decay into muon passing minimal selection + GenMuon_mask = (events['GenLepton_pt'].compute() >= 18) & (np.abs(events['GenLepton_eta'].compute()) <= 2.1) & (events['GenLepton_kind'].compute() == 4) + return GenMuon_mask + + def matching_L1Taus_obj(L1Taus, Obj, dR_matching_min = 0.5): obj_inpair, l1taus_inpair = ak.unzip(ak.cartesian([Obj, L1Taus], nested=True)) dR_obj_l1taus = delta_r(obj_inpair, l1taus_inpair) @@ -38,6 +84,22 @@ def matching_L1Taus_obj(L1Taus, Obj, dR_matching_min = 0.5): mask = ak.any(mask_obj_l1taus, axis=-1) return mask +def matching_L1Muons_obj(L1Muons, Obj, dR_matching_min = 0.5): + obj_inpair, l1muons_inpair = ak.unzip(ak.cartesian([Obj, L1Muons], nested=True)) + dR_obj_l1muons = delta_r(obj_inpair, l1muons_inpair) + mask_obj_l1muons = (dR_obj_l1muons < dR_matching_min) + + mask = ak.any(mask_obj_l1muons, axis=-1) + return mask + +def matching_L1Egamma_obj(L1Egamma, Obj, dR_matching_min = 0.5): + obj_inpair, l1egamma_inpair = ak.unzip(ak.cartesian([Obj, L1Egamma], nested=True)) + dR_obj_l1egamma = delta_r(obj_inpair, l1egamma_inpair) + mask_obj_l1egamma = (dR_obj_l1egamma < dR_matching_min) + + mask = ak.any(mask_obj_l1egamma, axis=-1) + return mask + def matching_Gentaus(L1Taus, Taus, GenTaus, dR_matching_min = 0.5): gentaus_inpair, l1taus_inpair = ak.unzip(ak.cartesian([GenTaus, L1Taus], nested=True)) dR_gentaus_l1taus = delta_r(gentaus_inpair, l1taus_inpair) @@ -50,6 +112,31 @@ def matching_Gentaus(L1Taus, Taus, GenTaus, dR_matching_min = 0.5): matching_mask = ak.any(mask_gentaus_l1taus, axis=-1) & ak.any(mask_gentaus_taus, axis=-1) # Gentau should match l1Taus and Taus return matching_mask +def matching_GenMuons(L1Muons, Muons, GenMuons, dR_matching_min = 0.5): + genmuons_inpair, l1muons_inpair = ak.unzip(ak.cartesian([GenMuons, L1Muons], nested=True)) + dR_genmuons_l1muons = delta_r(genmuons_inpair, l1muons_inpair) + mask_genmuons_l1muons = (dR_genmuons_l1muons < dR_matching_min) + + genmuons_inpair, muons_inpair = ak.unzip(ak.cartesian([GenMuons, Muons], nested=True)) + dR_genmuons_muons = delta_r(genmuons_inpair, muons_inpair) + mask_genmuons_muons = (dR_genmuons_muons < dR_matching_min) + + matching_mask = ak.any(mask_genmuons_l1muons, axis=-1) & ak.any(mask_genmuons_muons, axis=-1) # Genmuon should match l1muons and Muons + return matching_mask + +def matching_GenElectrons(L1Egamma, Electrons, GenElectrons, dR_matching_min = 0.5): + genelectrons_inpair, l1egamma_inpair = ak.unzip(ak.cartesian([GenElectrons, L1Egamma], nested=True)) + dR_genelectrons_l1egamma = delta_r(genelectrons_inpair, l1egamma_inpair) + mask_genelectrons_l1egamma = (dR_genelectrons_l1egamma < dR_matching_min) + + genelectrons_inpair, electrons_inpair = ak.unzip(ak.cartesian([GenElectrons, Electrons], nested=True)) + dR_genelectrons_electrons = delta_r(genelectrons_inpair, electrons_inpair) + mask_genelectrons_electrons = (dR_genelectrons_electrons < dR_matching_min) + + matching_mask = ak.any(mask_genelectrons_l1egamma, axis=-1) & ak.any(mask_genelectrons_electrons, axis=-1) # Genelectrons should match l1egamma and Electrons + return matching_mask + + def compute_PNet_charge_prob(probTauP, probTauM): return np.abs(np.ones(len(probTauP))*0.5 - probTauP/(probTauP + probTauM)) @@ -61,6 +148,8 @@ def iterable(arg): + + class Dataset: def __init__(self, fileName): self.fileName = fileName @@ -186,7 +275,24 @@ def Save_Event_Nden_Eff(self, events, GenLepton, evt_mask, tmp_file): 'Jet_PNet_ptcorr', 'Jet_pt', 'Jet_eta', - 'Jet_phi'] + 'Jet_phi', + 'L1Egamma_pt', + 'L1Egamma_eta', + 'L1Egamma_phi', + 'L1Egamma_hwIso', + 'Electron_pt', + 'Electron_eta', + 'Electron_phi', + 'Electron_passSingleElectron', + 'Electron_passETau', + 'L1Muon_pt', + 'L1Muon_eta', + 'L1Muon_phi', + 'Muon_pt', + 'Muon_phi', + 'Muon_eta', + 'Muon_passSingleMuon', + 'Muon_passMuTau'] saved_info_GenLepton = ['pt', 'eta', diff --git a/config.ini b/config.ini index 94f048c..ec3644e 100644 --- a/config.ini +++ b/config.ini @@ -2,44 +2,48 @@ ref_run = 362617 LumiSectionsRange_low = 0 LumiSectionsRange_up = 245 -Area = 2022G +Area = 2023D # for run ref_run and lumisectionRange, OMS values for the L1 rate (L1A physics) L1A_physics = 91374.04 [HLT] # HLT name to study -HLTname = HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1 +#HLTname = HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1 +#HLTname = HLT_Ele24_eta2p1_WPTight_Gsf_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1 +HLTname = HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1 +#HLTname = HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3 # rate in CMS OMS for HLTname (for run ref_run and within lumisectionRange) HLT_rate = 50.12 [DATA] # path to the samples ... -SamplesPath = /eos/cms/store/group/phys_tau/TauTrigger/Run3_HLT/prod_v3/ +#SamplesPath = /eos/cms/store/group/phys_tau/TauTrigger/Run3_HLT/prod_v3/ +SamplesPath = /eos/cms/store/group/phys_tau/TauTrigger/Run3_HLT/prod_2024_v1 +#SamplesPath = /afs/cern.ch/user/s/skeshri + # .. For rate computation -number_of_ephemeral_folder = 9 +number_of_ephemeral_folder = 2 # .. For eff computation MCDataFolderNames = - ZprimeToTauTau_M-4000 - VBFHToTauTau_M125 - GluGluHToTauTau_M-125 - GluGluHToTauTau_M-125_ext1 - GluGlutoHHto2B2Tau_kl-1p00_kt-1p00_c2-0p00 - VBFHHto2B2Tau_CV-1_C2V-1_C3-1 +# GluGluHToTauTau_M-125 + GluGlutoHHto2B2Tau_kl-1p00_kt-1p00_c2-3p00 # path where to store the results, another tag would be add depending on HLTname so no need to specify it here -result_rate = /afs/cern.ch/user/p/pdebryas/PnetAtHLT/PnetAtHLT/ComputeRate/result/ -result_eff = /afs/cern.ch/user/p/pdebryas/PnetAtHLT/PnetAtHLT/ComputeEfficiency/result/ +#result_rate = /afs/cern.ch/user/p/pdebryas/PnetAtHLT/PnetAtHLT/ComputeRate/result/ +#result_eff = /afs/cern.ch/user/p/pdebryas/PnetAtHLT/PnetAtHLT/ComputeEfficiency/result/ +result_rate = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/result/ +result_eff = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/result/ # path where to store anatuples (events which pass denominator selection in eff/rate computation) -EffDenPath = /afs/cern.ch/work/p/pdebryas/PNetAtHLT/data/EfficiencyDen/ -RateDenPath = /afs/cern.ch/work/p/pdebryas/PNetAtHLT/data/RateDen/ +EffDenPath = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/EfficiencyDen/ +RateDenPath = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/RateDen/ # path where to store tmp files -tmpPath = /afs/cern.ch/work/p/pdebryas/PNetAtHLT/tmp/ +tmpPath = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/PnetAtHLT/tmp/ [OPT] # Use of current DeepTau WP: PNet_mode = false; else PNet_mode = true -PNet_mode = false +PNet_mode = true # Pnet_tauhm+Pnet_tauhp > PNet_WP(PNet_t1, PNet_t2) PNet_t1 = 0.56 # Pnet_tauhm+Pnet_tauhp > PNet_WP(PNet_t1, PNet_t2) PNet_t2 = 0.47 # Pnet_chargeprob > PNet_t3 -PNet_t3 = 30 \ No newline at end of file +PNet_t3 = 30 From 7f1c945a5b264b4da8bf62984eff6c066c367999 Mon Sep 17 00:00:00 2001 From: skeshri Date: Mon, 22 Apr 2024 15:59:56 +0200 Subject: [PATCH 3/5] Update in MuTau and ETau scripts for the optimisation --- ComputeEfficiency/compute_eff_algo.py | 3 +- HLTClass/ETauDataset.py | 504 ++++++++++++++++++++++++++ HLTClass/MuTauDataset.py | 455 +++++++++++++++++++++++ HLTClass/dataset.py | 22 +- Optimisation/helpers_opt.py | 2 +- Optimisation/merger.py | 25 ++ Optimisation/plot_manualDiTau.py | 106 ++++++ Optimisation/run_optETau.py | 61 ++++ Optimisation/tasks.py | 397 ++++++++++++++++++++ 9 files changed, 1572 insertions(+), 3 deletions(-) create mode 100644 HLTClass/ETauDataset.py create mode 100644 HLTClass/MuTauDataset.py create mode 100644 Optimisation/merger.py create mode 100644 Optimisation/plot_manualDiTau.py create mode 100644 Optimisation/run_optETau.py create mode 100644 Optimisation/tasks.py diff --git a/ComputeEfficiency/compute_eff_algo.py b/ComputeEfficiency/compute_eff_algo.py index eac1988..af36e7f 100644 --- a/ComputeEfficiency/compute_eff_algo.py +++ b/ComputeEfficiency/compute_eff_algo.py @@ -16,7 +16,8 @@ # only one file otherwise it's too long -FileNameList_eff = f"/afs/cern.ch/work/s/skeshri/TauHLT/Braden/TauHLTOptimzation/PnetAtHLT/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" +#FileNameList_eff = f"/afs/cern.ch/work/s/skeshri/TauHLT/Braden/TauHLTOptimzation/PnetAtHLT/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" +FileNameList_eff = f"/afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/EfficiencyDen/{HLT_name}/GluGlutoHHto2B2Tau_kl-1p00_kt-1p00_c2-3p00.root" #FileNameList_eff = files_from_path(f"/afs/cern.ch/work/p/pdebryas/PNetAtHLT/data/EfficiencyDen/{HLT_name}/") if HLT_name == 'HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1': diff --git a/HLTClass/ETauDataset.py b/HLTClass/ETauDataset.py new file mode 100644 index 0000000..e38eb13 --- /dev/null +++ b/HLTClass/ETauDataset.py @@ -0,0 +1,504 @@ +import awkward as ak +import numpy as np +from HLTClass.dataset import Dataset +from HLTClass.dataset import get_L1Taus, get_Taus, get_Jets, get_GenTaus, hGenTau_selection, matching_Gentaus, matching_L1Taus_obj, compute_PNet_charge_prob, get_L1Egamma, get_Electrons,get_GenElectrons, GenElectron_selection, matching_GenElectrons, matching_L1Egamma_obj +from helpers import delta_r + + +# ------------------------------ functions for ETau with PNet --------------------------------------------------------------- + +def compute_PNet_WP_ETau(tau_pt, par): + # return PNet WP for ETau + t1 = par[0] + t2 = par[1] + t3 = 0.001 + t4 = 0 + x1 = 130 + x2 = 200 + x3 = 500 + x4 = 1000 + PNet_WP = tau_pt*0. + ones = tau_pt/tau_pt + PNet_WP = ak.where((tau_pt <= ones*x1) == False, PNet_WP, ones*t1) + PNet_WP = ak.where((tau_pt >= ones*x4) == False, PNet_WP, ones*t4) + PNet_WP = ak.where(((tau_pt < ones*x2) & (tau_pt > ones*x1)) == False, PNet_WP, (t2 - t1) / (x2 - x1) * (tau_pt - ones*x1) + ones*t1) + PNet_WP = ak.where(((tau_pt >= ones*x2) & (tau_pt < ones*x3))== False, PNet_WP, (t3 - t2) / (x3 - x2) * (tau_pt - ones*x2) + ones*t2) + PNet_WP = ak.where(((tau_pt >= ones*x3) & (tau_pt < ones*x4))== False, PNet_WP, (t4 - t3) / (x4 - x3) * (tau_pt - ones*x3) + ones*t3) + return PNet_WP + + +def Jet_selection_Tau(events, par, apply_PNET_WP = True): + # return mask for Jet passing selection for ETau path + Jet_pt_corr = events['Jet_pt'].compute()*events['Jet_PNet_ptcorr'].compute() + Jets_mask = (events['Jet_pt'].compute() >= 30) & (np.abs(events['Jet_eta'].compute()) <= 2.3) & (Jet_pt_corr >= 30) + if apply_PNET_WP: + probTauP = events['Jet_PNet_probtauhp'].compute() + probTauM = events['Jet_PNet_probtauhm'].compute() + Jets_mask = Jets_mask & ((probTauP + probTauM) >= compute_PNet_WP_ETau(Jet_pt_corr, par)) & (compute_PNet_charge_prob(probTauP, probTauM) >= 0) + return Jets_mask + +def Electron_selection(events): + Electron_mask = (events['Electron_pt'].compute() >= 20) & (np.abs(events['Electron_eta'].compute()) <= 2.3) & (events['Electron_passSingleElectron'].compute() | events['Electron_passETau'].compute()) + return Electron_mask + + +def matching_dR_Min0p3(L1Taus, Obj, dR_matching_min = 0.3): + obj_inpair, l1taus_inpair = ak.unzip(ak.cartesian([Obj, L1Taus], nested=True)) + dR_obj_l1taus = delta_r(obj_inpair, l1taus_inpair) + mask_obj_l1taus = (dR_obj_l1taus > dR_matching_min) + + mask = ak.any(mask_obj_l1taus, axis=-1) + return mask + + + +def evt_sel_ETau(events, par, is_gen = False): + # Selection of event passing condition of ETau with PNet HLT path + mask of objects passing those conditions + + L1_EG_mask = L1_LooseIsoEG22er2p1_selection(events) | L1_LooseIsoEG24er2p1_selcetion(events) + L1_Tau_mask = L1_IsoTau26er2p1_selection(events) | L1_Tau70er2p1_selection(events) | L1_IsoTau27er2p1_selcetion(events) + + Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = True) + + Electron_mask = Electron_selection(events) + ETau_evt_mask = (ak.sum(L1_EG_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Electron_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Egamma = get_L1Egamma(events) + L1Egamma_EG = L1Egamma[L1_EG_mask] + + Jets = get_Jets(events) + Jets_Tau = Jets[Tau_mask] + + Electrons = get_Electrons(events) + Selected_Electrons = Electrons[Electron_mask] + + matching_dR_Min0p3_L1 = matching_dR_Min0p3(L1Taus_Tau, L1Egamma_EG, dR_matching_min=0.3) + + ETau_evt_mask = ETau_evt_mask & matching_dR_Min0p3_L1 + + if is_gen: + # if MC data, at least 1 GenTau should also pass + GenTau_mask = hGenTau_selection(events) + GenElectron_mask = GenElectron_selection(events) + ETau_evt_mask = ETau_evt_mask & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenElectron_mask, axis=-1) >= 1) + + # matching + if is_gen: + GenTaus = get_GenTaus(events) + GenTaus_Tau = GenTaus[GenTau_mask] + matchingGentaus_mask = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Tau) + # at least 1 GenTau should match L1Tau/Taus + evt_Taumask_matching = (ak.sum(matchingGentaus_mask, axis=-1) >= 1) + + GenElectrons = get_GenElectrons(events) + GenElectrons_Electron = GenElectrons[GenElectron_mask] + matchingGenElectrons_mask = matching_GenElectrons(L1Egamma_EG, Selected_Electrons, GenElectrons_Electron) + # at least 1 GenTau should match L1Tau/Taus + evt_Electronmask_matching = (ak.sum(matchingGenElectrons_mask, axis=-1) >= 1) + + else: + matchingJets_mask = matching_L1Taus_obj(L1Taus_Tau, Jets_Tau) + # at least 1 Tau should match L1Tau + evt_Taumask_matching = (ak.sum(matchingJets_mask, axis=-1) >= 1) + + matchingElectrons_mask = matching_L1Egamma_obj(L1Egamma_EG, Selected_Electrons) + evt_Electronmask_matching = (ak.sum(matchingElectrons_mask, axis=-1) >= 1) + + ETau_evt_mask = ETau_evt_mask & evt_Taumask_matching & evt_Electronmask_matching + if is_gen: + return ETau_evt_mask, matchingGentaus_mask, matchingGenElectrons_mask + else: + return ETau_evt_mask, matchingJets_mask, matchingElectrons_mask + +# ------------------------------ functions for HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1 --------------------------------------------------- +def compute_DeepTau_WP_ETau(tau_pt): + # return DeepTau WP for LooseDeepTauPFTauHPS30_eta2p1 + t1 = 0.649 + t2 = 0.441 + t3 = 0.05 + x1 = 35 + x2 = 100 + x3 = 300 + Tau_WP = tau_pt*0. + ones = tau_pt/tau_pt + Tau_WP = ak.where((tau_pt <= ones*x1) == False, Tau_WP, ones*t1) + Tau_WP = ak.where((tau_pt >= ones*x3) == False, Tau_WP, ones*t3) + Tau_WP = ak.where(((tau_pt < ones*x2) & (tau_pt > ones*x1)) == False, Tau_WP, (t2 - t1) / (x2 - x1) * (tau_pt - ones*x1) + ones*t1) + Tau_WP = ak.where(((tau_pt >= ones*x2) & (tau_pt < ones*x3))== False, Tau_WP, (t3 - t2) / (x3 - x2) * (tau_pt - ones*x2) + ones*t2) + return Tau_WP + +def Tau_selection_Tau(events, apply_DeepTau_WP = True): + # return mask for Tau passing selection for HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1 + tau_pt = events['Tau_pt'].compute() + Tau_mask = (tau_pt >= 30) & (np.abs(events['Tau_eta'].compute()) <= 2.1) + if apply_DeepTau_WP: + Tau_mask = Tau_mask & (events['Tau_deepTauVSjet'].compute() >= compute_DeepTau_WP_ETau(tau_pt)) + return Tau_mask + +def evt_sel_HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1(events, is_gen = False): + # Selection of event passing condition of HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1 + mask of objects passing those conditions + + L1_EG_mask = L1_LooseIsoEG22er2p1_selection(events) | L1_LooseIsoEG24er2p1_selcetion(events) + L1_Tau_mask = L1_IsoTau26er2p1_selection(events) | L1_Tau70er2p1_selection(events) | L1_IsoTau27er2p1_selcetion(events) + + Tau_mask = Tau_selection_Tau(events) + + Electron_mask = Electron_selection(events) + # at least 1 L1tau/ recoJet should pass + ETau_evt_mask = (ak.sum(L1_EG_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Electron_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Egamma = get_L1Egamma(events) + L1Egamma_EG = L1Egamma[L1_EG_mask] + + Taus = get_Taus(events) + Taus_Tau = Taus[Tau_mask] + + Electrons = get_Electrons(events) + Selected_Electrons = Electrons[Electron_mask] + + matching_dR_Min0p3_L1 = matching_dR_Min0p3(L1Taus_Tau, L1Egamma_EG, dR_matching_min=0.3) + + ETau_evt_mask = ETau_evt_mask & matching_dR_Min0p3_L1 + print("ETau_evt_mask",ETau_evt_mask) + + if is_gen: + # if MC data, at least 1 GenTau should also pass + GenTau_mask = hGenTau_selection(events) + GenElectron_mask = GenElectron_selection(events) + ETau_evt_mask = ETau_evt_mask & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenElectron_mask, axis=-1) >= 1) + + # matching + if is_gen: + GenTaus = get_GenTaus(events) + GenTaus_Tau = GenTaus[GenTau_mask] + matchingGentaus_mask = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Tau) + # at least 1 GenTau should match L1Tau/Taus + evt_Taumask_matching = (ak.sum(matchingGentaus_mask, axis=-1) >= 1) + + GenElectrons = get_GenElectrons(events) + GenElectrons_Electron = GenElectrons[GenElectron_mask] + matchingGenElectrons_mask = matching_GenElectrons(L1Egamma_EG, Selected_Electrons, GenElectrons_Electron) + # at least 1 GenTau decay to electron should match L1Tau/Taus + evt_Electronmask_matching = (ak.sum(matchingGenElectrons_mask, axis=-1) >= 1) + + else: + matchingTaus_mask = matching_L1Taus_obj(L1Taus_Tau, Taus_Tau) + # at least 1 Tau should match L1Tau + evt_Taumask_matching = (ak.sum(matchingTaus_mask, axis=-1) >= 1) + + matchingElectrons_mask = matching_L1Egamma_obj(L1Egamma_EG, Selected_Electrons) + # at least 1 L1EG should match L1Tau + evt_Electronmask_matching = (ak.sum(matchingElectrons_mask, axis=-1) >= 1) + + ETau_evt_mask = ETau_evt_mask & evt_Taumask_matching & evt_Electronmask_matching + if is_gen: + return ETau_evt_mask, matchingGentaus_mask, matchingGenElectrons_mask + else: + return ETau_evt_mask, matchingTaus_mask, matchingElectrons_mask + +# ------------------------------ Common functions for Ditau path --------------------------------------------------------------- +def L1Tau_Tau130er2p1_selection(events): + L1Taus = get_L1Taus(events) + # return mask for L1tau passing Tau120er2p1 selection + L1_Tau120er2p1_mask = (events['L1Tau_pt'].compute() >= 130) & (events['L1Tau_eta'].compute() <= 2.131) & (events['L1Tau_eta'].compute() >= -2.131) + return L1_Tau130er2p1_mask + +def L1_LooseIsoEG22er2p1_selection(events): + L1_LooseIsoEG22er2p1 = (events['L1Egamma_pt'].compute() >= 22) & (events['L1Egamma_eta'].compute() <= 2.131) & (events['L1Egamma_eta'].compute() >= -2.131) + return L1_LooseIsoEG22er2p1 + +def L1_IsoTau26er2p1_selection(events): + L1_IsoTau26er2p1 = (events['L1Tau_pt'].compute() >= 26) & (events['L1Tau_eta'].compute() <= 2.131) & (events['L1Tau_eta'].compute() >= -2.131) & (events['L1Tau_hwIso'].compute() > 0) + return L1_IsoTau26er2p1 + +def L1_Tau70er2p1_selection(events): + L1_Tau70er2p1 = (events['L1Tau_pt'].compute() >= 70) & (events['L1Tau_eta'].compute() <= 2.131) & (events['L1Tau_eta'].compute() >= -2.131) + return L1_Tau70er2p1 + +def L1_LooseIsoEG24er2p1_selcetion(events): + L1_LooseIsoEG24er2p1 = (events['L1Egamma_pt'].compute() >= 24) & (events['L1Egamma_eta'].compute() <= 2.131) & (events['L1Egamma_eta'].compute() >= -2.131) + return L1_LooseIsoEG24er2p1 + +def L1_IsoTau27er2p1_selcetion(events): + L1_IsoTau27er2p1 = (events['L1Tau_pt'].compute() >= 27) & (events['L1Tau_eta'].compute() <= 2.131) & (events['L1Tau_eta'].compute() >= -2.131) & (events['L1Tau_hwIso'].compute() > 0) + return L1_IsoTau27er2p1 + + +def Denominator_Selection_ETau(GenLepton): + # return mask for event passing minimal Gen requirements for ETau HLT (1 hadronic Taus with min vis. pt and eta and 1 Tau decayed to electron) + mask_GenTau = (GenLepton['kind'] == 5) + mask_GenEle = (GenLepton['kind'] == 3) + ev_mask = (ak.sum(mask_GenTau, axis=-1) == 1) & (ak.sum(mask_GenEle, axis=-1) == 1) # 1 Gen tau and one Gen Ele should pass this requirements + return ev_mask + +class ETauDataset(Dataset): + def __init__(self, fileName): + Dataset.__init__(self, fileName) + + # ------------------------------ functions to Compute Rate --------------------------------------------------------------- + + def get_Nnum_Nden_ETauDeepNet(self): + print(f'Computing rate for HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1:') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + + ETau_evt_mask, matchingTaus_mask, matchingElectrons_mask = evt_sel_HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1(events, is_gen = False) + N_num = len(events[ETau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def get_Nnum_Nden_ETauPNet(self, par): + print(f'Computing Rate for ETau path with param: {par}') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + ETau_evt_mask, matchingJets_mask, matchingElectrons_mask = evt_sel_ETau(events, par, is_gen = False) + N_num = len(events[ETau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def get_Nnum_Nden_ETauDeepNet_nPV(self, nPV): + print(f'Computing rate for HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1 for nPV = {nPV}:') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + events = events[events['nPFPrimaryVertex'].compute() == nPV] + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + + ETau_evt_mask, matchingTaus_mask, matchingElectrons_mask = evt_sel_HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1(events, is_gen = False) + N_num = len(events[ETau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def get_Nnum_Nden_ETauPNet_nPV(self, par, nPV): + print(f'Computing Rate for ETau path with param: {par} and for nPV = {nPV}:') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + events = events[events['nPFPrimaryVertex'].compute() == nPV] + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + + ETau_evt_mask, matchingJets_mask, matchingElectrons_mask = evt_sel_ETau(events, par, is_gen = False) + N_num = len(events[ETau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + +# ------------------------------ functions for ComputeEfficiency --------------------------------------------------------------- + + def save_Event_Nden_eff_ETau(self, tmp_file): + #Save only needed informations for events passing minimal Gen requirements for ETau HLT (passing denominator selection for efficiency) + events = self.get_events() + print(f"Number of events in the file: {len(events)}") + GenLepton = self.get_GenLepton(events) + evt_mask = Denominator_Selection_ETau(GenLepton) + print(f"Number of events with exactly 1 hadronic Tau and one Electron (kind=TauDecayedToHadrons and TauDecayedTo Electron): {ak.sum(evt_mask)}") + self.Save_Event_Nden_Eff(events, GenLepton, evt_mask, tmp_file) + return + + def produceRoot_ETauDeepNet(self, out_file): + #load all events that pass denominator Selection + events = self.get_events() + + # To compute efficiency, we save in denominator GenTau which pass minimal selection + GenTau_mask = hGenTau_selection(events) + GenTaus = get_GenTaus(events) + Tau_Den = GenTaus[GenTau_mask] + + GenElectron_mask = GenElectron_selection(events) + GenElectrons = get_GenElectrons(events) + Electron_Den = GenElectrons[GenElectron_mask] + + mask_den_selection = ( ak.num(Tau_Den['pt']) >=1 ) & ( ak.num(Electron_Den['pt']) >=1 ) + Tau_Den = Tau_Den[mask_den_selection] + events = events[mask_den_selection] + + print(f"Number of GenTaus passing denominator selection: {len(ak.flatten(Tau_Den))}") + + ETau_evt_mask, matchingGentaus_mask, matchingElectrons_mask = evt_sel_HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1(events, is_gen = True) + Tau_Num = (Tau_Den[matchingGentaus_mask])[ETau_evt_mask] + print(f"Number of GenTaus passing numerator selection: {len(ak.flatten(Tau_Num))}") + events_Num = events[ETau_evt_mask] + + self.save_info(events, events_Num, Tau_Den, Tau_Num, out_file) + return + + def produceRoot_ETauPNet(self, out_file, par): + #load all events that pass denominator Selection + events = self.get_events() + + # To compute efficiency, we save in denominator GenTau which pass minimal selection + GenTau_mask = hGenTau_selection(events) + GenTaus = get_GenTaus(events) + Tau_Den = GenTaus[GenTau_mask] + + GenElectron_mask = GenElectron_selection(events) + GenElectrons = get_GenElectrons(events) + Electron_Den = GenElectrons[GenElectron_mask] + + mask_den_selection = ( ak.num(Tau_Den['pt']) >=1 ) & ( ak.num(Electron_Den['pt']) >=1 ) + Tau_Den = Tau_Den[mask_den_selection] + events = events[mask_den_selection] + + print(f"Number of GenTaus passing denominator selection: {len(ak.flatten(Tau_Den))}") + + ETau_evt_mask, matchingGentaus_mask, matchingElectrons_mask = evt_sel_ETau(events, par, is_gen = True) + + Tau_Num = (Tau_Den[matchingGentaus_mask])[ETau_evt_mask] + print(f"Number of GenTaus passing numerator selection: {len(ak.flatten(Tau_Num))}") + + #events = events[SingleTau_evt_mask] + events_Num = events[ETau_evt_mask] + + self.save_info(events, events_Num, Tau_Den, Tau_Num, out_file) + return + + # ------------------------------ functions to Compute Efficiency for opt --------------------------------------------------------------- + + def ComputeEffAlgo_ETauPNet(self, par): + + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + + # Selection of L1/Gen and Jets objects without PNET WP + L1_EG_mask = L1_LooseIsoEG22er2p1_selection(events) | L1_LooseIsoEG24er2p1_selcetion(events) + L1_Tau_mask = L1_IsoTau26er2p1_selection(events) | L1_Tau70er2p1_selection(events) | L1_IsoTau27er2p1_selcetion(events) + Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = False) + GenTau_mask = hGenTau_selection(events) + Electron_mask = Electron_selection(events) + GenElectron_mask = GenElectron_selection(events) + + # at least 1 L1tau/ Jet/ GenTau should pass + ETau_evt_mask = (ak.sum(L1_EG_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Electron_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Egamma = get_L1Egamma(events) + L1Egamma_EG = L1Egamma[L1_EG_mask] + + Jets = get_Jets(events) + Jets_Tau = Jets[Tau_mask] + + Electrons = get_Electrons(events) + Selected_Electrons = Electrons[Electron_mask] + + GenTaus = get_GenTaus(events) + GenTaus_Sel = GenTaus[GenTau_mask] + + GenElectrons = get_GenElectrons(events) + GenElectrons_Electrons = GenElectrons[GenElectron_mask] + + matching_dR_Min0p3_L1 = matching_dR_Min0p3(L1Taus_Tau, L1Egamma_EG, dR_matching_min=0.3) + + ETau_evt_mask = ETau_evt_mask & matching_dR_Min0p3_L1 + + #matching + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + + matchingGenElectrons_mask = matching_GenElectrons(L1Egamma_EG, Selected_Electrons, GenElectrons_Electron) + evt_Electronmask_matching = (ak.sum(matchingGenElectrons_mask, axis=-1) >= 1) + + ETau_evt_mask = ETau_evt_mask & evt_Taumask_matching & evt_Electronmask_matching + + N_den = len(events[ETau_evt_mask]) + print(f"Number of events in denominator: {N_den}") + + # Selection of L1/Gen and Jets objects with PNET WP + Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = True) + # at least 1 L1tau/ Jet/ GenTau should pass + ETau_evt_mask = (ak.sum(L1_EG_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Electron_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + ETau_evt_mask = ETau_evt_mask & matching_dR_Min0p3_L1 + + #matching + Jets_Tau = Jets[Tau_mask] + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + ETau_evt_mask = ETau_evt_mask & evt_Taumask_matching & evt_Electronmask_matching + + N_num = len(events[ETau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def ComputeEffAlgo_ETauDeepNet(self): + + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + + # Selection of L1/Gen and Jets objects without PNET WP + L1_EG_mask = L1_LooseIsoEG22er2p1_selection(events) | L1_LooseIsoEG24er2p1_selcetion(events) + L1_Tau_mask = L1_IsoTau26er2p1_selection(events) | L1_Tau70er2p1_selection(events) | L1_IsoTau27er2p1_selcetion(events) + Tau_mask = Tau_selection_Tau(events, apply_DeepTau_WP = False) + GenTau_mask = hGenTau_selection(events) + Electron_mask = Electron_selection(events) + GenElectron_mask = GenElectron_selection(events) + + # at least 1 L1tau/ Jet/ GenTau should pass + ETau_evt_mask = (ak.sum(L1_EG_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Electron_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Egamma = get_L1Egamma(events) + L1Egamma_EG = L1Egamma[L1_EG_mask] + + Taus = get_Taus(events) + Taus_Tau = Taus[Tau_mask] + + Electrons = get_Electrons(events) + Selected_Electrons = Electrons[Electron_mask] + + GenTaus = get_GenTaus(events) + GenTaus_Sel = GenTaus[GenTau_mask] + + GenElectrons = get_GenElectrons(events) + GenElectrons_Electrons = GenElectrons[GenElectron_mask] + + matching_dR_Min0p3_L1 = matching_dR_Min0p3(L1Taus_Tau, L1Egamma_EG, dR_matching_min=0.3) + + ETau_evt_mask = ETau_evt_mask & matching_dR_Min0p3_L1 + + #matching + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + + matchingGenElectrons_mask = matching_GenElectrons(L1Egamma_EG, Selected_Electrons, GenElectrons_Electron) + evt_Electronmask_matching = (ak.sum(matchingGenElectrons_mask, axis=-1) >= 1) + + ETau_evt_mask = ETau_evt_mask & evt_Taumask_matching & evt_Electronmask_matching + + N_den = len(events[ETau_evt_mask]) + print(f"Number of events in denominator: {N_den}") + + # Selection of L1/Gen and Jets objects with PNET WP + Tau_mask = Tau_selection_Tau(events, apply_DeepTau_WP = True) + # at least 1 L1tau/ Jet/ GenTau should pass + ETau_evt_mask = (ak.sum(L1_EG_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Electron_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + ETau_evt_mask = ETau_evt_mask & matching_dR_Min0p3_L1 + + #matching + Taus_Tau = Taus[Tau_mask] + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + ETau_evt_mask = ETau_evt_mask & evt_Taumask_matching & evt_Electronmask_matching + + N_num = len(events[ETau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + diff --git a/HLTClass/MuTauDataset.py b/HLTClass/MuTauDataset.py new file mode 100644 index 0000000..facb65c --- /dev/null +++ b/HLTClass/MuTauDataset.py @@ -0,0 +1,455 @@ +import awkward as ak +import numpy as np +from HLTClass.dataset import Dataset +from HLTClass.dataset import get_L1Taus, get_Taus, get_Jets, get_GenTaus, hGenTau_selection, matching_Gentaus, matching_L1Taus_obj, compute_PNet_charge_prob, get_L1Muons, get_Muons, get_GenMuons, GenMuon_selection, matching_GenMuons, matching_L1Muons_obj + + +# ------------------------------ functions for MuTau with PNet --------------------------------------------------------------- +def compute_PNet_WP_MuTau(tau_pt, par): + # return PNet WP for MuTau + t1 = par[0] + t2 = par[1] + t3 = 0.001 + t4 = 0 + x1 = 27 + x2 = 100 + x3 = 500 + x4 = 1000 + PNet_WP = tau_pt*0. + ones = tau_pt/tau_pt + PNet_WP = ak.where((tau_pt <= ones*x1) == False, PNet_WP, ones*t1) + PNet_WP = ak.where((tau_pt >= ones*x4) == False, PNet_WP, ones*t4) + PNet_WP = ak.where(((tau_pt < ones*x2) & (tau_pt > ones*x1)) == False, PNet_WP, (t2 - t1) / (x2 - x1) * (tau_pt - ones*x1) + ones*t1) + PNet_WP = ak.where(((tau_pt >= ones*x2) & (tau_pt < ones*x3))== False, PNet_WP, (t3 - t2) / (x3 - x2) * (tau_pt - ones*x2) + ones*t2) + PNet_WP = ak.where(((tau_pt >= ones*x3) & (tau_pt < ones*x4))== False, PNet_WP, (t4 - t3) / (x4 - x3) * (tau_pt - ones*x3) + ones*t3) + return PNet_WP + + +def Jet_selection_Tau(events, par, apply_PNET_WP = True): + # return mask for Jet passing selection for MuTau path + Jet_pt_corr = events['Jet_pt'].compute()*events['Jet_PNet_ptcorr'].compute() + Jets_mask = (events['Jet_pt'].compute() >= 30) & (np.abs(events['Jet_eta'].compute()) <= 2.3) & (Jet_pt_corr >= 30) + if apply_PNET_WP: + probTauP = events['Jet_PNet_probtauhp'].compute() + probTauM = events['Jet_PNet_probtauhm'].compute() + Jets_mask = Jets_mask & ((probTauP + probTauM) >= compute_PNet_WP_MuTau(Jet_pt_corr, par)) & (compute_PNet_charge_prob(probTauP, probTauM) >= 0) + return Jets_mask + +def Muon_selection(events): + Muon_mask = (events['Muon_pt'].compute() >= 20) & (np.abs(events['Muon_eta'].compute()) <= 2.3) & (events['Muon_passSingleMuon'].compute() | events['Muon_passMuTau'].compute()) + return Muon_mask + +def evt_sel_MuTau(events, par, is_gen = False): + # Selection of event passing condition of MuTau with PNet HLT path + mask of objects passing those conditions + + L1_Mu_mask = L1_Mu18er2p1_selection(events) + L1_Tau_mask = L1_Tau24er2p1_OR_Tau26er2p1_selection(events) + + Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = True) + + Muon_mask = Muon_selection(events) + # at least 1 L1tau/ recoJet should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Muons = get_L1Muons(events) + L1Muons_Muon = L1Muons[L1_Mu_mask] + + Jets = get_Jets(events) + Jets_Tau = Jets[Tau_mask] + + Muons = get_Muons(events) + Selected_Muons = Muons[Muon_mask] + + + if is_gen: + # if MC data, at least 1 GenTau should also pass + GenTau_mask = hGenTau_selection(events) + GenMuon_mask = GenMuon_selection(events) + MuTau_evt_mask = MuTau_evt_mask & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) + + # matching + if is_gen: + GenTaus = get_GenTaus(events) + GenTaus_Tau = GenTaus[GenTau_mask] + matchingGentaus_mask = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Tau) + # at least 1 GenTau should match L1Tau/Taus + evt_Taumask_matching = (ak.sum(matchingGentaus_mask, axis=-1) >= 1) + + GenMuons = get_GenMuons(events) + GenMuons_Muon = GenMuons[GenMuon_mask] + matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) + # at least 1 GenTau should match L1Tau/Taus + evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) + else: + matchingJets_mask = matching_L1Taus_obj(L1Taus_Tau, Jets_Tau) + # at least 1 Tau should match L1Tau + evt_Taumask_matching = (ak.sum(matchingJets_mask, axis=-1) >= 1) + + matchingMuons_mask = matching_L1Muons_obj(L1Muons_Muon, Selected_Muons) + evt_Muonmask_matching = (ak.sum(matchingMuons_mask, axis=-1) >= 1) + + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + if is_gen: + return MuTau_evt_mask, matchingGentaus_mask, matchingGenMuons_mask + else: + return MuTau_evt_mask, matchingJets_mask, matchingMuons_mask + +# ------------------------------ functions for DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1 --------------------------------------------------- +def compute_DeepTau_WP_MuTau(tau_pt): + # return DeepTau WP for LooseDeepTauPFTauHPS27_eta2p1 + t1 = 0.649 + t2 = 0.441 + t3 = 0.05 + x1 = 35 + x2 = 100 + x3 = 300 + + Tau_WP = tau_pt*0. + ones = tau_pt/tau_pt + Tau_WP = ak.where((tau_pt <= ones*x1) == False, Tau_WP, ones*t1) + Tau_WP = ak.where((tau_pt >= ones*x3) == False, Tau_WP, ones*t3) + Tau_WP = ak.where(((tau_pt < ones*x2) & (tau_pt > ones*x1)) == False, Tau_WP, (t2 - t1) / (x2 - x1) * (tau_pt - ones*x1) + ones*t1) + Tau_WP = ak.where(((tau_pt >= ones*x2) & (tau_pt < ones*x3))== False, Tau_WP, (t3 - t2) / (x3 - x2) * (tau_pt - ones*x2) + ones*t2) + return Tau_WP + +def Tau_selection_Tau(events, apply_DeepTau_WP = True): + # return mask for Tau passing selection for LooseDeepTauPFTauHPS27_eta2p1 + tau_pt = events['Tau_pt'].compute() + Tau_mask = (tau_pt >= 30) & (np.abs(events['Tau_eta'].compute()) <= 2.1) + if apply_DeepTau_WP: + Tau_mask = Tau_mask & (events['Tau_deepTauVSjet'].compute() >= compute_DeepTau_WP_MuTau(tau_pt)) + return Tau_mask + +def evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is_gen = False): + # Selection of event passing condition of IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1 + mask of objects passing those conditions + + L1_Mu_mask = L1_Mu18er2p1_selection(events) + L1_Tau_mask = L1_Tau24er2p1_OR_Tau26er2p1_selection(events) + + Tau_mask = Tau_selection_Tau(events) + + Muon_mask = Muon_selection(events) + # at least 1 L1tau/ recoJet should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Muons = get_L1Muons(events) + L1Muons_Muon = L1Muons[L1_Mu_mask] + + Taus = get_Taus(events) + Taus_Tau = Taus[Tau_mask] + + Muons = get_Muons(events) + Selected_Muons = Muons[Muon_mask] + + + if is_gen: + # if MC data, at least 1 GenTau should also pass + GenTau_mask = hGenTau_selection(events) + GenMuon_mask = GenMuon_selection(events) + MuTau_evt_mask = MuTau_evt_mask & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) + + # matching + if is_gen: + GenTaus = get_GenTaus(events) + GenTaus_Tau = GenTaus[GenTau_mask] + matchingGentaus_mask = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Tau) + # at least 1 GenTau should match L1Tau/Taus + evt_Taumask_matching = (ak.sum(matchingGentaus_mask, axis=-1) >= 1) + + GenMuons = get_GenMuons(events) + GenMuons_Muon = GenMuons[GenMuon_mask] + matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) + # at least 1 GenTau should match L1Tau/Taus + evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) + else: + matchingTaus_mask = matching_L1Taus_obj(L1Taus_Tau, Taus_Tau) + # at least 1 Tau should match L1Tau + evt_Taumask_matching = (ak.sum(matchingTaus_mask, axis=-1) >= 1) + + matchingMuons_mask = matching_L1Muons_obj(L1Muons_Muon, Selected_Muons) + evt_Muonmask_matching = (ak.sum(matchingMuons_mask, axis=-1) >= 1) + + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + if is_gen: + return MuTau_evt_mask, matchingGentaus_mask, matchingGenMuons_mask + else: + return MuTau_evt_mask, matchingTaus_mask, matchingMuons_mask + +# ------------------------------ Common functions for Mutau path --------------------------------------------------------------- + +def L1_Mu18er2p1_selection(events): + L1_Mu18er2p1 = (events['L1Muon_pt'].compute() >= 18) & (events['L1Muon_eta'].compute() <= 2.131) & (events['L1Muon_eta'].compute() >= -2.131) + return L1_Mu18er2p1 + +def L1_Tau24er2p1_OR_Tau26er2p1_selection(events): + L1_Tau24er2p1_OR_Tau26er2p1 = ((events['L1Tau_pt'].compute() >=24) | (events['L1Tau_pt'].compute() >= 26)) & (events['L1Tau_eta'].compute() <= 2.131) & (events['L1Tau_eta'].compute() >= -2.131) + return L1_Tau24er2p1_OR_Tau26er2p1 + +def Denominator_Selection_MuTau(GenLepton): + # return mask for event passing minimal Gen requirements for MuTau HLT + mask_GenTau = (GenLepton['kind'] == 5) + mask_GenMuon = (GenLepton['kind'] == 4) + ev_mask = (ak.sum(mask_GenTau, axis=-1) == 1) & (ak.sum(mask_GenMuon, axis=-1) == 1) # 1 Gen tau and one Gen Muon should pass this requirements + return ev_mask + +class MuTauDataset(Dataset): + def __init__(self, fileName): + Dataset.__init__(self, fileName) + + # ------------------------------ functions to Compute Rate --------------------------------------------------------------- + + def get_Nnum_Nden_MuTauDeepNet(self): + print(f'Computing rate for HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1:') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + + MuTau_evt_mask, matchingTaus_mask, matchingMuons_mask = evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is_gen = False) + N_num = len(events[MuTau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def get_Nnum_Nden_MuTauPNet(self, par): + print(f'Computing Rate for MuTau path with param: {par}') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + MuTau_evt_mask, matchingJets_mask, matchingMuons_mask = evt_sel_MuTau(events, par, is_gen = False) + N_num = len(events[MuTau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def get_Nnum_Nden_MuTauDeepNet_nPV(self, nPV): + print(f'Computing rate for HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1 for nPV = {nPV}:') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + events = events[events['nPFPrimaryVertex'].compute() == nPV] + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + + MuTau_evt_mask, matchingTaus_mask, matchingMuons_mask = evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is_gen = False) + N_num = len(events[MuTau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def get_Nnum_Nden_MuTauPNet_nPV(self, par, nPV): + print(f'Computing Rate for MuTau path with param: {par} and for nPV = {nPV}:') + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + events = events[events['nPFPrimaryVertex'].compute() == nPV] + N_den = len(events) + print(f"Number of events in denominator: {N_den}") + + MuTau_evt_mask, matchingJets_mask, matchingMuons_mask = evt_sel_MuTau(events, par, is_gen = False) + N_num = len(events[MuTau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + +# ------------------------------ functions for ComputeEfficiency --------------------------------------------------------------- + + def save_Event_Nden_eff_MuTau(self, tmp_file): + #Save only needed informations for events passing minimal Gen requirements for diTau HLT (passing denominator selection for efficiency) + events = self.get_events() + print(f"Number of events in the file: {len(events)}") + GenLepton = self.get_GenLepton(events) + evt_mask = Denominator_Selection_MuTau(GenLepton) + print(f"Number of events with exactly 1 hadronic Tau and one Muon (kind=TauDecayedToHadrons and TauDecayedTo Muon): {ak.sum(evt_mask)}") + self.Save_Event_Nden_Eff(events, GenLepton, evt_mask, tmp_file) + return + + def produceRoot_MuTauDeepNet(self, out_file): + #load all events that pass denominator Selection + events = self.get_events() + + # To compute efficiency, we save in denominator GenTau which pass minimal selection + GenTau_mask = hGenTau_selection(events) + GenTaus = get_GenTaus(events) + Tau_Den = GenTaus[GenTau_mask] + + GenMuon_mask = GenMuon_selection(events) + GenMuons = get_GenMuons(events) + Muon_Den = GenMuons[GenMuon_mask] + + mask_den_selection = ( ak.num(Tau_Den['pt']) >=1 ) & ( ak.num(Muon_Den['pt']) >=1 ) + Tau_Den = Tau_Den[mask_den_selection] + events = events[mask_den_selection] + + print(f"Number of GenTaus passing denominator selection: {len(ak.flatten(Tau_Den))}") + + MuTau_evt_mask, matchingGentaus_mask, matchingMuons_mask = evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is_gen = True) + Tau_Num = (Tau_Den[matchingGentaus_mask])[MuTau_evt_mask] + print(f"Number of GenTaus passing numerator selection: {len(ak.flatten(Tau_Num))}") + events_Num = events[MuTau_evt_mask] + + self.save_info(events, events_Num, Tau_Den, Tau_Num, out_file) + return + + def produceRoot_MuTauPNet(self, out_file, par): + #load all events that pass denominator Selection + events = self.get_events() + + # To compute efficiency, we save in denominator GenTau which pass minimal selection + GenTau_mask = hGenTau_selection(events) + GenTaus = get_GenTaus(events) + Tau_Den = GenTaus[GenTau_mask] + + GenMuon_mask = GenMuon_selection(events) + GenMuons = get_GenMuons(events) + Muon_Den = GenMuons[GenMuon_mask] + + mask_den_selection = ( ak.num(Tau_Den['pt']) >=1 ) & ( ak.num(Muon_Den['pt']) >=1 ) + Tau_Den = Tau_Den[mask_den_selection] + events = events[mask_den_selection] + + print(f"Number of GenTaus passing denominator selection: {len(ak.flatten(Tau_Den))}") + + MuTau_evt_mask, matchingGentaus_mask, matchingMuons_mask = evt_sel_MuTau(events, par, is_gen = True) + + Tau_Num = (Tau_Den[matchingGentaus_mask])[MuTau_evt_mask] + print(f"Number of GenTaus passing numerator selection: {len(ak.flatten(Tau_Num))}") + + events_Num = events[MuTau_evt_mask] + + self.save_info(events, events_Num, Tau_Den, Tau_Num, out_file) + return + + # ------------------------------ functions to Compute Efficiency for opt --------------------------------------------------------------- + + def ComputeEffAlgo_MuTauPNet(self, par): + + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + + L1_Mu_mask = L1_Mu18er2p1_selection(events) + L1_Tau_mask = L1_Tau24er2p1_OR_Tau26er2p1_selection(events) + Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = False) + GenTau_mask = hGenTau_selection(events) + GenMuon_mask = GenMuon_selection(events) + + Muon_mask = Muon_selection(events) + + # at least 1 L1tau/ recoJet should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Muons = get_L1Muons(events) + L1Muons_Muon = L1Muons[L1_Mu_mask] + + Jets = get_Jets(events) + Jets_Tau = Jets[Tau_mask] + + Muons = get_Muons(events) + Selected_Muons = Muons[Muon_mask] + + GenTaus = get_GenTaus(events) + GenTaus_Sel = GenTaus[GenTau_mask] + + GenMuons = get_GenMuons(events) + GenMuons_Muon = GenMuons[GenMuon_mask] + + #matching + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + + matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) + evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) + + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + + N_den = len(events[MuTau_evt_mask]) + print(f"Number of events in denominator: {N_den}") + + # Selection of L1/Gen and Jets objects with PNET WP + Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = True) + # at least 1 L1tau/ Jet/ GenTau should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + + #matching + Jets_Tau = Jets[Tau_mask] + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + + N_num = len(events[MuTau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num + + def ComputeEffAlgo_MuTauDeepNet(self): + + + #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator + events = self.get_events() + + L1_Mu_mask = L1_Mu18er2p1_selection(events) + L1_Tau_mask = L1_Tau24er2p1_OR_Tau26er2p1_selection(events) + Tau_mask = Tau_selection_Tau(events, apply_DeepTau_WP = False) + GenTau_mask = hGenTau_selection(events) + GenMuon_mask = GenMuon_selection(events) + + Muon_mask = Muon_selection(events) + + # at least 1 L1tau/ recoJet should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + + L1Taus = get_L1Taus(events) + L1Taus_Tau = L1Taus[L1_Tau_mask] + + L1Muons = get_L1Muons(events) + L1Muons_Muon = L1Muons[L1_Mu_mask] + + Taus = get_Taus(events) + Taus_Tau = Taus[Tau_mask] + + Muons = get_Muons(events) + Selected_Muons = Muons[Muon_mask] + + GenTaus = get_GenTaus(events) + GenTaus_Sel = GenTaus[GenTau_mask] + + GenMuons = get_GenMuons(events) + GenMuons_Muon = GenMuons[GenMuon_mask] + + #matching + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + + matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) + evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) + + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + + N_den = len(events[MuTau_evt_mask]) + print(f"Number of events in denominator: {N_den}") + + # Selection of L1/Gen and Jets objects with PNET WP + Tau_mask = Tau_selection_Tau(events, apply_DeepTau_WP = True) + # at least 1 L1tau/ Jet/ GenTau should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + + #matching + Taus_Tau = Taus[Tau_mask] + matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Sel) + evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) + + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + + N_num = len(events[MuTau_evt_mask]) + print(f"Number of events in numerator: {N_num}") + print('') + return N_den, N_num diff --git a/HLTClass/dataset.py b/HLTClass/dataset.py index 9988a74..af161dd 100644 --- a/HLTClass/dataset.py +++ b/HLTClass/dataset.py @@ -220,7 +220,10 @@ def Save_Event_Nden_Rate(self, tmp_file, run, lumiSections_range): 'nPFPrimaryVertex', 'nPFSecondaryVertex', 'L1_DoubleIsoTau34er2p1', + 'L1_Mu18er2p1_Tau24er2p1', + 'L1_Mu18er2p1_Tau26er2p1', 'HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1', + 'HLT_IsoMu20_eta2p1_PNetTauhPFJet27_Loose_eta2p3_CrossL1', 'Tau_pt', 'Tau_eta', 'Tau_phi', @@ -237,7 +240,24 @@ def Save_Event_Nden_Rate(self, tmp_file, run, lumiSections_range): 'Jet_PNet_ptcorr', 'Jet_pt', 'Jet_eta', - 'Jet_phi'] + 'Jet_phi', + 'L1Egamma_pt', + 'L1Egamma_eta', + 'L1Egamma_phi', + 'L1Egamma_hwIso', + 'Electron_pt', + 'Electron_eta', + 'Electron_phi', + 'Electron_passSingleElectron', + 'Electron_passETau', + 'L1Muon_pt', + 'L1Muon_eta', + 'L1Muon_phi', + 'Muon_pt', + 'Muon_phi', + 'Muon_eta', + 'Muon_passSingleMuon', + 'Muon_passMuTau'] if len(events)!= 0: print('Saving info in tmp file') diff --git a/Optimisation/helpers_opt.py b/Optimisation/helpers_opt.py index 83d323b..95809be 100644 --- a/Optimisation/helpers_opt.py +++ b/Optimisation/helpers_opt.py @@ -6,4 +6,4 @@ def loss(rate, rate_budget): return 0.5-(1./(2*rate_budget))*rate if rate > (rate_budget + 1.): return 1 - return math.exp(k * (rate - rate_budget)) - 1 \ No newline at end of file + return math.exp(k * (rate - rate_budget)) - 1 diff --git a/Optimisation/merger.py b/Optimisation/merger.py new file mode 100644 index 0000000..909355f --- /dev/null +++ b/Optimisation/merger.py @@ -0,0 +1,25 @@ +import json +import os +import pandas as pd + +#p = "/eos/home-j/jleonhol/www/PNetAtHLT/2024_v1/PNetAtHLT/Optimisation/result/" + +#p = "/eos/home-j/jleonhol/www/PNetAtHLT/2024_v1/PNetAtHLT/OptimisationDiTau/result_new/" + +p = "/afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/Optimisation/result/" + +#p = "/eos/home-j/jleonhol/www/PNetAtHLT/2024_v1/PNetAtHLT/OptimisationDiTauSingleTau/result/" + +# p = "/eos/home-j/jleonhol/www/PNetAtHLT/2024_v1/PNetAtHLT/OptimisationDiTauSingleTauDenDouble/result/" + +l = [] +for f in os.listdir(p): + with open(p + f) as f: + l.append(json.load(f)) + +df = pd.DataFrame(l) + +#df.to_pickle("results_ditaujet.pickle") +df.to_pickle("results_mutau.pickle") +#df.to_pickle("results_ditau_singletau.pickle") +# df.to_pickle("results_ditau_singletau_doubleden.pickle") diff --git a/Optimisation/plot_manualDiTau.py b/Optimisation/plot_manualDiTau.py new file mode 100644 index 0000000..f500646 --- /dev/null +++ b/Optimisation/plot_manualDiTau.py @@ -0,0 +1,106 @@ +import pandas as pd +import json +import numpy as np + +def plot(xaxis, yaxis, parameters, x_title, y_title, min_x, max_x, min_y, max_y, save_path, + data=None, params_to_mark=[], plot_text=False): + import matplotlib + matplotlib.use("PDF") + from matplotlib import pyplot as plt + plt.figure() + ax = plt.subplot() + #if min_x and max_x: + # ax.set_xlim(min_x, max_x) + #if min_y and max_y: + # ax.set_ylim(min_y, max_y) + # data_final = data[data["params"].astype(str) == "[26, 0.64, 0.46]"] + try: + ax = data.plot(x=xaxis, y=yaxis, kind="scatter", ax=ax) + if min_x and max_x: + ax.set_xlim(min_x, max_x) + if min_y and max_y: + ax.set_ylim(min_y, max_y) + except: + plt.scatter(xaxis, yaxis, marker="o") + for (x, y, label) in zip(data[xaxis], data[yaxis], parameters): + if label in params_to_mark: + plt.scatter(x, y, marker="o", color="r") + if not plot_text and label not in params_to_mark: + continue + plt.annotate(label, # this is the text + (x, y), # this is the point to label + textcoords="offset points", # how to position the text + xytext=(0, 10), # distance from text to points (x,y) + ha='center', # horizontal alignment can be left, right or center + size=5) + plt.xlabel(x_title) + plt.ylabel(y_title) + + x_text=0.05 + y_text=0.9 + plt.text(x_text, 1.02, "CMS", fontsize='large', fontweight='bold', + transform=ax.transAxes) + upper_text = "private work" + plt.text(x_text + 0.1, 1.02, upper_text, transform=ax.transAxes) + # text = [self.dataset.process.label.latex, self.category.label.latex] + # for t in text: + # plt.text(x_text, y_text, t, transform=ax.transAxes) + # y_text -= 0.05 + + plt.savefig(save_path) + plt.close('all') + +if __name__ == '__main__': + + # with open("results_tri_optimised.json") as f: + # d = json.load(f) + + + # import glob + # read_files = glob.glob("/eos/home-j/jleonhol/www/PNetAtHLT/PNetAtHLT/Optimisation/result/*.json") + # d = [] + # for f in read_files: + # with open(f, "rb") as infile: + # d.append(json.load(infile)) + + df = pd.read_pickle("/afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/Optimisation/results_mutau.pickle") + + print(df) + + # df = pd.DataFrame(d) + # df.convert_dtypes() + + # with open("results.json") as f: + # d_noditaujet = json.load(f) + + # # # # with open("results_deeptau.json") as f: + # # # # dd = json.load(f) + + # # # # results = { + # # # # "method": [ + # # # # "DeepTau", + # # # # "DeepTau (no DiTauJet)", + # # # # "PNet (no DiTauJet)", + # # # # ], + # # # # "rate": [ + # # # # dd["rate"], + # # # # dd["noditaujet_rate"], + # # # # df["rate_noditaujet"][0], + # # # # ], + # # # # "efficiency": [ + # # # # dd["efficiency"], + # # # # dd["noditaujet_eff"], + # # # # df["eff_noditaujet"][0], + # # # # ] + # # # # } + # # # # results = pd.DataFrame(results) + # # # # print(results) + + # plot("rate", "efficiency", df["params"], "Rate (Hz)", "Efficiency", 17.80, 18.20, 0.84, 0.86, "plot.pdf", data=df) + # plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 60, 61, 0.675, 0.68, "plot_tri_pt_optimised_cut.pdf", data=df) + # plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 60, 61, 0.679, 0.68, "plot_optimised_cut.pdf", data=df) + plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 0, 100, 0., 1., "plot_optimised_ditau_new.pdf", data=df, params_to_mark=[["deeptau"], [0.56, 0.47], [0.6, 0.5]]) + plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 48, 50, 0.7975, 0.81, "plot_optimised_ditau_medium_new.pdf", data=df, params_to_mark=[[0.56, 0.47]], plot_text=True) + plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 39, 40, 0.77, 0.79, "plot_optimised_ditau_tight_new.pdf", data=df, params_to_mark=[[0.60, 0.50]], plot_text=True) + # plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 0., 1., 0., 1., "plot_tri_pt_optimised_cut.pdf", data=df) + #plot([15], [0.95], [[1,2]], "Rate (Hz)", "Efficiency", 10, 20, 0.9, 1., "plot.pdf") diff --git a/Optimisation/run_optETau.py b/Optimisation/run_optETau.py new file mode 100644 index 0000000..0cfb30f --- /dev/null +++ b/Optimisation/run_optETau.py @@ -0,0 +1,61 @@ +from helpers import files_from_path, load_cfg_file +import os +from scipy.optimize import minimize +from HLTClass.ETauDataset import ETauDataset +import numpy as np +from Optimisation.helpers_opt import loss + +def run_optimization(dataset_rate, dataset_eff, rate_budget): + config = load_cfg_file() + PNetparam_t3 = float(config['OPT']['PNet_t3']) + + def f(a): + a = np.append(a, [PNetparam_t3]) + + N_den, N_num = dataset_rate.get_Nnum_Nden_ETauPNet(a) + rate = (N_num/N_den)*L1A_physics + + N_den, N_num = dataset_eff.ComputeEffAlgo_ETauPNet(a) + eff_algo = (N_num/N_den) + + print(f'Rate: {rate}') + print(f'Algo Eff: {eff_algo}') + print(f'Output score: {- eff_algo + loss(rate, rate_budget)}') + print('---------------------------------------------------------------------------------------------------------') + return - eff_algo + loss(rate, rate_budget) + + res = minimize(f, [0.95, 0.95], bounds=((0.8, 0.99), (0.8, 0.99)), method="L-BFGS-B", options={"eps": [0.01, 0.01]}) + return res.x + +if __name__ == '__main__': + + config = load_cfg_file() + RefRun = int(config['RUNINFO']['ref_run']) + HLT_name = config['HLT']['HLTname'] + L1A_physics = float(config['RUNINFO']['L1A_physics']) + Rate_path = os.path.join(config['DATA']['RateDenPath'], f'Run_{RefRun}') + Eff_path = os.path.join(config['DATA']['EffDenPath'], HLT_name) + PNetparam_t3 = float(config['OPT']['PNet_t3']) # optimise only t1 and t2 + + FileNameList_eff = f"/afs/cern.ch/work/s/skeshri/TauHLT/Braden/TauHLTOptimzation/PnetAtHLT/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" # optimisation with ZprimeToTauTau_M-4000 data + FileNameList_rate = files_from_path(Rate_path)[0] # only one otherwise it's too long (gives good aprosimation for the rate) + + dataset_eff = ETauDataset(FileNameList_eff) + dataset_rate = ETauDataset(FileNameList_rate) + + optim_param = run_optimization(dataset_rate, dataset_eff, rate_budget = 15) + optim_param = np.append(optim_param, [PNetparam_t3]) + + print('Optimisation finished:') + print("Optimized parameters:", optim_param) + + N_den, N_num = dataset_rate.get_Nnum_Nden_ETauPNet(optim_param) + rate_opt = (N_num/N_den)*L1A_physics + + print("Approx. Rate:", rate_opt) + + N_den, N_num = dataset_eff.ComputeEffAlgo_ETauPNet(optim_param) + algo_eff = (N_num/N_den) + + print("Approx. Algo Eff:", algo_eff) + diff --git a/Optimisation/tasks.py b/Optimisation/tasks.py new file mode 100644 index 0000000..ba35ff0 --- /dev/null +++ b/Optimisation/tasks.py @@ -0,0 +1,397 @@ +#!/usr/bin/env python +import law +import os +import shutil +import json + +from law_customizations import Task, HTCondorWorkflow +from helpers import files_from_path, load_cfg_file, hadd_anatuple +from HLTClass.MuTauDataset import MuTauDataset +#from HLTClass.DiTauJetDataset import DiTauJetDataset +#from HLTClass.DiTauDataset import DiTauDataset +#from HLTClass.DoubleORSingleTauDataset import DoubleORSingleTauDataset +#from Optimisation.run_fullcombo import Threshold_optimiser +#from Optimisation.run_fulldeeptau import Threshold_optimiser as Th_opt_deeptau + +""" +class SaveOptimisationResults(Task, HTCondorWorkflow, law.LocalWorkflow): + ''' + Produce root file where Events passing denominator selection are saved + ''' + + config = load_cfg_file() + output_path = config['DATA']['result_opt'] + + def create_branch_map(self): + branches = {} + os.makedirs(self.output_path, exist_ok=True) + + max_param_1 = { + 26: 0.7, + 27: 0.65, + 28: 0.6, + 29: 0.55, + 30: 0.5 + } + + max_param_2 = { + # 26: 0.7, + 26: 0.5, + 27: 0.65, + 28: 0.6, + 29: 0.55, + 30: 0.5 + } + + min_param_1 = { + # 26: 0.5, + 26: 0.6, + 27: 0.45, + 28: 0.4, + 29: 0.35, + 30: 0.3 + } + + # min_param_2 = 0.1 + min_param_2 = 0.4 + + # min_param_1 = 0.2 + # # min_param_1 = 0.59 + # min_param_2 = 0.2 + # # min_param_2 = 0.59 + # max_param = 1. + + step = 0.01 + pt_cuts = [26, 27, 28] + # pt_cuts = [26, 27, 28, 29, 30] + index = 0 + for pt in pt_cuts: + for i in range(int(round((max_param_1[pt] - min_param_1[pt]) / step)) + 1): + for j in range(int(round((max_param_2[pt] - min_param_2) / step)) + 1): + param = (round(min_param_1[pt] + i * step, 2), round(min_param_2 + j * step, 2)) + if param[1] > param[0]: + continue + branches[index] = (pt, param[0], param[1]) + index += 1 + + branches[index] = ("deeptau", ) + return branches + + def output(self): + p = self.branch_data + path = os.path.join(self.output_path, f'results_{"_".join([str(elem) for elem in p])}.json') + return law.LocalFileTarget(path) + + def run(self): + HLT_name = "HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60" + L1A_physics = float(self.config['RUNINFO']['L1A_physics']) + + RefRun = list(filter(None, (x.strip() for x in self.config['RUNINFO']['ref_run'].splitlines()))) + tag = f'Run_{RefRun[0]}' if len(RefRun) == 1 else f'Run_{"_".join(RefRun)}' + Rate_path = os.path.join(self.config['DATA']['RateDenPath'], tag) + Eff_path = os.path.join(self.config['DATA']['EffDenPath'], HLT_name) + + # FileNameList_eff = f"/eos/user/j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" + FileNameList_eff = [ + f"/eos/user/j/jleonhol/www/PNetAtHLT/2024_v1/data/EfficiencyDen/{HLT_name}/VBFHToTauTau_M125.root" + ] + FileNameList_rate = files_from_path(Rate_path) + FileNameList_rate = FileNameList_rate[0] + + params = self.branch_map[self.branch] + + if "deeptau" not in params: + th = Threshold_optimiser() + else: + th = Th_opt_deeptau() + th.dataset_eff = DiTauJetDataset(FileNameList_eff) + th.dataset_rate = DiTauJetDataset(FileNameList_rate) + + if "deeptau" not in params: + eff, rate = th.f((params[1], params[2]), params[0]) + eff_noditaujet, rate_noditaujet = th.f((params[1], params[2]), params[0], False) + else: + eff, rate = th.f() + eff_noditaujet, rate_noditaujet = th.f(False) + + d = { + "params": params, "eff": eff, "rate": rate, + "eff_noditaujet": eff_noditaujet, "rate_noditaujet": rate_noditaujet + } + + with open(self.output().path, "w+") as out: + json.dump(d, out) + + +class SaveOptimisationResultsDiTau(SaveOptimisationResults): + ''' + Produce root file where Events passing denominator selection are saved + ''' + + config = load_cfg_file() + output_path = config['DATA']['result_opt_ditau'] + add_deeptau = True + def create_branch_map(self): + branches = {} + os.makedirs(self.output_path, exist_ok=True) + + param_ranges = [ + (0.4, 0.65), + (0.4, 0.55) + ] + + # min_param_1 = 0.2 + # # min_param_1 = 0.59 + # min_param_2 = 0.2 + # # min_param_2 = 0.59 + # max_param = 1. + + step = 0.01 + # pt_cuts = [26, 27, 28, 29, 30] + index = 0 + for i in range(int(round((param_ranges[0][1] - param_ranges[0][0]) / step)) + 1): + for j in range(int(round((param_ranges[1][1] - param_ranges[1][0]) / step)) + 1): + param = (round(param_ranges[0][0] + i * step, 2), round(param_ranges[1][0] + j * step, 2)) + if param[1] > param[0]: + continue + branches[index] = (param[0], param[1]) + index += 1 + + if self.add_deeptau: + branches[index] = ("deeptau", ) + + return branches + + def run(self): + config = load_cfg_file() + HLT_name = "HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1" + L1A_physics = float(self.config['RUNINFO']['L1A_physics']) + + RefRun = list(filter(None, (x.strip() for x in self.config['RUNINFO']['ref_run'].splitlines()))) + tag = f'Run_{RefRun[0]}' if len(RefRun) == 1 else f'Run_{"_".join(RefRun)}' + Rate_path = os.path.join(self.config['DATA']['RateDenPath'], tag) + Eff_path = os.path.join(self.config['DATA']['EffDenPath'], HLT_name) + + # FileNameList_eff = f"/eos/user/j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" + FileNameList_eff = [ + f"/eos/user/j/jleonhol/www/PNetAtHLT/2024_v1/data/EfficiencyDen/{HLT_name}/VBFHToTauTau_M125.root" + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/GluGlutoHHto2B2Tau_kl-1p00_kt-1p00_c2-0p00.root", + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/VBFHHto2B2Tau_CV-1_C2V-1_C3-1.root", + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/VBFHToTauTau_M125.root", + ] + FileNameList_rate = files_from_path(Rate_path) + FileNameList_rate = FileNameList_rate[0] + dataset_eff = DiTauDataset(FileNameList_eff) + dataset_rate = DiTauDataset(FileNameList_rate) + + params = self.branch_map[self.branch] + if not "deeptau" in params: + N_den, N_num = dataset_rate.get_Nnum_Nden_DiTauPNet(params) + rate = (N_num/N_den)*L1A_physics + + N_den, N_num = dataset_eff.ComputeEffAlgo_DiTauPNet(params) + eff = (N_num/N_den) + else: + N_den, N_num = dataset_rate.get_Nnum_Nden_HLTDoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1() + rate = (N_num/N_den)*L1A_physics + + N_den, N_num = dataset_eff.ComputeEffAlgo_HLTDoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1() + eff = (N_num/N_den) + + d = { + "params": params, "eff": eff, "rate": rate, + } + + with open(self.output().path, "w+") as out: + json.dump(d, out) + + +class SaveOptimisationResultsDiTauSingleTau(SaveOptimisationResultsDiTau): + ''' + Produce root file where Events passing denominator selection are saved + ''' + + config = load_cfg_file() + output_path = config['DATA']['result_opt_ditau_singletau'] + add_deeptau = False + par_double = [0.57, 0.44, 0] + + den_double = False + + def create_branch_map(self): + branches = {} + os.makedirs(self.output_path, exist_ok=True) + + param_ranges = [ + (0.9, 1.1), + (0.8, 1.) + ] + + # min_param_1 = 0.2 + # # min_param_1 = 0.59 + # min_param_2 = 0.2 + # # min_param_2 = 0.59 + # max_param = 1. + + step = 0.01 + # pt_cuts = [26, 27, 28, 29, 30] + index = 0 + for i in range(int(round((param_ranges[0][1] - param_ranges[0][0]) / step)) + 1): + for j in range(int(round((param_ranges[1][1] - param_ranges[1][0]) / step)) + 1): + param = (round(param_ranges[0][0] + i * step, 2), round(param_ranges[1][0] + j * step, 2)) + if param[1] > param[0]: + continue + branches[index] = (param[0], param[1]) + index += 1 + + if self.add_deeptau: + branches[index] = ("deeptau", ) + + return branches + + def run(self): + HLT_name = "HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1" + L1A_physics = float(self.config['RUNINFO']['L1A_physics']) + + RefRun = list(filter(None, (x.strip() for x in self.config['RUNINFO']['ref_run'].splitlines()))) + tag = f'Run_{RefRun[0]}' if len(RefRun) == 1 else f'Run_{"_".join(RefRun)}' + Rate_path = os.path.join(self.config['DATA']['RateDenPath'], tag) + Eff_path = os.path.join(self.config['DATA']['EffDenPath'], HLT_name) + + # FileNameList_eff = f"/eos/user/j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" + FileNameList_eff = [ + f"/eos/user/j/jleonhol/www/PNetAtHLT/2024_v1/data/EfficiencyDen/{HLT_name}/VBFHToTauTau_M125.root" + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/GluGlutoHHto2B2Tau_kl-1p00_kt-1p00_c2-0p00.root", + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/VBFHHto2B2Tau_CV-1_C2V-1_C3-1.root", + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/VBFHToTauTau_M125.root", + ] + FileNameList_rate = files_from_path(Rate_path) + FileNameList_rate = FileNameList_rate[0] + dataset_eff = DoubleORSingleTauDataset(FileNameList_eff) + dataset_rate = DoubleORSingleTauDataset(FileNameList_rate) + + params = self.branch_map[self.branch] + if not "deeptau" in params: + N_den, N_num = dataset_rate.get_Nnum_Nden_DoubleORSinglePNet(params) + rate = (N_num/N_den)*L1A_physics + + N_den, N_num = dataset_eff.ComputeEffAlgo_DoubleORSinglePNet(params) + eff = (N_num/N_den) + else: + raise ValueError("Not implemented") + N_den, N_num = dataset_rate.get_Nnum_Nden_HLTDoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1() + rate = (N_num/N_den)*L1A_physics + + N_den, N_num = dataset_eff.ComputeEffAlgo_HLTDoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1() + eff = (N_num/N_den) + + d = { + "params": params, "eff": eff, "rate": rate, + } + + with open(self.output().path, "w+") as out: + json.dump(d, out) + + +class SaveOptimisationResultsDiTauSingleTauDenDouble(SaveOptimisationResultsDiTauSingleTau): + config = load_cfg_file() + output_path = config['DATA']['result_opt_ditau_singletau_dendouble'] + add_deeptau = False + par_double = [0.57, 0.44, 0] + den_double = True +""" + + + +class SaveOptimisationResultsMuTau(Task, HTCondorWorkflow, law.LocalWorkflow): + ''' + Produce root file where Events passing denominator selection are saved + ''' + + config = load_cfg_file() + output_path = config['DATA']['result_opt_Etau'] + add_deeptau = False + par_double = [0.57, 0.44, 0] + + den_double = False + + def create_branch_map(self): + branches = {} + os.makedirs(self.output_path, exist_ok=True) + + param_ranges = [ + (0.45, 0.60), + (0.35, 0.50) + ] + + # min_param_1 = 0.2 + # # min_param_1 = 0.59 + # min_param_2 = 0.2 + # # min_param_2 = 0.59 + # max_param = 1. + + step = 0.01 + # pt_cuts = [26, 27, 28, 29, 30] + index = 0 + for i in range(int(round((param_ranges[0][1] - param_ranges[0][0]) / step)) + 1): + for j in range(int(round((param_ranges[1][1] - param_ranges[1][0]) / step)) + 1): + param = (round(param_ranges[0][0] + i * step, 2), round(param_ranges[1][0] + j * step, 2)) + if param[1] > param[0]: + continue + branches[index] = (param[0], param[1]) + index += 1 + + if self.add_deeptau: + branches[index] = ("deeptau", ) + + return branches + + def run(self): + #HLT_name = "HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1" + HLT_name = "HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1" + L1A_physics = float(self.config['RUNINFO']['L1A_physics']) + + RefRun = list(filter(None, (x.strip() for x in self.config['RUNINFO']['ref_run'].splitlines()))) + tag = f'Run_{RefRun[0]}' if len(RefRun) == 1 else f'Run_{"_".join(RefRun)}' + Rate_path = os.path.join(self.config['DATA']['RateDenPath'], tag) + Eff_path = os.path.join(self.config['DATA']['EffDenPath'], HLT_name) + + # FileNameList_eff = f"/eos/user/j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/{HLT_name}/ZprimeToTauTau_M-4000.root" + FileNameList_eff = [ + f"/afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/EfficiencyDen/{HLT_name}/GluGlutoHHto2B2Tau_kl-1p00_kt-1p00_c2-3p00.root" + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/GluGlutoHHto2B2Tau_kl-1p00_kt-1p00_c2-0p00.root", + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/VBFHHto2B2Tau_CV-1_C2V-1_C3-1.root", + # "/eos/home-j/jleonhol/www/PNetAtHLT/data/EfficiencyDen/HLT_DoubleMediumDeepTauPFTauHPS30_L2NN_eta2p1_PFJet60/VBFHToTauTau_M125.root", + ] + FileNameList_rate = files_from_path(Rate_path) + FileNameList_rate = FileNameList_rate[0] + dataset_eff = MuTauDataset(FileNameList_eff) + dataset_rate = MuTauDataset(FileNameList_rate) + + params = self.branch_map[self.branch] + if not "deeptau" in params: + N_den, N_num = dataset_rate.get_Nnum_Nden_MuTauPNet( params) + rate = (N_num/N_den)*L1A_physics + + N_den, N_num = dataset_eff.ComputeEffAlgo_MuTauPNet(params) + eff = (N_num/N_den) + else: + raise ValueError("Not implemented") + N_den, N_num = dataset_rate.get_Nnum_Nden_MuTauDeepNet() + rate = (N_num/N_den)*L1A_physics + + N_den, N_num = dataset_eff.ComputeEffAlgo_MuTauDeepNet() + eff = (N_num/N_den) + + d = { + "params": params, "eff": eff, "rate": rate, + } + + with open(self.output().path, "w+") as out: + json.dump(d, out) + + + def output(self): + p = self.branch_data + path = os.path.join(self.output_path, f'results_{"_".join([str(elem) for elem in p])}.json') + return law.LocalFileTarget(path) From d203a3960647ea95e786ffb0cffc7b167c3e06bf Mon Sep 17 00:00:00 2001 From: skeshri Date: Thu, 25 Apr 2024 11:08:47 +0200 Subject: [PATCH 4/5] update for MuTau optimasation --- HLTClass/MuTauDataset.py | 29 ++++++++++++++++++++++++++--- Optimisation/plot_manualDiTau.py | 2 +- Optimisation/tasks.py | 5 +++-- config.ini | 32 +++++++++++++++++++++----------- env.sh | 5 +++-- law.cfg | 1 + 6 files changed, 55 insertions(+), 19 deletions(-) diff --git a/HLTClass/MuTauDataset.py b/HLTClass/MuTauDataset.py index facb65c..1fd7f59 100644 --- a/HLTClass/MuTauDataset.py +++ b/HLTClass/MuTauDataset.py @@ -2,6 +2,16 @@ import numpy as np from HLTClass.dataset import Dataset from HLTClass.dataset import get_L1Taus, get_Taus, get_Jets, get_GenTaus, hGenTau_selection, matching_Gentaus, matching_L1Taus_obj, compute_PNet_charge_prob, get_L1Muons, get_Muons, get_GenMuons, GenMuon_selection, matching_GenMuons, matching_L1Muons_obj +from helpers import delta_r + + +def matching_dR_Min0p5(Obj1, Obj2, dR_matching_min = 0.5): + obj1_inpair, obj2_inpair = ak.unzip(ak.cartesian([Obj2, Obj1], nested=True)) + dR_obj2_obj1 = delta_r(obj2_inpair, obj1_inpair) + mask_obj2_obj1 = (dR_obj2_obj1 > dR_matching_min) + + mask = ak.any(mask_obj2_obj1, axis=-1) + return mask # ------------------------------ functions for MuTau with PNet --------------------------------------------------------------- @@ -28,7 +38,7 @@ def compute_PNet_WP_MuTau(tau_pt, par): def Jet_selection_Tau(events, par, apply_PNET_WP = True): # return mask for Jet passing selection for MuTau path Jet_pt_corr = events['Jet_pt'].compute()*events['Jet_PNet_ptcorr'].compute() - Jets_mask = (events['Jet_pt'].compute() >= 30) & (np.abs(events['Jet_eta'].compute()) <= 2.3) & (Jet_pt_corr >= 30) + Jets_mask = (events['Jet_pt'].compute() >= 27) & (np.abs(events['Jet_eta'].compute()) <= 2.3) & (Jet_pt_corr >= 27) if apply_PNET_WP: probTauP = events['Jet_PNet_probtauhp'].compute() probTauM = events['Jet_PNet_probtauhm'].compute() @@ -63,6 +73,8 @@ def evt_sel_MuTau(events, par, is_gen = False): Muons = get_Muons(events) Selected_Muons = Muons[Muon_mask] + #overlap_mask1 = matching_dR_Min0p5(L1Muons_Muon, L1Taus_Tau) + #overlap_mask2 = matching_dR_Min0p5(Selected_Muons, Jets_Tau) if is_gen: # if MC data, at least 1 GenTau should also pass @@ -81,6 +93,8 @@ def evt_sel_MuTau(events, par, is_gen = False): GenMuons = get_GenMuons(events) GenMuons_Muon = GenMuons[GenMuon_mask] matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) + #overlap_mask3 = matching_dR_Min0p5(GenMuons_Muon,GenTaus_Tau) + #overlap_mask = overlap_mask1 & overlap_mask2 & overlap_mask3 # at least 1 GenTau should match L1Tau/Taus evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) else: @@ -90,8 +104,10 @@ def evt_sel_MuTau(events, par, is_gen = False): matchingMuons_mask = matching_L1Muons_obj(L1Muons_Muon, Selected_Muons) evt_Muonmask_matching = (ak.sum(matchingMuons_mask, axis=-1) >= 1) + #overlap_mask = overlap_mask1 & overlap_mask2 - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + #MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching if is_gen: return MuTau_evt_mask, matchingGentaus_mask, matchingGenMuons_mask else: @@ -147,6 +163,9 @@ def evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is Muons = get_Muons(events) Selected_Muons = Muons[Muon_mask] + #overlap_mask1 = matching_dR_Min0p5(L1Muons_Muon, L1Taus_Tau) + #overlap_mask2 = matching_dR_Min0p5(Selected_Muons, Taus_Tau) + if is_gen: # if MC data, at least 1 GenTau should also pass @@ -165,6 +184,8 @@ def evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is GenMuons = get_GenMuons(events) GenMuons_Muon = GenMuons[GenMuon_mask] matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) + #overlap_mask3 = matching_dR_Min0p5(GenMuons_Muon,GenTaus_Tau) + #overlap_mask = overlap_mask1 & overlap_mask2 & overlap_mask3 # at least 1 GenTau should match L1Tau/Taus evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) else: @@ -174,8 +195,10 @@ def evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is matchingMuons_mask = matching_L1Muons_obj(L1Muons_Muon, Selected_Muons) evt_Muonmask_matching = (ak.sum(matchingMuons_mask, axis=-1) >= 1) + #overlap_mask = overlap_mask1 & overlap_mask2 - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + #MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching if is_gen: return MuTau_evt_mask, matchingGentaus_mask, matchingGenMuons_mask else: diff --git a/Optimisation/plot_manualDiTau.py b/Optimisation/plot_manualDiTau.py index f500646..c07fc87 100644 --- a/Optimisation/plot_manualDiTau.py +++ b/Optimisation/plot_manualDiTau.py @@ -99,7 +99,7 @@ def plot(xaxis, yaxis, parameters, x_title, y_title, min_x, max_x, min_y, max_y, # plot("rate", "efficiency", df["params"], "Rate (Hz)", "Efficiency", 17.80, 18.20, 0.84, 0.86, "plot.pdf", data=df) # plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 60, 61, 0.675, 0.68, "plot_tri_pt_optimised_cut.pdf", data=df) # plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 60, 61, 0.679, 0.68, "plot_optimised_cut.pdf", data=df) - plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 0, 100, 0., 1., "plot_optimised_ditau_new.pdf", data=df, params_to_mark=[["deeptau"], [0.56, 0.47], [0.6, 0.5]]) + plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 0, 100, 0., 1., "plot_optimised_ditau_new.pdf", data=df, params_to_mark=[["deeptau"],[0.48,0.4],[0.52,0.42],[0.56,0.47]]) plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 48, 50, 0.7975, 0.81, "plot_optimised_ditau_medium_new.pdf", data=df, params_to_mark=[[0.56, 0.47]], plot_text=True) plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 39, 40, 0.77, 0.79, "plot_optimised_ditau_tight_new.pdf", data=df, params_to_mark=[[0.60, 0.50]], plot_text=True) # plot("rate", "eff", df["params"], "Rate (Hz)", "Efficiency", 0., 1., 0., 1., "plot_tri_pt_optimised_cut.pdf", data=df) diff --git a/Optimisation/tasks.py b/Optimisation/tasks.py index ba35ff0..95030d9 100644 --- a/Optimisation/tasks.py +++ b/Optimisation/tasks.py @@ -310,7 +310,7 @@ class SaveOptimisationResultsMuTau(Task, HTCondorWorkflow, law.LocalWorkflow): config = load_cfg_file() output_path = config['DATA']['result_opt_Etau'] - add_deeptau = False + add_deeptau = True par_double = [0.57, 0.44, 0] den_double = False @@ -319,6 +319,7 @@ def create_branch_map(self): branches = {} os.makedirs(self.output_path, exist_ok=True) + param_ranges = [ (0.45, 0.60), (0.35, 0.50) @@ -376,7 +377,7 @@ def run(self): N_den, N_num = dataset_eff.ComputeEffAlgo_MuTauPNet(params) eff = (N_num/N_den) else: - raise ValueError("Not implemented") + #raise ValueError("Not implemented") N_den, N_num = dataset_rate.get_Nnum_Nden_MuTauDeepNet() rate = (N_num/N_den)*L1A_physics diff --git a/config.ini b/config.ini index ec3644e..ead2d11 100644 --- a/config.ini +++ b/config.ini @@ -1,16 +1,24 @@ [RUNINFO] -ref_run = 362617 -LumiSectionsRange_low = 0 -LumiSectionsRange_up = 245 +#ref_run = 362617 +#LumiSectionsRange_low = 0 +#LumiSectionsRange_up = 245 +#Area = 2023D +## for run ref_run and lumisectionRange, OMS values for the L1 rate (L1A physics) +#L1A_physics = 91374.04 + +ref_run = 370293 +LumiSectionsRange_low = 174 +LumiSectionsRange_up = 265 Area = 2023D -# for run ref_run and lumisectionRange, OMS values for the L1 rate (L1A physics) -L1A_physics = 91374.04 +L1A_physics = 96216.75 + [HLT] # HLT name to study #HLTname = HLT_DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1 -#HLTname = HLT_Ele24_eta2p1_WPTight_Gsf_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1 +#HLTname = HLT_Ele24_eta2p1_WPTight_Gsf_LooseDeepTauPFTauHPS30_eta2p1_CrossL1 HLTname = HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1 +#HLTname = HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1 #HLTname = HLT_LooseDeepTauPFTauHPS180_L2NN_eta2p1_v3 # rate in CMS OMS for HLTname (for run ref_run and within lumisectionRange) HLT_rate = 50.12 @@ -22,7 +30,7 @@ SamplesPath = /eos/cms/store/group/phys_tau/TauTrigger/Run3_HLT/prod_2024_v1 #SamplesPath = /afs/cern.ch/user/s/skeshri # .. For rate computation -number_of_ephemeral_folder = 2 +number_of_ephemeral_folder = 8 # .. For eff computation MCDataFolderNames = # GluGluHToTauTau_M-125 @@ -32,6 +40,8 @@ MCDataFolderNames = #result_eff = /afs/cern.ch/user/p/pdebryas/PnetAtHLT/PnetAtHLT/ComputeEfficiency/result/ result_rate = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/result/ result_eff = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/result/ + +result_opt_Etau = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/Optimisation/result # path where to store anatuples (events which pass denominator selection in eff/rate computation) EffDenPath = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/EfficiencyDen/ RateDenPath = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/RateDen/ @@ -40,10 +50,10 @@ tmpPath = /afs/cern.ch/work/s/skeshri/TauHLT/Braden/Forked/TauTriggerDev/PnetAtH [OPT] # Use of current DeepTau WP: PNet_mode = false; else PNet_mode = true -PNet_mode = true +PNet_mode = false # Pnet_tauhm+Pnet_tauhp > PNet_WP(PNet_t1, PNet_t2) -PNet_t1 = 0.56 +PNet_t1 = 0.48 # Pnet_tauhm+Pnet_tauhp > PNet_WP(PNet_t1, PNet_t2) -PNet_t2 = 0.47 +PNet_t2 = 0.4 # Pnet_chargeprob > PNet_t3 -PNet_t3 = 30 +PNet_t3 = 0.001 diff --git a/env.sh b/env.sh index cd5190e..7b568aa 100644 --- a/env.sh +++ b/env.sh @@ -12,7 +12,7 @@ action() { export RUN_PATH="$this_dir" #setup private conda installation for env activation - PRIVATE_CONDA_INSTALL=/afs/cern.ch/work/p/pdebryas/miniconda3 + PRIVATE_CONDA_INSTALL=/eos/user/s/skeshri/miniconda3 __conda_setup="$($PRIVATE_CONDA_INSTALL/bin/conda shell.${SHELL##*/} hook)" if [ $? -eq 0 ]; then eval "$__conda_setup" @@ -26,7 +26,8 @@ action() { unset __conda_setup conda activate PnetEnv + #conda activate /afs/cern.ch/work/s/skeshri/.conda/pkgs source "$( law completion )" "" } -action \ No newline at end of file +action diff --git a/law.cfg b/law.cfg index 3ed5898..7c53f50 100644 --- a/law.cfg +++ b/law.cfg @@ -1,6 +1,7 @@ [modules] ComputeRate.tasks ComputeEfficiency.tasks +Optimisation.tasks [job] job_file_dir: $RUN_PATH/jobs From 0d8c61000ac7a1d995efba1dc24e58989937a3a9 Mon Sep 17 00:00:00 2001 From: skeshri Date: Tue, 21 May 2024 22:44:59 +0200 Subject: [PATCH 5/5] updated MuTauDataset class --- HLTClass/MuTauDataset.py | 217 ++++++++++++++++++++++++++++----------- 1 file changed, 159 insertions(+), 58 deletions(-) diff --git a/HLTClass/MuTauDataset.py b/HLTClass/MuTauDataset.py index 1fd7f59..e81bb3c 100644 --- a/HLTClass/MuTauDataset.py +++ b/HLTClass/MuTauDataset.py @@ -3,16 +3,53 @@ from HLTClass.dataset import Dataset from HLTClass.dataset import get_L1Taus, get_Taus, get_Jets, get_GenTaus, hGenTau_selection, matching_Gentaus, matching_L1Taus_obj, compute_PNet_charge_prob, get_L1Muons, get_Muons, get_GenMuons, GenMuon_selection, matching_GenMuons, matching_L1Muons_obj from helpers import delta_r - +import math +import numba as nb def matching_dR_Min0p5(Obj1, Obj2, dR_matching_min = 0.5): - obj1_inpair, obj2_inpair = ak.unzip(ak.cartesian([Obj2, Obj1], nested=True)) + obj1_inpair, obj2_inpair = ak.unzip(ak.cartesian([Obj2, Obj1], nested=False)) dR_obj2_obj1 = delta_r(obj2_inpair, obj1_inpair) - mask_obj2_obj1 = (dR_obj2_obj1 > dR_matching_min) - - mask = ak.any(mask_obj2_obj1, axis=-1) + mask_obj2_obj1 = (dR_obj2_obj1 < dR_matching_min) + #mask = ak.any(mask_obj2_obj1, axis=-1) + mask = mask_obj2_obj1 return mask +#@nb.jit(nopython=True) +def phi_mpi_pi(x: float) -> float: + # okay + while (x >= 3.14159): + x -= (2 * 3.14159) + while (x < -3.14159): + x += (2 * 3.14159) + return x +#@nb.jit(nopython=True) +def deltaR(eta1: float, phi1: float, eta2: float, phi2: float) -> float: + deta = eta1 - eta2 + dphi = phi_mpi_pi(phi1 - phi2) + return float(math.sqrt(deta * deta + dphi * dphi)) + + +#@nb.jit(nopython=True) +def apply_ovrm(builder, obj1_eta, obj1_phi, obj2_eta, obj2_phi): + for iev in range(len(obj1_eta)): + builder.begin_list() + for j_eta, j_phi in zip(obj1_eta[iev], obj1_phi[iev]): + num_matches = 0 + dR = 999 + good_jet = False + for t_eta, t_phi in zip(obj2_eta[iev], obj2_phi[iev]): + # only save on last tau, so set a boolean here and fill it outside of the loop :) + dR = deltaR(j_eta, j_phi, t_eta, t_phi) + if dR > 0.5: + good_jet = True + builder.append(good_jet) + builder.end_list() + print("shape: ",len(builder.snapshot())) + ov_rm_mask = ak.sum(builder.snapshot(),axis=-1)>=1 + print("shape ovrm: ",len(ov_rm_mask)) + + return ov_rm_mask + # ------------------------------ functions for MuTau with PNet --------------------------------------------------------------- def compute_PNet_WP_MuTau(tau_pt, par): @@ -58,7 +95,7 @@ def evt_sel_MuTau(events, par, is_gen = False): Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = True) Muon_mask = Muon_selection(events) - # at least 1 L1tau/ recoJet should pass + # at least 1 L1tau & L1Mu and 1 recoJet & 1 recoMu should pass MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) L1Taus = get_L1Taus(events) @@ -73,11 +110,12 @@ def evt_sel_MuTau(events, par, is_gen = False): Muons = get_Muons(events) Selected_Muons = Muons[Muon_mask] - #overlap_mask1 = matching_dR_Min0p5(L1Muons_Muon, L1Taus_Tau) - #overlap_mask2 = matching_dR_Min0p5(Selected_Muons, Jets_Tau) + # OverlapMask (om) + om1 = apply_ovrm(ak.ArrayBuilder(), L1Muons_Muon["eta"], L1Muons_Muon["phi"], L1Taus_Tau["eta"], L1Taus_Tau["phi"]) + om2 = apply_ovrm(ak.ArrayBuilder(), Selected_Muons["eta"], Selected_Muons["phi"], Jets_Tau["eta"], Jets_Tau["phi"]) if is_gen: - # if MC data, at least 1 GenTau should also pass + # if MC data, at least 1 GenTau and 1 GenMuon should also pass GenTau_mask = hGenTau_selection(events) GenMuon_mask = GenMuon_selection(events) MuTau_evt_mask = MuTau_evt_mask & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) @@ -87,42 +125,41 @@ def evt_sel_MuTau(events, par, is_gen = False): GenTaus = get_GenTaus(events) GenTaus_Tau = GenTaus[GenTau_mask] matchingGentaus_mask = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Tau) - # at least 1 GenTau should match L1Tau/Taus + # at least 1 GenTau and 1 GenMuon should match L1Tau/Taus and L1Muon/Muon respectively evt_Taumask_matching = (ak.sum(matchingGentaus_mask, axis=-1) >= 1) GenMuons = get_GenMuons(events) GenMuons_Muon = GenMuons[GenMuon_mask] matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) - #overlap_mask3 = matching_dR_Min0p5(GenMuons_Muon,GenTaus_Tau) - #overlap_mask = overlap_mask1 & overlap_mask2 & overlap_mask3 - # at least 1 GenTau should match L1Tau/Taus + om3 = apply_ovrm(ak.ArrayBuilder(), GenMuons_Muon["eta"], GenMuons_Muon["phi"], GenTaus_Tau["eta"], GenTaus_Tau["phi"]) + overlap_mask = om1 & om2 & om3 evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) else: matchingJets_mask = matching_L1Taus_obj(L1Taus_Tau, Jets_Tau) - # at least 1 Tau should match L1Tau + # at least 1 Tau should match L1Tau evt_Taumask_matching = (ak.sum(matchingJets_mask, axis=-1) >= 1) - + + # at least 1 Muon should match L1Muon matchingMuons_mask = matching_L1Muons_obj(L1Muons_Muon, Selected_Muons) evt_Muonmask_matching = (ak.sum(matchingMuons_mask, axis=-1) >= 1) - #overlap_mask = overlap_mask1 & overlap_mask2 + overlap_mask = om1 & om2 - #MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask if is_gen: return MuTau_evt_mask, matchingGentaus_mask, matchingGenMuons_mask else: return MuTau_evt_mask, matchingJets_mask, matchingMuons_mask -# ------------------------------ functions for DoubleMediumDeepTauPFTauHPS35_L2NN_eta2p1 --------------------------------------------------- +# ------------------------------ functions for HLT_IsoMu20_eta2p1_LooseDeepTauPFTauHPS27_eta2p1_CrossL1 --------------------------------------------------- def compute_DeepTau_WP_MuTau(tau_pt): # return DeepTau WP for LooseDeepTauPFTauHPS27_eta2p1 - t1 = 0.649 - t2 = 0.441 - t3 = 0.05 - x1 = 35 + + t1 = 0.5419 + t2 = 0.4837 + t3 = 0.050 + x1 = 27 x2 = 100 x3 = 300 - Tau_WP = tau_pt*0. ones = tau_pt/tau_pt Tau_WP = ak.where((tau_pt <= ones*x1) == False, Tau_WP, ones*t1) @@ -134,7 +171,7 @@ def compute_DeepTau_WP_MuTau(tau_pt): def Tau_selection_Tau(events, apply_DeepTau_WP = True): # return mask for Tau passing selection for LooseDeepTauPFTauHPS27_eta2p1 tau_pt = events['Tau_pt'].compute() - Tau_mask = (tau_pt >= 30) & (np.abs(events['Tau_eta'].compute()) <= 2.1) + Tau_mask = (tau_pt >= 27) & (np.abs(events['Tau_eta'].compute()) <= 2.1) if apply_DeepTau_WP: Tau_mask = Tau_mask & (events['Tau_deepTauVSjet'].compute() >= compute_DeepTau_WP_MuTau(tau_pt)) return Tau_mask @@ -144,13 +181,24 @@ def evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is L1_Mu_mask = L1_Mu18er2p1_selection(events) L1_Tau_mask = L1_Tau24er2p1_OR_Tau26er2p1_selection(events) - Tau_mask = Tau_selection_Tau(events) Muon_mask = Muon_selection(events) - # at least 1 L1tau/ recoJet should pass - MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) + + print("L1_Mu_mask: ",ak.sum(ak.sum(L1_Mu_mask,axis=-1)>=1)) + print("L1_Tau_mask: ",ak.sum(ak.sum(L1_Tau_mask,axis=-1)>=1)) + print("L1_Mu_mask & L1_Tau_mask: ",ak.sum((ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1))) + print("L1_MuTau direct from branch: ",ak.sum(events["L1_Mu18er2p1_Tau24er2p1"].compute()>=1,axis=-1)) + + + # at least 1 L1tau & L1Moun and at least 1 recoTau & recoMuon should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) + + print("offline Tau_mask: ",ak.sum((ak.sum(Tau_mask, axis=-1) >= 1))) + print("offline Muon_mask: ",ak.sum((ak.sum(Muon_mask, axis=-1) >= 1))) + print("offline Muon and Tau mask: ",ak.sum((ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1))) + L1Taus = get_L1Taus(events) L1Taus_Tau = L1Taus[L1_Tau_mask] @@ -163,42 +211,51 @@ def evt_sel_HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1(events, is Muons = get_Muons(events) Selected_Muons = Muons[Muon_mask] - #overlap_mask1 = matching_dR_Min0p5(L1Muons_Muon, L1Taus_Tau) - #overlap_mask2 = matching_dR_Min0p5(Selected_Muons, Taus_Tau) - + #OverlapMask (om) + om1 = apply_ovrm(ak.ArrayBuilder(), L1Muons_Muon["eta"], L1Muons_Muon["phi"], L1Taus_Tau["eta"], L1Taus_Tau["phi"]) + om2 = apply_ovrm(ak.ArrayBuilder(), Selected_Muons["eta"], Selected_Muons["phi"], Taus_Tau["eta"], Taus_Tau["phi"]) if is_gen: - # if MC data, at least 1 GenTau should also pass + # if MC data, at least 1 GenTau and 1 GenMuon should also pass GenTau_mask = hGenTau_selection(events) GenMuon_mask = GenMuon_selection(events) MuTau_evt_mask = MuTau_evt_mask & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) + print("# of Gen Taus: ",ak.sum(GenTau_mask)) + print("# of Gen Muons: ",ak.sum(GenMuon_mask)) + print("MuTau_evt_mask & Gen Muon + Gen Tau mask: ",ak.sum(MuTau_evt_mask)) # matching if is_gen: GenTaus = get_GenTaus(events) GenTaus_Tau = GenTaus[GenTau_mask] matchingGentaus_mask = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Tau) - # at least 1 GenTau should match L1Tau/Taus + # at least 1 GenTau should match L1Tau/Taus and 1 GenMuon should match with L1Muon/Moun evt_Taumask_matching = (ak.sum(matchingGentaus_mask, axis=-1) >= 1) GenMuons = get_GenMuons(events) GenMuons_Muon = GenMuons[GenMuon_mask] matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) - #overlap_mask3 = matching_dR_Min0p5(GenMuons_Muon,GenTaus_Tau) - #overlap_mask = overlap_mask1 & overlap_mask2 & overlap_mask3 - # at least 1 GenTau should match L1Tau/Taus + om3 = apply_ovrm(ak.ArrayBuilder(), GenMuons_Muon["eta"], GenMuons_Muon["phi"], GenTaus_Tau["eta"], GenTaus_Tau["phi"]) + overlap_mask = om1 & om2 & om3 evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) + + print("L1, offline matching with Gen Muon: ",ak.sum(evt_Muonmask_matching)) + print("L1, offline matching with Gen Tau: ",ak.sum(evt_Taumask_matching)) + else: matchingTaus_mask = matching_L1Taus_obj(L1Taus_Tau, Taus_Tau) # at least 1 Tau should match L1Tau evt_Taumask_matching = (ak.sum(matchingTaus_mask, axis=-1) >= 1) + # at least 1 Muon should match L1Muon matchingMuons_mask = matching_L1Muons_obj(L1Muons_Muon, Selected_Muons) evt_Muonmask_matching = (ak.sum(matchingMuons_mask, axis=-1) >= 1) - #overlap_mask = overlap_mask1 & overlap_mask2 + overlap_mask = om1 & om2 + print("L1 and offline Muon matching: ",ak.sum(evt_Muonmask_matching)) + print("L1 and offline Tau matching: ",ak.sum(evt_Taumask_matching)) - #MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask + print("Final Muon and Tau mask: ",ak.sum(MuTau_evt_mask)) if is_gen: return MuTau_evt_mask, matchingGentaus_mask, matchingGenMuons_mask else: @@ -211,14 +268,14 @@ def L1_Mu18er2p1_selection(events): return L1_Mu18er2p1 def L1_Tau24er2p1_OR_Tau26er2p1_selection(events): - L1_Tau24er2p1_OR_Tau26er2p1 = ((events['L1Tau_pt'].compute() >=24) | (events['L1Tau_pt'].compute() >= 26)) & (events['L1Tau_eta'].compute() <= 2.131) & (events['L1Tau_eta'].compute() >= -2.131) + L1_Tau24er2p1_OR_Tau26er2p1 = ((events['L1Tau_pt'].compute() >=24)) & (events['L1Tau_eta'].compute() <= 2.131) & (events['L1Tau_eta'].compute() >= -2.131) return L1_Tau24er2p1_OR_Tau26er2p1 def Denominator_Selection_MuTau(GenLepton): # return mask for event passing minimal Gen requirements for MuTau HLT mask_GenTau = (GenLepton['kind'] == 5) mask_GenMuon = (GenLepton['kind'] == 4) - ev_mask = (ak.sum(mask_GenTau, axis=-1) == 1) & (ak.sum(mask_GenMuon, axis=-1) == 1) # 1 Gen tau and one Gen Muon should pass this requirements + ev_mask = (ak.sum(mask_GenTau, axis=-1) >= 1) & (ak.sum(mask_GenMuon, axis=-1) >= 1) # 1 Gen tau and one Gen Muon should pass this requirements return ev_mask class MuTauDataset(Dataset): @@ -227,7 +284,7 @@ def __init__(self, fileName): # ------------------------------ functions to Compute Rate --------------------------------------------------------------- - def get_Nnum_Nden_MuTauDeepNet(self): + def get_Nnum_Nden_HLT_MuTauDeepNet(self): print(f'Computing rate for HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1:') #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator events = self.get_events() @@ -252,7 +309,7 @@ def get_Nnum_Nden_MuTauPNet(self, par): print('') return N_den, N_num - def get_Nnum_Nden_MuTauDeepNet_nPV(self, nPV): + def get_Nnum_Nden_HLT_MuTauDeepNet_nPV(self, nPV): print(f'Computing rate for HLT_IsoMu20_eta2p1_PNetTauhPFJetPt30_Loose_eta2p3_CrossL1 for nPV = {nPV}:') #load all events in the file that belong to the run and lumiSections_range, save the number of events in Denominator events = self.get_events() @@ -363,8 +420,15 @@ def ComputeEffAlgo_MuTauPNet(self, par): Muon_mask = Muon_selection(events) - # at least 1 L1tau/ recoJet should pass - MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + print("len GenTauMask: ", ak.sum((ak.sum(GenTau_mask, axis=-1) >= 1) )) + print("len GenMuonMask: ", ak.sum((ak.sum(GenMuon_mask, axis=-1) >= 1) )) + print("len GenMask: ", ak.sum((ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1))) + print("len GenMask & tau mask: ", ak.sum((ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1))) + print("len L1 & Offline Tau and Muon mask: ", ak.sum((ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) )) + + + # at least 1 L1tau & L1Muon and at least 1 recoJet & recoMuon and at least 1 GenTau and GenMuon should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) L1Taus = get_L1Taus(events) L1Taus_Tau = L1Taus[L1_Tau_mask] @@ -390,23 +454,39 @@ def ComputeEffAlgo_MuTauPNet(self, par): matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) + print("len evt_Taumask_matching: ", ak.sum(evt_Taumask_matching )) + print("len evt_Muonmask_matching: ", ak.sum(evt_Muonmask_matching )) + print("len Tau and Muon mask matching: ",ak.sum(evt_Taumask_matching & evt_Muonmask_matching)) + + # OverlapMask (om) + om1 = apply_ovrm(ak.ArrayBuilder(), L1Muons_Muon["eta"], L1Muons_Muon["phi"], L1Taus_Tau["eta"], L1Taus_Tau["phi"]) + om2 = apply_ovrm(ak.ArrayBuilder(), Selected_Muons["eta"], Selected_Muons["phi"], Jets_Tau["eta"], Jets_Tau["phi"]) + om3 = apply_ovrm(ak.ArrayBuilder(), GenMuons_Muon["eta"], GenMuons_Muon["phi"], GenTaus_Sel["eta"], GenTaus_Sel["phi"]) + overlap_mask = om1 & om2 & om3 + + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + print("=========== the PNet Disabled ================") + print("len evt_Taumask_matching: ",len(events[evt_Taumask_matching])) + print("len evt_Muonmask_matching: ",len(events[evt_Muonmask_matching])) + print("len overlap_mask: ",len(events[overlap_mask])) + N_den = len(events[MuTau_evt_mask]) print(f"Number of events in denominator: {N_den}") - # Selection of L1/Gen and Jets objects with PNET WP Tau_mask = Jet_selection_Tau(events, par, apply_PNET_WP = True) - # at least 1 L1tau/ Jet/ GenTau should pass - MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) #matching Jets_Tau = Jets[Tau_mask] matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Jets_Tau, GenTaus_Sel) evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask + print("=========== After enabling the PNet ================") + print("len evt_Taumask_matching: ",len(events[evt_Taumask_matching])) + print("len evt_Muonmask_matching: ",len(events[evt_Muonmask_matching])) + print("len overlap_mask: ",len(events[overlap_mask])) N_num = len(events[MuTau_evt_mask]) print(f"Number of events in numerator: {N_num}") @@ -427,8 +507,14 @@ def ComputeEffAlgo_MuTauDeepNet(self): Muon_mask = Muon_selection(events) - # at least 1 L1tau/ recoJet should pass - MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) + print("len GenTauMask: ", ak.sum((ak.sum(GenTau_mask, axis=-1) >= 1) )) + print("len GenMuonMask: ", ak.sum((ak.sum(GenMuon_mask, axis=-1) >= 1) )) + print("len GenMask: ", ak.sum((ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1))) + print("len GenMask & tau mask: ", ak.sum((ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1))) + print("len L1 & Offline Tau and Muon mask: ", ak.sum((ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) )) + + # at least 1 L1tau & L1Muon and at least 1 recoJet & recoMuon and at least 1 GenTau and GenMuon should pass + MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) & (ak.sum(GenMuon_mask, axis=-1) >= 1) L1Taus = get_L1Taus(events) L1Taus_Tau = L1Taus[L1_Tau_mask] @@ -454,25 +540,40 @@ def ComputeEffAlgo_MuTauDeepNet(self): matchingGenMuons_mask = matching_GenMuons(L1Muons_Muon, Selected_Muons, GenMuons_Muon) evt_Muonmask_matching = (ak.sum(matchingGenMuons_mask, axis=-1) >= 1) - - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + + print("len evt_Taumask_matching: ", ak.sum(evt_Taumask_matching )) + print("len evt_Muonmask_matching: ", ak.sum(evt_Muonmask_matching )) + print("len Tau and Muon mask matching: ",ak.sum(evt_Taumask_matching & evt_Muonmask_matching)) + + # OverlapMask (om) + om1 = apply_ovrm(ak.ArrayBuilder(), L1Muons_Muon["eta"], L1Muons_Muon["phi"], L1Taus_Tau["eta"], L1Taus_Tau["phi"]) + om2 = apply_ovrm(ak.ArrayBuilder(), Selected_Muons["eta"], Selected_Muons["phi"], Taus_Tau["eta"], Taus_Tau["phi"]) + om3 = apply_ovrm(ak.ArrayBuilder(), GenMuons_Muon["eta"], GenMuons_Muon["phi"], GenTaus_Sel["eta"], GenTaus_Sel["phi"]) + overlap_mask = om1 & om2 & om3 + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask + N_den = len(events[MuTau_evt_mask]) + print("=========== the deeptau Disabled ================") + print("len evt_Taumask_matching: ",len(events[evt_Taumask_matching])) + print("len evt_Muonmask_matching: ",len(events[evt_Muonmask_matching])) + print("len overlap_mask: ",len(events[overlap_mask])) print(f"Number of events in denominator: {N_den}") - # Selection of L1/Gen and Jets objects with PNET WP Tau_mask = Tau_selection_Tau(events, apply_DeepTau_WP = True) - # at least 1 L1tau/ Jet/ GenTau should pass - MuTau_evt_mask = (ak.sum(L1_Mu_mask, axis=-1) >= 1) & (ak.sum(L1_Tau_mask, axis=-1) >= 1) & (ak.sum(Tau_mask, axis=-1) >= 1) & (ak.sum(Muon_mask, axis=-1) >= 1) & (ak.sum(GenTau_mask, axis=-1) >= 1) #matching Taus_Tau = Taus[Tau_mask] matching_GenTaus_mask_Tau = matching_Gentaus(L1Taus_Tau, Taus_Tau, GenTaus_Sel) evt_Taumask_matching = (ak.sum(matching_GenTaus_mask_Tau, axis=-1) >= 1) - MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching + MuTau_evt_mask = MuTau_evt_mask & evt_Taumask_matching & evt_Muonmask_matching & overlap_mask N_num = len(events[MuTau_evt_mask]) + print("=========== After enabling the deeptau ================") + print("len evt_Taumask_matching: ",len(events[evt_Taumask_matching])) + print("len evt_Muonmask_matching: ",len(events[evt_Muonmask_matching])) + print("len overlap_mask: ",len(events[overlap_mask])) print(f"Number of events in numerator: {N_num}") print('') return N_den, N_num