From 7b0f67996df811b462279e2d2de25d3c6cd3fa3f Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Tue, 8 Apr 2025 11:38:38 -0700 Subject: [PATCH 1/7] Commit the stuff from the bbzz checks for ref --- analysis/wwz/get_bbzz_yields.py | 122 +++++++++++++++++++++++ analysis/wwz/wwz4l.py | 98 +++++++++++++++--- ewkcoffea/modules/selection_wwz.py | 43 ++++++++ input_samples/cfgs/wwz_analysis/bbzz.cfg | 116 +++++++++++++++++++++ 4 files changed, 367 insertions(+), 12 deletions(-) create mode 100644 analysis/wwz/get_bbzz_yields.py create mode 100644 input_samples/cfgs/wwz_analysis/bbzz.cfg diff --git a/analysis/wwz/get_bbzz_yields.py b/analysis/wwz/get_bbzz_yields.py new file mode 100644 index 00000000..8888950b --- /dev/null +++ b/analysis/wwz/get_bbzz_yields.py @@ -0,0 +1,122 @@ +import pickle +import hist +import gzip + +import ewkcoffea.modules.sample_groupings as sg + + +################################## +# Make sample grouping dict +GROUPING_DICT_R2= { + #"ZZ" : ["ZZTo4l", "ggToZZTo2e2mu", "ggToZZTo2e2tau", "ggToZZTo2mu2tau", "ggToZZTo4e", "ggToZZTo4mu", "ggToZZTo4tau"], + #"ttZ" : [ "TTZToLL_M_1to10", "TTZToLLNuNu_M_10", ], + "higgs" : [ "ggHToZZ4L", ], + #"HH" : [ "HHTobbVV", ], + #"WZ" : [ "WZTo3LNu", ], +} + +GROUPING_DICT_R3= { + "ZZ" : ["ZZTo4l", "ggToZZTo2e2mu", "ggToZZTo2e2tau", "ggToZZTo2mu2tau", "ggToZZTo4e", "ggToZZTo4mu", "ggToZZTo4tau"], + "ttZ" : [ "TTZToLL_M_4to50", "TTZToLL_M_50", ], + "higgs" : [ "ggHToZZ4L", ], +} + +def attach_years(year): + out_dict = {} + if year == "run2": + years = ["UL16APV","UL16","UL17","UL18"] + #years = ["UL18"] + sample_dict_base = GROUPING_DICT_R2 + elif year == "run3": + years = ["2022","2022EE","2023","2023BPix"] + sample_dict_base = GROUPING_DICT_R3 + else: + raise Exception("Unknown year") + for proc_group in sample_dict_base.keys(): + out_dict[proc_group] = [] + for proc_base_name in sample_dict_base[proc_group]: + for year_str in years: + out_dict[proc_group].append(f"{year_str}_{proc_base_name}") + return out_dict + + +################################## + + +################################## +# Get the list of categories on the sparese axis +def get_axis_cats(histo,axis_name): + process_list = [x for x in histo.axes[axis_name]] + return process_list +# Regroup categories (e.g. processes) +def group(h, oldname, newname, grouping): + + # Build up a grouping dict that drops any proc that is not in our h + grouping_slim = {} + proc_lst = get_axis_cats(h,oldname) + for grouping_name in grouping.keys(): + for proc in grouping[grouping_name]: + if proc in proc_lst: + if grouping_name not in grouping_slim: + grouping_slim[grouping_name] = [] + grouping_slim[grouping_name].append(proc) + #else: + # print(f"WARNING: process {proc} not in this hist") + + # From Nick: https://github.com/CoffeaTeam/coffea/discussions/705#discussioncomment-4604211 + hnew = hist.Hist( + hist.axis.StrCategory(grouping_slim, name=newname), + *(ax for ax in h.axes if ax.name != oldname), + storage=h.storage_type(), + ) + for i, indices in enumerate(grouping_slim.values()): + hnew.view(flow=True)[i] = h[{oldname: indices}][{oldname: sum}].view(flow=True) + + return hnew +################################## + +def main(): + + year = "run2" + pkl_file_path = "histos/tmp.pkl.gz" + #pkl_file_path = "histos/r2_zz_h_ttz_noSyst_1.pkl.gz" # Standard selection + #pkl_file_path = "histos/r2_zz_h_ttz_noSyst_lowpt.pkl.gz" # Lower lep pt + #year = "run3" + #pkl_file_path = "histos/r3_zz_h_ttz_noSyst.pkl.gz" # Standard selection + + # Make grouping dict + GOUPING_DICT = attach_years(year) + + histo_dict = pickle.load(gzip.open(pkl_file_path)) + + cat_lst = [ + "all_events", + "bbzz", + "bbzz_zzcand", + "bbzz_zzcand_bb", + "bbzz_zzcand_bb_m4l", + + "bbzz_zzcand_bb_m4l_4m", + "bbzz_zzcand_bb_m4l_4e", + "bbzz_zzcand_bb_m4l_2e2m", + "sr_4l_bdt_of_trn", + ] + + histo = histo_dict["njets"][{"systematic":"nominal"}] + histo_grouped = group(histo,"process","process_grp",GOUPING_DICT) + + for proc in GOUPING_DICT: + + print("\n-------",proc,"-------") + + for cat in cat_lst: + histo = histo_grouped[{"category":cat,"process_grp":proc}] + print("\n",cat) + val = sum(histo.values(flow=True)) + err = (sum(histo.variances(flow=True)))**0.5 + perr = round(100*err/val,1) + print(f"{val} +- {err}, {perr}%") + + + +main() diff --git a/analysis/wwz/wwz4l.py b/analysis/wwz/wwz4l.py index 9448092c..fd60e43c 100644 --- a/analysis/wwz/wwz4l.py +++ b/analysis/wwz/wwz4l.py @@ -67,17 +67,24 @@ def __init__(self, samples, wc_names_lst=[], hist_lst=None, ecut_threshold=None, "scalarptsum_jet" : axis.Regular(180, 0, 700, name="scalarptsum_jet", label="H_T"), "mll_01": axis.Regular(180, 0, 200, name="mll_01", label="mll_l0_l1"), "mllll": axis.Regular(180, 0, 1000, name="mllll", label="mllll"), - "l0pt" : axis.Regular(180, 0, 500, name="l0pt", label="l0pt"), - "j0pt" : axis.Regular(180, 0, 500, name="j0pt", label="j0 pt"), + "l0pt" : axis.Regular(180, 0, 300, name="l0pt", label="l0pt"), + "l1pt" : axis.Regular(180, 0, 200, name="l1pt", label="l1pt"), + "l2pt" : axis.Regular(180, 0, 100, name="l2pt", label="l2pt"), + "l3pt" : axis.Regular(180, 0, 100, name="l3pt", label="l3pt"), + "j0pt" : axis.Regular(180, 0, 300, name="j0pt", label="j0 pt"), "j0eta" : axis.Regular(180, -3, 3, name="j0eta", label="j0 eta"), "j0phi" : axis.Regular(180, -3.1416, 3.1416, name="j0phi", label="j0 phi"), + "j1pt" : axis.Regular(180, 0, 300, name="j1pt", label="j1 pt"), + "j1eta" : axis.Regular(180, -3, 3, name="j1eta", label="j1 eta"), + "j1phi" : axis.Regular(180, -3.1416, 3.1416, name="j1phi", label="j1 phi"), + "mjj" : axis.Regular(180, 0, 300, name="mjj", label="mjj"), "abs_pdgid_sum" : axis.Regular(15, 40, 55, name="abs_pdgid_sum", label="Sum of abs pdgId for all 4 lep"), - "w_lep0_pt" : axis.Regular(180, 0, 300, name="w_lep0_pt", label="Leading W lep pt"), - "w_lep1_pt" : axis.Regular(180, 0, 300, name="w_lep1_pt", label="Subleading W lep pt"), - "z_lep0_pt" : axis.Regular(180, 0, 300, name="z_lep0_pt", label="Leading Z lep pt"), - "z_lep1_pt" : axis.Regular(180, 0, 300, name="z_lep1_pt", label="Subleading Z lep pt"), + "w_lep0_pt" : axis.Regular(180, 0, 100, name="w_lep0_pt", label="Leading W lep pt"), + "w_lep1_pt" : axis.Regular(180, 0, 100, name="w_lep1_pt", label="Subleading W lep pt"), + "z_lep0_pt" : axis.Regular(180, 0, 200, name="z_lep0_pt", label="Leading Z lep pt"), + "z_lep1_pt" : axis.Regular(180, 0, 200, name="z_lep1_pt", label="Subleading Z lep pt"), "w_lep0_eta" : axis.Regular(180, -3, 3, name="w_lep0_eta", label="Leading W lep eta"), "w_lep1_eta" : axis.Regular(180, -3, 3, name="w_lep1_eta", label="Subleading W lep eta"), "z_lep0_eta" : axis.Regular(180, -3, 3, name="z_lep0_eta", label="Leading Z lep eta"), @@ -388,6 +395,7 @@ def process(self, events): # Put njets and l_fo_conept_sorted into events and get 4l event selection mask events["l_wwz_t"] = l_wwz_t es_ec.add4lmask_wwz(events, year, isData, histAxisName, is2022,is2023) + es_ec.add4lmask_bbzz(events, year, isData, histAxisName, is2022,is2023) ######### Normalization and weights ########### @@ -515,8 +523,13 @@ def process(self, events): # Count jets njets = ak.num(goodJets) ht = ak.sum(goodJets.pt,axis=-1) - j0 = goodJets[ak.argmax(goodJets.pt,axis=-1,keepdims=True)] + #j0 = goodJets[ak.argmax(goodJets.pt,axis=-1,keepdims=True)] + goodJets_ptordered = goodJets[ak.argsort(goodJets.pt,axis=-1,ascending=False)] + goodJets_ptordered_padded = ak.pad_none(goodJets_ptordered, 2) + j0 = goodJets_ptordered_padded[:,0] + j1 = goodJets_ptordered_padded[:,1] + mjj = (j0+j1).mass # Loose DeepJet WP btagger = "btag" # For deep flavor WPs @@ -694,9 +707,12 @@ def process(self, events): abs_pdgid_sum = (abs(l0.pdgId) + abs(l1.pdgId) + abs(l2.pdgId) + abs(l3.pdgId)) l0pt = l0.pt - j0pt = ak.flatten(j0.pt) # Flatten to go from [[j0pt],[j0pt],...] -> [j0pt,j0pt,...] - j0eta = ak.flatten(j0.eta) # Flatten to go from [[j0pt],[j0pt],...] -> [j0pt,j0pt,...] - j0phi = ak.flatten(j0.phi) # Flatten to go from [[j0pt],[j0pt],...] -> [j0pt,j0pt,...] + j0pt = j0.pt + j0eta = j0.eta + j0phi = j0.phi + j1pt = j1.pt + j1eta = j1.eta + j1phi = j1.phi mll_01 = (l0+l1).mass mllll = (l0+l1+l2+l3).mass scalarptsum_lep = l0.pt + l1.pt + l2.pt + l3.pt @@ -759,9 +775,16 @@ def process(self, events): "mll_01" : mll_01, "mllll" : mllll, "l0pt" : l0pt, + "l1pt" : l1.pt, + "l2pt" : l2.pt, + "l3pt" : l3.pt, "j0pt" : j0pt, "j0eta" : j0eta, "j0phi" : j0phi, + "j1pt" : j1pt, + "j1eta" : j1eta, + "j1phi" : j1phi, + "mjj" : mjj, "abs_pdgid_sum": abs_pdgid_sum, @@ -818,8 +841,8 @@ def process(self, events): "mll_min_afos" : mll_min_afos, "mll_min_sfos" : mll_min_sfos, - "mlb_min" : mlb_min, - "mlb_max" : mlb_max, + #"mlb_min" : mlb_min, + #"mlb_max" : mlb_max, } # Include the genPartFlav, though this is only defined for MC, so just fill with 1 if data @@ -829,6 +852,17 @@ def process(self, events): dense_variables_dict["z_lep1_genPartFlav"] = z_lep1.genPartFlav if not isData else events.nom + ############### bbzz selection ############### + z1_candidate_sfos = w_lep0.pdgId == -(w_lep1.pdgId) + z1_candidate_good_pt = ((w_lep0+w_lep1).mass > 12) & ((w_lep0+w_lep1).mass<120) + + has_z0_candidate = (z_lep0+z_lep1).mass > 40 + has_z1_candidate = z1_candidate_sfos & z1_candidate_good_pt + + has_zz_candidate = has_z0_candidate & has_z1_candidate + #################################### + + ######### Evaluate the BDTs (get WWZ, ZH, and WZ scores for SF and OF) ######### # Get the list of variables for the BDTs (and fill None with -9999 to not cause problems), and eval @@ -1014,6 +1048,36 @@ def process(self, events): selections.add("sr_4l_sf_incl", (veto_map_mask & pass_trg & events.is4lWWZ & bmask_exactly0loose & events.wwz_presel_sf & w_candidates_mll_far_from_z & (met.pt >= 65.0))) # Inclusive over SF sr (only applying cuts that are applied to all SF SRs), just use for visualization selections.add("sr_4l_of_incl", (veto_map_mask & pass_trg & events.is4lWWZ & bmask_exactly0loose & events.wwz_presel_of)) # Inclusive over OF sr (only applying cuts that are applied to all OF SRs), just use for visualization + ########## For bbzz SRs ########## + z0_ee = ((abs(z_lep0.pdgId) == 11) & (abs(z_lep1.pdgId) == 11)) + z1_ee = ((abs(w_lep0.pdgId) == 11) & (abs(w_lep1.pdgId) == 11)) + z0_mm = ((abs(z_lep0.pdgId) == 13) & (abs(z_lep1.pdgId) == 13)) + z1_mm = ((abs(w_lep0.pdgId) == 13) & (abs(w_lep1.pdgId) == 13)) + + mllll_bbzz = (mllll > 115) & (mllll < 135) + bbzz_4e = z0_ee & z1_ee + bbzz_4m = z0_mm & z1_mm + bbzz_2e2m = (z0_ee & z1_mm) | (z0_mm & z1_ee) + selections.add("bbzz", (veto_map_mask & pass_trg & events.is4lbbzz)) + selections.add("bbzz_zzcand", (veto_map_mask & pass_trg & events.is4lbbzz & has_zz_candidate)) + selections.add("bbzz_zzcand_bb", (veto_map_mask & pass_trg & events.is4lbbzz & has_zz_candidate & (njets>=2))) + selections.add("bbzz_zzcand_bb_m4l", (veto_map_mask & pass_trg & events.is4lbbzz & has_zz_candidate & (njets>=2) & mllll_bbzz)) + selections.add("bbzz_zzcand_bb_m4l_4e", (veto_map_mask & pass_trg & events.is4lbbzz & has_zz_candidate & (njets>=2) & mllll_bbzz & bbzz_4e)) + selections.add("bbzz_zzcand_bb_m4l_4m", (veto_map_mask & pass_trg & events.is4lbbzz & has_zz_candidate & (njets>=2) & mllll_bbzz & bbzz_4m)) + selections.add("bbzz_zzcand_bb_m4l_2e2m", (veto_map_mask & pass_trg & events.is4lbbzz & has_zz_candidate & (njets>=2) & mllll_bbzz & bbzz_2e2m)) + bbzz_cat_lst = [ + "all_events", + "sr_4l_bdt_of_trn", + "bbzz", + "bbzz_zzcand", + "bbzz_zzcand_bb", + "bbzz_zzcand_bb_m4l", + "bbzz_zzcand_bb_m4l_4e", + "bbzz_zzcand_bb_m4l_4m", + "bbzz_zzcand_bb_m4l_2e2m", + ] + ################################## + # For BDT SRs sr_4l_bdt_sf_presel = (veto_map_mask & pass_trg & events.is4lWWZ & bmask_exactly0loose & events.wwz_presel_sf & w_candidates_mll_far_from_z) @@ -1113,6 +1177,9 @@ def process(self, events): cat_dict["lep_chan_lst"] = cat_dict["lep_chan_lst"] + bdt_sr_names + bdt_misc_names + # TMP: Just hard code list of channels + cat_dict["lep_chan_lst"] = bbzz_cat_lst + ######### Fill histos ######### # List the hists that are only defined for some categories @@ -1124,7 +1191,14 @@ def process(self, events): "j0pt" : ["all_events", "4l_presel", "sr_4l_sf", "sr_4l_of", "sr_4l_sf_incl", "sr_4l_of_incl", "cr_4l_sf", "cr_4l_sf_higgs"] + analysis_cats + bdt_misc_names, "j0eta" : ["all_events", "4l_presel", "sr_4l_sf", "sr_4l_of", "sr_4l_sf_incl", "sr_4l_of_incl", "cr_4l_sf", "cr_4l_sf_higgs"] + analysis_cats + bdt_misc_names, "j0phi" : ["all_events", "4l_presel", "sr_4l_sf", "sr_4l_of", "sr_4l_sf_incl", "sr_4l_of_incl", "cr_4l_sf", "cr_4l_sf_higgs"] + analysis_cats + bdt_misc_names, + "j1pt" : ["all_events", "4l_presel", "sr_4l_sf", "sr_4l_of", "sr_4l_sf_incl", "sr_4l_of_incl", "cr_4l_sf", "cr_4l_sf_higgs"] + analysis_cats + bdt_misc_names, + "j1eta" : ["all_events", "4l_presel", "sr_4l_sf", "sr_4l_of", "sr_4l_sf_incl", "sr_4l_of_incl", "cr_4l_sf", "cr_4l_sf_higgs"] + analysis_cats + bdt_misc_names, + "j1phi" : ["all_events", "4l_presel", "sr_4l_sf", "sr_4l_of", "sr_4l_sf_incl", "sr_4l_of_incl", "cr_4l_sf", "cr_4l_sf_higgs"] + analysis_cats + bdt_misc_names, + "mjj" : ["all_events", "4l_presel", "sr_4l_sf", "sr_4l_of", "sr_4l_sf_incl", "sr_4l_of_incl", "cr_4l_sf", "cr_4l_sf_higgs"] + analysis_cats + bdt_misc_names, "l0pt" : ["all_events"], + "l1pt" : ["all_events"], + "l2pt" : ["all_events"], + "l3pt" : ["all_events"], "mll_01" : ["all_events"], "mllll" : ["all_events"], "scalarptsum_lep" : ["all_events"], diff --git a/ewkcoffea/modules/selection_wwz.py b/ewkcoffea/modules/selection_wwz.py index 269d6b2e..bce0dc61 100644 --- a/ewkcoffea/modules/selection_wwz.py +++ b/ewkcoffea/modules/selection_wwz.py @@ -366,6 +366,49 @@ def add4lmask_wwz(events, year, isData, sample_name,is2022,is2023): events['is4lWWZ'] = ak.fill_none(mask,False) +# For the HH bbZZ 4l analysis checks +def add4lmask_bbzz(events, year, isData, sample_name,is2022,is2023): + + # Leptons and padded leptons + leps = events.l_wwz_t + leps_padded = ak.pad_none(leps,4) + + # Filters + filter_flags = events.Flag + if (is2022 or is2023): + filters = filter_flags.goodVertices & filter_flags.globalSuperTightHalo2016Filter & filter_flags.EcalDeadCellTriggerPrimitiveFilter & filter_flags.BadPFMuonFilter & filter_flags.ecalBadCalibFilter & filter_flags.BadPFMuonDzFilter & filter_flags.hfNoisyHitsFilter & filter_flags.eeBadScFilter + elif year in ["2016","2016APV"]: + filters = filter_flags.goodVertices & filter_flags.globalSuperTightHalo2016Filter & filter_flags.HBHENoiseFilter & filter_flags.HBHENoiseIsoFilter & filter_flags.EcalDeadCellTriggerPrimitiveFilter & filter_flags.BadPFMuonFilter & filter_flags.BadPFMuonDzFilter & filter_flags.eeBadScFilter + else: + filters = filter_flags.goodVertices & filter_flags.globalSuperTightHalo2016Filter & filter_flags.HBHENoiseFilter & filter_flags.HBHENoiseIsoFilter & filter_flags.EcalDeadCellTriggerPrimitiveFilter & filter_flags.BadPFMuonFilter & filter_flags.BadPFMuonDzFilter & filter_flags.eeBadScFilter & filter_flags.ecalBadCalibFilter + + # Lep multiplicity + nlep_4 = (ak.num(leps) == 4) + + # Remove low mass resonances + cleanup = (events.min_mll_afos > 12) + + #mask = filters & nlep_4 & on_z & cleanup + mask = filters & nlep_4 & cleanup + + # Do gen cleanups + if sample_name in get_ec_param("vh_list"): + genparts = events.GenPart + is_zh = (abs(genparts[:,2].pdgId) == 23) # 3rd genparticle should be v for these samples + is_w_from_h = ((abs(genparts.pdgId)==24) & (abs(genparts.distinctParent.pdgId) == 25)) + gen_mask = ~(is_zh & ak.any(is_w_from_h,axis=-1)) + mask = mask & gen_mask + + # SFs: + events['sf_4l_muon'] = leps_padded[:,0].sf_nom_muon*leps_padded[:,1].sf_nom_muon*leps_padded[:,2].sf_nom_muon*leps_padded[:,3].sf_nom_muon + events['sf_4l_elec'] = leps_padded[:,0].sf_nom_elec*leps_padded[:,1].sf_nom_elec*leps_padded[:,2].sf_nom_elec*leps_padded[:,3].sf_nom_elec + events['sf_4l_hi_muon'] = leps_padded[:,0].sf_hi_muon*leps_padded[:,1].sf_hi_muon*leps_padded[:,2].sf_hi_muon*leps_padded[:,3].sf_hi_muon + events['sf_4l_hi_elec'] = leps_padded[:,0].sf_hi_elec*leps_padded[:,1].sf_hi_elec*leps_padded[:,2].sf_hi_elec*leps_padded[:,3].sf_hi_elec + events['sf_4l_lo_muon'] = leps_padded[:,0].sf_lo_muon*leps_padded[:,1].sf_lo_muon*leps_padded[:,2].sf_lo_muon*leps_padded[:,3].sf_lo_muon + events['sf_4l_lo_elec'] = leps_padded[:,0].sf_lo_elec*leps_padded[:,1].sf_lo_elec*leps_padded[:,2].sf_lo_elec*leps_padded[:,3].sf_lo_elec + + events['is4lbbzz'] = ak.fill_none(mask,False) + # Takes as input the lep collection # Finds SFOS pair that is closest to the Z peak diff --git a/input_samples/cfgs/wwz_analysis/bbzz.cfg b/input_samples/cfgs/wwz_analysis/bbzz.cfg new file mode 100644 index 00000000..5c6b81e8 --- /dev/null +++ b/input_samples/cfgs/wwz_analysis/bbzz.cfg @@ -0,0 +1,116 @@ +# Jsons for 4l skim + +# Prefix for UCSD UAF +/ceph/cms/ + +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggHToZZ4L_unskimmed.json + +# Unskimmed +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_HHTobbVV.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_HHTobbVV.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_HHTobbVV.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_HHTobbVV.json + +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_WZTo3LNu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_WZTo3LNu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_WZTo3LNu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_WZTo3LNu.json + +################### Run 2 ################### + +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_TTZToLL_M_1to10.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_TTZToLLNuNu_M_10.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_TTZToLL_M_1to10.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_TTZToLLNuNu_M_10.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_TTZToLL_M_1to10.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_TTZToLLNuNu_M_10.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_TTZToLL_M_1to10.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_TTZToLLNuNu_M_10.json + +../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ggHToZZ4L.json +../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ggHToZZ4L.json +../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ggHToZZ4L.json +../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggHToZZ4L.json + +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16APV_ggToZZTo4tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL16_ggToZZTo4tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL17_ggToZZTo4tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/UL18_ggToZZTo4tau.json + + +################### Run 3 ################### + +## 2022 and 2023 ttZ +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_TTZToLL_M_4to50.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_TTZToLL_M_50.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_TTZToLL_M_4to50.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_TTZToLL_M_50.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_TTZToLL_M_4to50.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_TTZToLL_M_50.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_TTZToLL_M_4to50.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_TTZToLL_M_50.json +# +## 2022 ggH and ZZ +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ggHToZZ4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022_ggToZZTo4tau.json +## 2022EE ggH and ZZ +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ggHToZZ4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2022EE_ggToZZTo4tau.json +# +## 2023 ggH and ZZ +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ggHToZZ4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023_ggToZZTo4tau.json +# +## 2023Bpix ggH and ZZ +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ggHToZZ4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ZZTo4L.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ggToZZTo2e2mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ggToZZTo2e2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ggToZZTo2mu2tau.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ggToZZTo4e.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ggToZZTo4mu.json +#../../input_samples/sample_jsons/wwz_analysis_4l_skims_v1/bkg_samples/2023BPix_ggToZZTo4tau.json + + From 61236f2221492c1a5ddef85cb94a3e10714e669c Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Thu, 10 Apr 2025 12:28:10 -0700 Subject: [PATCH 2/7] In progress bdt for vvv --- analysis/bdt_training/binary_ingress.py | 356 ++++++++++++++++++++++++ analysis/bdt_training/binary_train.py | 169 +++++++++++ 2 files changed, 525 insertions(+) create mode 100644 analysis/bdt_training/binary_ingress.py create mode 100644 analysis/bdt_training/binary_train.py diff --git a/analysis/bdt_training/binary_ingress.py b/analysis/bdt_training/binary_ingress.py new file mode 100644 index 00000000..b5ec8239 --- /dev/null +++ b/analysis/bdt_training/binary_ingress.py @@ -0,0 +1,356 @@ +import argparse +import numpy as np +import pickle +import gzip +import cloudpickle +import ewkcoffea.modules.sample_groupings as sg +from ewkcoffea.modules.paths import ewkcoffea_path as ewkcoffea_path +from topcoffea.modules.get_param_from_jsons import GetParam + +_VAR_LST_TMP = [ + "njets", + "j0pt", + "nleps", +] + +VAR_LST_TMP = [ + #"absdphi_4l_met", + #"absdphi_min_afas", + #"absdphi_min_afos", + #"absdphi_min_sfos", + #"absdphi_wl0_met", + #"absdphi_wl0_wl1", + #"absdphi_wl1_met", + #"absdphi_wleps_met", + #"absdphi_zl0_zl1", + #"absdphi_zleps_met", + "absdphi_z_ww", + #"abs_pdgid_sum", + #"cos_helicity_x", + #"dr_wl0_j_min", + #"dr_wl0_wl1", + #"dr_wl1_j_min", + #"dr_wleps_zleps", + #"dr_zl0_zl1", + #"j0eta", + #"j0phi", + "j0pt", + #"j1eta", + #"j1phi", + "j1pt", + "l0pt", + "l1pt", + "l2pt", + "l3pt", + #"metphi", + #"met", + "mjj_j0j1", + "mjj_nearz", + #"mll_01", + "mllll", + #"mll_min_afas", + #"mll_min_afos", + #"mll_min_sfos", + #"mll_wl0_wl1", + #"mll_zl0_zl1", + #"mt2", + "mt_4l_met", + #"mt_wl0_met", + #"mt_wl1_met", + #"mt_wleps_met", + "nbtagsl", + "njets", + #"nleps", + #"ptl4", + "pt_wl0_wl1", + "pt_zl0_zl1", + "scalarptsum_jet", + #"scalarptsum_lepmetjet", + #"scalarptsum_lepmet", + #"scalarptsum_lep", + #"w_lep0_eta", + #"w_lep0_genPartFlav", + #"w_lep0_phi", + #"w_lep0_pt", + #"w_lep1_eta", + #"w_lep1_genPartFlav", + #"w_lep1_phi", + #"w_lep1_pt", + #"z_lep0_eta", + #"z_lep0_genPartFlav", + #"z_lep0_phi", + #"z_lep0_pt", + #"z_lep1_eta", + #"z_lep1_genPartFlav", + #"z_lep1_phi", + #"z_lep1_pt", +] + +def prepare_bdt_training_data(pklfilepath): + + # Open up the pickle file + d = pickle.load(gzip.open(pklfilepath)) + + # Dictionary where we will store our master data + myd = {} + + list_output_names = [ + "list_bdt_base_proc", + "list_bdt_base_wgt", + "list_bdt_base_evt", + ] + + get_ec_param = GetParam(ewkcoffea_path("params/params.json")) + + # Parse the pickle file and store to myd as numpy arrays + for list_output_name in list_output_names: + myd[list_output_name] = np.array(d[list_output_name]) + for var in VAR_LST_TMP: + myd["base_bdt_" + var] = np.array(d["list_base_bdt_" + var]) + + # Sample groupings to split the data into signals, and backgrounds + sample_dict_mc = sg.create_mc_sample_dict("run2") + + # List of categories that will use BDT to label + #categories = ["WWZ", "ZH", "Bkg"] + #categories = ["ZZZ", "WZZ", "Bkg"] + categories = ["Sig", "Bkg"] + + # Sample splittings + sample_splittings = ["train", "test", "all"] # later we can add validate if needed + + # Dictionary where we will store our training and testing data + bdt_data = {} + + # Creating data structure in bdt_data + # axis 0: categories + # axis 1: split + # axis 2: bdt variables and weight + for cat in categories: + bdt_data[cat] = {} + for split in sample_splittings: + bdt_data[cat][split] = {} + # BDT inputs + for var in VAR_LST_TMP: bdt_data[cat][split]["base_bdt_" + var] = np.array([]) + # also weight + bdt_data[cat][split]["base_bdt_weight"] = np.array([]) + #bdt_data[cat][split]["sf_bdt_weight"] = np.array([]) + # also event number + bdt_data[cat][split]["base_bdt_event"] = np.array([]) + #bdt_data[cat][split]["sf_bdt_event"] = np.array([]) + + + # Loop over sample groupings + for proc in sample_dict_mc: + + # Grouping handling for background (all backgrounds are grouped as "Bkg") + # proc_cat will be either "WWZ", "ZH", or "Bkg" + #proc_cat = proc + #if proc in ["ZZ", "ttZ", "tWZ", "WZ", "other"]: + if proc in ["ZZ", "ttZ", "VHnobb", "other"]: + proc_cat = "Bkg" + if proc in ["ZZZ", "WZZ"]: + proc_cat = "Sig" + + # Loop over all the groupings and fill the empty numpy array with the actual inputs + for sample_name in sample_dict_mc[proc]: + + # Processing OF BDT inputs + for var in VAR_LST_TMP: + + bdt_vname = "base_bdt_" + var + + # Get the original numpy array + a = bdt_data[proc_cat]["all"][bdt_vname] + + # Get more variables for this sample_name + # masking is done via "bdt_of/sf_proc_list" which holds sample_name (that is how it was saved in wwz4l.py) + b = myd[bdt_vname][myd["list_bdt_base_proc"] == sample_name] + + # Concatenate + c = np.concatenate((a, b)) + + # Set it back to bdt_data where we will store + bdt_data[proc_cat]["all"][bdt_vname] = c + + + # We need to save weights as well + bdt_data[proc_cat]["all"]["base_bdt_weight"] = np.concatenate((bdt_data[proc_cat]["all"]["base_bdt_weight"], myd["list_bdt_base_wgt"][myd["list_bdt_base_proc"] == sample_name])) + + # We need to save event numbers as well + bdt_data[proc_cat]["all"]["base_bdt_event"] = np.concatenate((bdt_data[proc_cat]["all"]["base_bdt_event"], myd["list_bdt_base_evt"][myd["list_bdt_base_proc"] == sample_name])) + + + + # Split the testing and training by event number (even event number vs. odd event number) + for cat in categories: + for key in bdt_data[cat]["all"].keys(): + if "base" in key: + test = bdt_data[cat]["all"][key][bdt_data[cat]["all"]["base_bdt_event"] % 2 == 1] # Getting events with event number that is even + train = bdt_data[cat]["all"][key][bdt_data[cat]["all"]["base_bdt_event"] % 2 == 0] # Getting events with event number that is odd + bdt_data[cat]["test"][key] = test + bdt_data[cat]["train"][key] = train + + + return bdt_data + +def main(): + + # Set up the command line parser + parser = argparse.ArgumentParser() + parser.add_argument("pkl_file_path", help = "The path to the pkl file") + args = parser.parse_args() + + bdt_data = prepare_bdt_training_data(args.pkl_file_path) + + # Format: + #{ + # "Sig" : { + # "all" : { + # "tag_weight" : [], + # "tag_event" : [], + # "tag_var1" : [], + # "tag_var2" : [], + # }, + # "test" : { + # ... + # } + # "train" : { + # ... + # }, + # "Bkg" : { + # ... + # }, + #} + + #print(bdt_data.keys()) + #for k in bdt_data: + # print("\n",k) + # print(k,bdt_data[k].keys()) + # for kk in bdt_data[k]: + # print("\n",kk) + # for kkk in bdt_data[k][kk]: + # print("\nHere:",k,kk,kkk) + # #print(type(bdt_data[k][kk][kkk])) + # #print(len(bdt_data[k][kk][kkk]),bdt_data[k][kk][kkk]) + #exit() + + # Gathering categories and variable names + categories = bdt_data.keys() + base_variables = [] + #sf_variables = [] + for cat in bdt_data.keys(): + for key in bdt_data[cat]["train"].keys(): + if "weight" in key: # skip weight + continue + if "event" in key: # skip weight + continue + if "base_bdt_" in key: + base_variables.append(key) + + # Integer labeling (WWZ = 0, ZH = 1, Bkg = 2) + label_d = {} + #label_d["ZZZ"] = 0 + #label_d["WZZ"] = 1 + #label_d["Bkg"] = 2 + + label_d["Sig"] = 0 + label_d["Bkg"] = 1 + + # Remove duplicates and get only unique variable names + # Since we looped over multiple categories same variable names were put in + base_variables = list(set(base_variables)) + #sf_variables = list(set(sf_variables)) + + # Building XGBoost input (which needs to be in a flat matrix) + X_train_of = [] # Input variables + y_train_of = np.array([]) # Output labels + w_train_of = np.array([]) # Output labels + X_test_of = [] # Input variables + y_test_of = np.array([]) # Output labels + w_test_of = np.array([]) # Output labels + + # Allocate some empty arrays where we will store things + for var in VAR_LST_TMP : X_train_of.append(np.array([])) # Allocate N variable worth of lists + for var in VAR_LST_TMP : X_test_of.append(np.array([])) # Allocate N variable worth of lists + + # Loop over the variables and store + print("bdt_data.keys()",bdt_data.keys()) + for cat in bdt_data.keys(): + # Store the variables + for ivar, var in enumerate(base_variables): + X_train_of[ivar] = np.concatenate((X_train_of[ivar], bdt_data[cat]["train"][var])) + X_test_of[ivar] = np.concatenate((X_test_of[ivar], bdt_data[cat]["test"][var])) + # Take the last variable and get length and create labels + y_train_of = np.concatenate((y_train_of, np.full(len(bdt_data[cat]["train"][var]), label_d[cat]))) + y_test_of = np.concatenate((y_test_of, np.full(len(bdt_data[cat]["test"][var]), label_d[cat]))) + w_train_of = np.concatenate((w_train_of, bdt_data[cat]["train"]["base_bdt_weight"])) + w_test_of = np.concatenate((w_test_of, bdt_data[cat]["test"]["base_bdt_weight"])) + #for ivar, var in enumerate(sf_variables): + #X_train_sf[ivar] = np.concatenate((X_train_sf[ivar], bdt_data[cat]["train"][var])) + #X_test_sf[ivar] = np.concatenate((X_test_sf[ivar], bdt_data[cat]["test"][var])) + + # Turn them into numpy array + X_train_of = np.array(X_train_of) + X_test_of = np.array(X_test_of) + + # Transpose to have the events in rows and variables in columns + X_train_of = X_train_of.T + X_test_of = X_test_of.T + + print("Printing OF training sample numpy array shapes") + print(X_train_of.shape) + print(y_train_of.shape) + print(w_train_of.shape) + print("Printing OF testing sample numpy array shapes") + print(X_test_of.shape) + print(y_test_of.shape) + print(w_test_of.shape) + + ##### TMVA + + ## Building TMVA output numpy array (which needs to be in a flat matrix) + #y_train_tmva_of = [] + #y_test_tmva_of = [] + + ## Allocate some empty arrays where we will store things (there are 3 scores in tmva) + ##tmva_scores_categories = ["wwz", "zh", "bkg"] + #tmva_scores_categories = ["zzz", "wzz", "bkg"] + #for var in range(len(tmva_scores_categories)): y_train_tmva_of.append(np.array([])) # Allocate N variable worth of lists + #for var in range(len(tmva_scores_categories)): y_test_tmva_of.append(np.array([])) # Allocate N variable worth of lists + + ## Loop over the variables and store + #for cat in bdt_data.keys(): + # # Store the variables + # for ivar, var in enumerate(tmva_scores_categories): + # y_train_tmva_of[ivar] = np.concatenate((y_train_tmva_of[ivar], bdt_data[cat]["train"][f"tmva_bdt_base_{var}"])) + # y_test_tmva_of[ivar] = np.concatenate((y_test_tmva_of[ivar], bdt_data[cat]["test"][f"tmva_bdt_base_{var}"])) + # #y_train_tmva_sf[ivar] = np.concatenate((y_train_tmva_sf[ivar], bdt_data[cat]["train"][f"tmva_bdt_sf_{var}"])) + # #y_test_tmva_sf[ivar] = np.concatenate((y_test_tmva_sf[ivar], bdt_data[cat]["test"][f"tmva_bdt_sf_{var}"])) + + #y_train_tmva_of = np.array(y_train_tmva_of) + #y_test_tmva_of = np.array(y_test_tmva_of) + #y_train_tmva_sf = np.array(y_train_tmva_sf) + #y_test_tmva_sf = np.array(y_test_tmva_sf) + + #y_train_tmva_of = y_train_tmva_of.T + #y_test_tmva_of = y_test_tmva_of.T + #y_train_tmva_sf = y_train_tmva_sf.T + #y_test_tmva_sf = y_test_tmva_sf.T + + dd = {} + dd["X_train_of"] = X_train_of + dd["y_train_of"] = y_train_of + dd["w_train_of"] = w_train_of + dd["X_test_of"] = X_test_of + dd["y_test_of"] = y_test_of + dd["w_test_of"] = w_test_of + #dd["y_train_tmva_of"] = y_train_tmva_of + #dd["y_test_tmva_of"] = y_test_tmva_of + + with gzip.open("bdt.pkl.gz", "wb") as fout: + cloudpickle.dump(dd, fout) + +if __name__ == "__main__": + + main() diff --git a/analysis/bdt_training/binary_train.py b/analysis/bdt_training/binary_train.py new file mode 100644 index 00000000..c7ff449d --- /dev/null +++ b/analysis/bdt_training/binary_train.py @@ -0,0 +1,169 @@ +import matplotlib.pyplot as plt +import numpy as np +import mpltern +import pickle +import gzip +from sklearn.metrics import roc_curve, auc +from sklearn.preprocessing import label_binarize +import xgboost as xgb +import os +import shutil + +from topcoffea.scripts.make_html import make_html + +# Avoid import but not used error for mpltern: +# - We need to import it in order to register the "ternary" projection +# - But mpltern is not explicitly used anywhere, making flake8 upset +# - Thus just add an assert line, so that flake8 will not catch it +assert mpltern + +out_dir_name = "outdir_train" +if not os.path.exists(out_dir_name): + os.mkdir(out_dir_name) + +#dd = pickle.load(gzip.open("bdt.pkl.gz")) +#dd = pickle.load(gzip.open("bdt_v0.pkl.gz")) +#dd = pickle.load(gzip.open("bdt_zzz_siphon02.pkl.gz")) +#dd = pickle.load(gzip.open("bdt_zzz_siphon02_18vars.pkl.gz")) +#dd = pickle.load(gzip.open("bdt_siphon04_18vars_j00.pkl.gz")) +dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars.pkl.gz")) + +X_train_of = dd["X_train_of"] +y_train_of = dd["y_train_of"] +w_train_of = dd["w_train_of"] +X_test_of = dd["X_test_of"] +y_test_of = dd["y_test_of"] +w_test_of = dd["w_test_of"] + +print("X_train_of",X_train_of,len(X_train_of)) +print("y_train_of",y_train_of,len(y_train_of)) +print("w_train_of",w_train_of,len(w_train_of)) +print("X_test_of",X_train_of,len(X_test_of)) +print("y_test_of",y_train_of,len(y_test_of)) +print("w_test_of",w_train_of,len(w_test_of)) + + +#xgb_clf = xgb.XGBClassifier(objective='multi:softmax', +xgb_clf = xgb.XGBClassifier( + objective='binary:logistic', + #num_class=3, + #num_class=2, + missing=1, + booster="gbtree", + grow_policy="depthwise", + #learning_rate=2.0, + learning_rate=0.5, + #n_estimators=5, + n_estimators=1000, + #eval_metric=['merror','mlogloss'], + eval_metric=['error','logloss'], + device="cuda", + seed=42, + #n_jobs=32 + n_jobs=64 +) + +params = xgb_clf.get_params() +print("XGBoost Classifier Parameters:\n") +for param, value in params.items(): + print(f"{param}: {value}") + + +xgb_clf.fit(X_train_of, + y_train_of, + sample_weight=np.maximum(w_train_of, 0), + verbose=1, # set to 1 to see xgb training round intermediate results + eval_set=[(X_train_of, y_train_of), (X_test_of, y_test_of)] +) + + +xgb_clf.save_model(f"{out_dir_name}/bdt.json") + +############# +# Hist +p_test_0 = xgb_clf.predict_proba(X_test_of)[:,0] +p_test_1 = xgb_clf.predict_proba(X_test_of)[:,1] +p_train_0 = xgb_clf.predict_proba(X_train_of)[:,0] +p_train_1 = xgb_clf.predict_proba(X_train_of)[:,1] + +is_sig_test = label_binarize(y_test_of, classes=np.unique(y_test_of))[:,0] == 0 +is_bkg_test = label_binarize(y_test_of, classes=np.unique(y_test_of))[:,0] == 1 +is_sig_train = label_binarize(y_train_of, classes=np.unique(y_train_of))[:,0] == 0 +is_bkg_train = label_binarize(y_train_of, classes=np.unique(y_train_of))[:,0] == 1 + +print("p0",p_test_0) +print("p1",p_test_1) +print("is_sig_test",is_sig_test) +print("is_bkg_test",is_bkg_test) +print("s",p_test_0[is_sig_test],max(p_test_0[is_sig_test])) +print("b",p_test_0[is_bkg_test],max(p_test_0[is_bkg_test])) +fig, ax = plt.subplots(figsize=(9,5)) +plt.hist(p_test_1[is_sig_test],bins=100,histtype="step",label="test is_sig",density=True) +plt.hist(p_test_0[is_bkg_test],bins=100,histtype="step",label="test is_bkg",density=True) +plt.hist(p_train_1[is_sig_train],bins=100,histtype="step",label="train is_sig",density=True) +plt.hist(p_train_0[is_bkg_train],bins=100,histtype="step",label="train is_bkg",density=True) +plt.legend() +ax.set_xlim(0,1) +plt.xlabel('prob') +plt.title('Hist of prob') +plt.savefig(f"{out_dir_name}/hist.png") +plt.savefig(f"{out_dir_name}/hist.pdf") +plt.clf() +############# + + +# preparing evaluation metric plots +results = xgb_clf.evals_result() +epochs = len(results['validation_0']['logloss']) +x_axis = range(0, epochs) + +# xgboost 'logloss' plot +fig, ax = plt.subplots(figsize=(9,5)) +ax.plot(x_axis, results['validation_0']['logloss'], label='Train') +ax.plot(x_axis, results['validation_1']['logloss'], label='Test') +ax.legend() +plt.ylabel('logloss') +plt.title('GridSearchCV XGBoost logloss') +plt.savefig(f"{out_dir_name}/logloss.png") +plt.savefig(f"{out_dir_name}/logloss.pdf") +plt.clf() + +# xgboost 'error' plot +fig, ax = plt.subplots(figsize=(9,5)) +ax.plot(x_axis, results['validation_0']['error'], label='Train') +ax.plot(x_axis, results['validation_1']['error'], label='Test') +ax.legend() +plt.ylabel('error') +plt.title('GridSearchCV XGBoost error') +plt.savefig(f"{out_dir_name}/error.png") +plt.savefig(f"{out_dir_name}/error.pdf") +plt.clf() + +### ROC +y_prob_of = xgb_clf.predict_proba(X_test_of) +n_classes = len(np.unique(y_test_of)) +fpr = {} +tpr = {} +roc_auc = {} +y_test_of_bin = label_binarize(y_test_of, classes=np.unique(y_test_of)) + +fpr, tpr, _ = roc_curve(label_binarize(y_test_of,classes=np.unique(y_test_of)), p_test_1, sample_weight=np.abs(w_test_of)) +roc_auc = auc(fpr, tpr) +plt.plot(fpr, tpr, lw=2, label=f'ROC curve test (AUC = {roc_auc})') +fpr, tpr, _ = roc_curve(label_binarize(y_train_of,classes=np.unique(y_train_of)), p_train_1, sample_weight=np.abs(w_train_of)) +roc_auc = auc(fpr, tpr) +plt.plot(fpr, tpr, lw=2, label=f'ROC curve train (AUC = {roc_auc})') + +plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--') +plt.xlim([0.0, 1.0]) +plt.ylim([0.0, 1.05]) +plt.xlabel('False Positive Rate') +plt.ylabel('True Positive Rate') +plt.title('Receiver Operating Characteristic (ROC) Curve') +plt.legend(loc="lower right") +plt.savefig(f"{out_dir_name}/roc.png") +plt.savefig(f"{out_dir_name}/roc.pdf") +### + +shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) + From cfa937e7eeb12a245359041bbd7c02ca365dad6a Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Thu, 10 Apr 2025 13:04:14 -0700 Subject: [PATCH 3/7] Clean up --- analysis/bdt_training/binary_train.py | 316 +++++++++++++------------- 1 file changed, 159 insertions(+), 157 deletions(-) diff --git a/analysis/bdt_training/binary_train.py b/analysis/bdt_training/binary_train.py index c7ff449d..f44ef919 100644 --- a/analysis/bdt_training/binary_train.py +++ b/analysis/bdt_training/binary_train.py @@ -1,6 +1,5 @@ import matplotlib.pyplot as plt import numpy as np -import mpltern import pickle import gzip from sklearn.metrics import roc_curve, auc @@ -11,159 +10,162 @@ from topcoffea.scripts.make_html import make_html -# Avoid import but not used error for mpltern: -# - We need to import it in order to register the "ternary" projection -# - But mpltern is not explicitly used anywhere, making flake8 upset -# - Thus just add an assert line, so that flake8 will not catch it -assert mpltern - -out_dir_name = "outdir_train" -if not os.path.exists(out_dir_name): - os.mkdir(out_dir_name) - -#dd = pickle.load(gzip.open("bdt.pkl.gz")) -#dd = pickle.load(gzip.open("bdt_v0.pkl.gz")) -#dd = pickle.load(gzip.open("bdt_zzz_siphon02.pkl.gz")) -#dd = pickle.load(gzip.open("bdt_zzz_siphon02_18vars.pkl.gz")) -#dd = pickle.load(gzip.open("bdt_siphon04_18vars_j00.pkl.gz")) -dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars.pkl.gz")) - -X_train_of = dd["X_train_of"] -y_train_of = dd["y_train_of"] -w_train_of = dd["w_train_of"] -X_test_of = dd["X_test_of"] -y_test_of = dd["y_test_of"] -w_test_of = dd["w_test_of"] - -print("X_train_of",X_train_of,len(X_train_of)) -print("y_train_of",y_train_of,len(y_train_of)) -print("w_train_of",w_train_of,len(w_train_of)) -print("X_test_of",X_train_of,len(X_test_of)) -print("y_test_of",y_train_of,len(y_test_of)) -print("w_test_of",w_train_of,len(w_test_of)) - - -#xgb_clf = xgb.XGBClassifier(objective='multi:softmax', -xgb_clf = xgb.XGBClassifier( - objective='binary:logistic', - #num_class=3, - #num_class=2, - missing=1, - booster="gbtree", - grow_policy="depthwise", - #learning_rate=2.0, - learning_rate=0.5, - #n_estimators=5, - n_estimators=1000, - #eval_metric=['merror','mlogloss'], - eval_metric=['error','logloss'], - device="cuda", - seed=42, - #n_jobs=32 - n_jobs=64 -) - -params = xgb_clf.get_params() -print("XGBoost Classifier Parameters:\n") -for param, value in params.items(): - print(f"{param}: {value}") - - -xgb_clf.fit(X_train_of, - y_train_of, - sample_weight=np.maximum(w_train_of, 0), - verbose=1, # set to 1 to see xgb training round intermediate results - eval_set=[(X_train_of, y_train_of), (X_test_of, y_test_of)] -) - - -xgb_clf.save_model(f"{out_dir_name}/bdt.json") - -############# -# Hist -p_test_0 = xgb_clf.predict_proba(X_test_of)[:,0] -p_test_1 = xgb_clf.predict_proba(X_test_of)[:,1] -p_train_0 = xgb_clf.predict_proba(X_train_of)[:,0] -p_train_1 = xgb_clf.predict_proba(X_train_of)[:,1] - -is_sig_test = label_binarize(y_test_of, classes=np.unique(y_test_of))[:,0] == 0 -is_bkg_test = label_binarize(y_test_of, classes=np.unique(y_test_of))[:,0] == 1 -is_sig_train = label_binarize(y_train_of, classes=np.unique(y_train_of))[:,0] == 0 -is_bkg_train = label_binarize(y_train_of, classes=np.unique(y_train_of))[:,0] == 1 - -print("p0",p_test_0) -print("p1",p_test_1) -print("is_sig_test",is_sig_test) -print("is_bkg_test",is_bkg_test) -print("s",p_test_0[is_sig_test],max(p_test_0[is_sig_test])) -print("b",p_test_0[is_bkg_test],max(p_test_0[is_bkg_test])) -fig, ax = plt.subplots(figsize=(9,5)) -plt.hist(p_test_1[is_sig_test],bins=100,histtype="step",label="test is_sig",density=True) -plt.hist(p_test_0[is_bkg_test],bins=100,histtype="step",label="test is_bkg",density=True) -plt.hist(p_train_1[is_sig_train],bins=100,histtype="step",label="train is_sig",density=True) -plt.hist(p_train_0[is_bkg_train],bins=100,histtype="step",label="train is_bkg",density=True) -plt.legend() -ax.set_xlim(0,1) -plt.xlabel('prob') -plt.title('Hist of prob') -plt.savefig(f"{out_dir_name}/hist.png") -plt.savefig(f"{out_dir_name}/hist.pdf") -plt.clf() -############# - - -# preparing evaluation metric plots -results = xgb_clf.evals_result() -epochs = len(results['validation_0']['logloss']) -x_axis = range(0, epochs) - -# xgboost 'logloss' plot -fig, ax = plt.subplots(figsize=(9,5)) -ax.plot(x_axis, results['validation_0']['logloss'], label='Train') -ax.plot(x_axis, results['validation_1']['logloss'], label='Test') -ax.legend() -plt.ylabel('logloss') -plt.title('GridSearchCV XGBoost logloss') -plt.savefig(f"{out_dir_name}/logloss.png") -plt.savefig(f"{out_dir_name}/logloss.pdf") -plt.clf() - -# xgboost 'error' plot -fig, ax = plt.subplots(figsize=(9,5)) -ax.plot(x_axis, results['validation_0']['error'], label='Train') -ax.plot(x_axis, results['validation_1']['error'], label='Test') -ax.legend() -plt.ylabel('error') -plt.title('GridSearchCV XGBoost error') -plt.savefig(f"{out_dir_name}/error.png") -plt.savefig(f"{out_dir_name}/error.pdf") -plt.clf() - -### ROC -y_prob_of = xgb_clf.predict_proba(X_test_of) -n_classes = len(np.unique(y_test_of)) -fpr = {} -tpr = {} -roc_auc = {} -y_test_of_bin = label_binarize(y_test_of, classes=np.unique(y_test_of)) - -fpr, tpr, _ = roc_curve(label_binarize(y_test_of,classes=np.unique(y_test_of)), p_test_1, sample_weight=np.abs(w_test_of)) -roc_auc = auc(fpr, tpr) -plt.plot(fpr, tpr, lw=2, label=f'ROC curve test (AUC = {roc_auc})') -fpr, tpr, _ = roc_curve(label_binarize(y_train_of,classes=np.unique(y_train_of)), p_train_1, sample_weight=np.abs(w_train_of)) -roc_auc = auc(fpr, tpr) -plt.plot(fpr, tpr, lw=2, label=f'ROC curve train (AUC = {roc_auc})') - -plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--') -plt.xlim([0.0, 1.0]) -plt.ylim([0.0, 1.05]) -plt.xlabel('False Positive Rate') -plt.ylabel('True Positive Rate') -plt.title('Receiver Operating Characteristic (ROC) Curve') -plt.legend(loc="lower right") -plt.savefig(f"{out_dir_name}/roc.png") -plt.savefig(f"{out_dir_name}/roc.pdf") -### - -shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) - +def main(): + + ### Get input data ### + + out_dir_name = "outdir_train" + if not os.path.exists(out_dir_name): + os.mkdir(out_dir_name) + + #dd = pickle.load(gzip.open("bdt.pkl.gz")) + #dd = pickle.load(gzip.open("bdt_v0.pkl.gz")) + #dd = pickle.load(gzip.open("bdt_zzz_siphon02.pkl.gz")) + #dd = pickle.load(gzip.open("bdt_zzz_siphon02_18vars.pkl.gz")) + #dd = pickle.load(gzip.open("bdt_siphon04_18vars_j00.pkl.gz")) + dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars.pkl.gz")) + + X_train_of = dd["X_train_of"] + y_train_of = dd["y_train_of"] + w_train_of = dd["w_train_of"] + X_test_of = dd["X_test_of"] + y_test_of = dd["y_test_of"] + w_test_of = dd["w_test_of"] + + print("X_train_of",X_train_of,len(X_train_of)) + print("y_train_of",y_train_of,len(y_train_of)) + print("w_train_of",w_train_of,len(w_train_of)) + print("X_test_of",X_train_of,len(X_test_of)) + print("y_test_of",y_train_of,len(y_test_of)) + print("w_test_of",w_train_of,len(w_test_of)) + + + ### Define model and fit ### + + xgb_clf = xgb.XGBClassifier( + objective='binary:logistic', + missing=1, + booster="gbtree", + grow_policy="depthwise", + #learning_rate=2.0, + learning_rate=0.5, + n_estimators=5, + #n_estimators=1000, + eval_metric=['error','logloss'], + device="cuda", + seed=42, + #n_jobs=32 + n_jobs=64 + ) + params = xgb_clf.get_params() + print("XGBoost Classifier Parameters:\n") + for param, value in params.items(): + print(f" {param}: {value}") + + # Train + xgb_clf.fit( + X_train_of, + y_train_of, + sample_weight=np.maximum(w_train_of, 0), + verbose=1, # set to 1 to see xgb training round intermediate results + eval_set=[(X_train_of, y_train_of), (X_test_of, y_test_of)] + ) + + xgb_clf.save_model(f"{out_dir_name}/bdt.json") + + + ############# Make resutls plots ############# + + shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) + + ### Make hists ### + + p_test_0 = xgb_clf.predict_proba(X_test_of)[:,0] + p_test_1 = xgb_clf.predict_proba(X_test_of)[:,1] + p_train_0 = xgb_clf.predict_proba(X_train_of)[:,0] + p_train_1 = xgb_clf.predict_proba(X_train_of)[:,1] + + is_sig_test = label_binarize(y_test_of, classes=np.unique(y_test_of))[:,0] == 0 + is_bkg_test = label_binarize(y_test_of, classes=np.unique(y_test_of))[:,0] == 1 + is_sig_train = label_binarize(y_train_of, classes=np.unique(y_train_of))[:,0] == 0 + is_bkg_train = label_binarize(y_train_of, classes=np.unique(y_train_of))[:,0] == 1 + + print("p0",p_test_0) + print("p1",p_test_1) + print("is_sig_test",is_sig_test) + print("is_bkg_test",is_bkg_test) + print("s",p_test_0[is_sig_test],max(p_test_0[is_sig_test])) + print("b",p_test_0[is_bkg_test],max(p_test_0[is_bkg_test])) + + fig, ax = plt.subplots(figsize=(9,5)) + plt.hist(p_train_1[is_sig_train],bins=100,histtype="step",label="train is_sig",density=True) + plt.hist(p_test_1[is_sig_test],bins=100,histtype="step",label="test is_sig",density=True) + plt.hist(p_train_0[is_bkg_train],bins=100,histtype="step",label="train is_bkg",density=True) + plt.hist(p_test_0[is_bkg_test],bins=100,histtype="step",label="test is_bkg",density=True) + + plt.legend() + ax.set_xlim(0,1) + plt.xlabel('prob') + plt.title('Hist of prob') + plt.savefig(f"{out_dir_name}/hist.png") + plt.savefig(f"{out_dir_name}/hist.pdf") + plt.clf() + + + ### Make metric plots ### + + # preparing evaluation metric plots + results = xgb_clf.evals_result() + epochs = len(results['validation_0']['logloss']) + x_axis = range(0, epochs) + + # xgboost 'logloss' plot + fig, ax = plt.subplots(figsize=(9,5)) + ax.plot(x_axis, results['validation_0']['logloss'], label='Train') + ax.plot(x_axis, results['validation_1']['logloss'], label='Test') + ax.legend() + plt.ylabel('logloss') + plt.title('GridSearchCV XGBoost logloss') + plt.savefig(f"{out_dir_name}/logloss.png") + plt.savefig(f"{out_dir_name}/logloss.pdf") + plt.clf() + + # xgboost 'error' plot + fig, ax = plt.subplots(figsize=(9,5)) + ax.plot(x_axis, results['validation_0']['error'], label='Train') + ax.plot(x_axis, results['validation_1']['error'], label='Test') + ax.legend() + plt.ylabel('error') + plt.title('GridSearchCV XGBoost error') + plt.savefig(f"{out_dir_name}/error.png") + plt.savefig(f"{out_dir_name}/error.pdf") + plt.clf() + + ### Make ROC plots ### + + y_prob_of = xgb_clf.predict_proba(X_test_of) + n_classes = len(np.unique(y_test_of)) + fpr = {} + tpr = {} + roc_auc = {} + y_test_of_bin = label_binarize(y_test_of, classes=np.unique(y_test_of)) + + fpr, tpr, _ = roc_curve(label_binarize(y_train_of,classes=np.unique(y_train_of)), p_train_1, sample_weight=np.abs(w_train_of)) + roc_auc = auc(fpr, tpr) + plt.plot(fpr, tpr, lw=2, label=f'ROC curve train (AUC = {roc_auc})') + fpr, tpr, _ = roc_curve(label_binarize(y_test_of,classes=np.unique(y_test_of)), p_test_1, sample_weight=np.abs(w_test_of)) + roc_auc = auc(fpr, tpr) + plt.plot(fpr, tpr, lw=2, label=f'ROC curve test (AUC = {roc_auc})') + + plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--') + plt.xlim([0.0, 1.0]) + plt.ylim([0.0, 1.05]) + plt.xlabel('False Positive Rate') + plt.ylabel('True Positive Rate') + plt.title('Receiver Operating Characteristic (ROC) Curve') + plt.legend(loc="lower right") + plt.savefig(f"{out_dir_name}/roc.png") + plt.savefig(f"{out_dir_name}/roc.pdf") + + +main() From d175766c1bd97651159fdbd2427b4335fc881ab6 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Thu, 10 Apr 2025 18:06:18 -0700 Subject: [PATCH 4/7] Adding code for variable plotting --- analysis/bdt_training/binary_train.py | 84 ++++++++++++++++++++++++++- 1 file changed, 81 insertions(+), 3 deletions(-) diff --git a/analysis/bdt_training/binary_train.py b/analysis/bdt_training/binary_train.py index f44ef919..c8a345a7 100644 --- a/analysis/bdt_training/binary_train.py +++ b/analysis/bdt_training/binary_train.py @@ -10,9 +10,30 @@ from topcoffea.scripts.make_html import make_html +# Take an array find a cap where 99.5% of vals are below +def get_cap_val(in_arr): + + cap_val = None + min_val = min(in_arr) + max_val = max(in_arr) + nsteps = 10 + step_size = abs(max_val-min_val)/nsteps + for iterator in range(nsteps): + val = min_val + iterator*step_size + tot = len(in_arr) + sub_tot = len(in_arr[in_arr 0.999: + cap_val = val + break + + if cap_val is None: cap_val = max(in_arr) + + return cap_val + + def main(): - ### Get input data ### + ############# Get input data ############# out_dir_name = "outdir_train" if not os.path.exists(out_dir_name): @@ -40,7 +61,62 @@ def main(): print("w_test_of",w_train_of,len(w_test_of)) - ### Define model and fit ### + ############# Plot input vars ############# + + shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) + + nvars = 16 + sig_key = 0 + bkg_key = 1 + + for i in range(nvars): + print(i) + + ### Get the capped arrays ### + # Train + var_train_i = X_train_of[:,i] + cap_val_train = get_cap_val(var_train_i) + # Test + var_test_i = X_test_of[:,i] + cap_val_test = get_cap_val(var_test_i) + # Cap + cap_val = max(cap_val_test,cap_val_train) + var_train_i = np.where(var_train_i Date: Thu, 10 Apr 2025 20:22:42 -0700 Subject: [PATCH 5/7] var lst --- analysis/bdt_training/binary_ingress.py | 81 +++++++++++-------------- analysis/bdt_training/binary_train.py | 59 +++++++++--------- 2 files changed, 64 insertions(+), 76 deletions(-) diff --git a/analysis/bdt_training/binary_ingress.py b/analysis/bdt_training/binary_ingress.py index b5ec8239..7284aadf 100644 --- a/analysis/bdt_training/binary_ingress.py +++ b/analysis/bdt_training/binary_ingress.py @@ -12,8 +12,28 @@ "j0pt", "nleps", ] - VAR_LST_TMP = [ + "ptl4", + "scalarptsum_lepmetjet", + "l0pt", + "l1pt", + "l2pt", + "l3pt", + "j0pt", + "j1pt", + "absdphi_zl0_zl1", + "absdphi_wl0_wl1", + "mll_wl0_wl1", + "mll_zl0_zl1", + "mll_min_sfos", + "met", + "njets", + "nbtagsl", + "mjj_j0j1", + "mjj_nearz", +] + +_VAR_LST_TMP = [ #"absdphi_4l_met", #"absdphi_min_afas", #"absdphi_min_afos", @@ -86,6 +106,11 @@ #"z_lep1_pt", ] +VAR_LST_TMP_FULL = [] +for var_base in VAR_LST_TMP: + VAR_LST_TMP_FULL.append("base_bdt_"+var_base) + + def prepare_bdt_training_data(pklfilepath): # Open up the pickle file @@ -104,15 +129,14 @@ def prepare_bdt_training_data(pklfilepath): # Parse the pickle file and store to myd as numpy arrays for list_output_name in list_output_names: - myd[list_output_name] = np.array(d[list_output_name]) + myd[list_output_name] = np.array(d[list_output_name]) for var in VAR_LST_TMP: - myd["base_bdt_" + var] = np.array(d["list_base_bdt_" + var]) + myd["base_bdt_" + var] = np.array(d["list_base_bdt_" + var]) # Sample groupings to split the data into signals, and backgrounds sample_dict_mc = sg.create_mc_sample_dict("run2") # List of categories that will use BDT to label - #categories = ["WWZ", "ZH", "Bkg"] #categories = ["ZZZ", "WZZ", "Bkg"] categories = ["Sig", "Bkg"] @@ -144,9 +168,6 @@ def prepare_bdt_training_data(pklfilepath): for proc in sample_dict_mc: # Grouping handling for background (all backgrounds are grouped as "Bkg") - # proc_cat will be either "WWZ", "ZH", or "Bkg" - #proc_cat = proc - #if proc in ["ZZ", "ttZ", "tWZ", "WZ", "other"]: if proc in ["ZZ", "ttZ", "VHnobb", "other"]: proc_cat = "Bkg" if proc in ["ZZZ", "WZZ"]: @@ -185,6 +206,7 @@ def prepare_bdt_training_data(pklfilepath): # Split the testing and training by event number (even event number vs. odd event number) for cat in categories: for key in bdt_data[cat]["all"].keys(): + print("key",key) if "base" in key: test = bdt_data[cat]["all"][key][bdt_data[cat]["all"]["base_bdt_event"] % 2 == 1] # Getting events with event number that is even train = bdt_data[cat]["all"][key][bdt_data[cat]["all"]["base_bdt_event"] % 2 == 0] # Getting events with event number that is odd @@ -275,10 +297,13 @@ def main(): for var in VAR_LST_TMP : X_test_of.append(np.array([])) # Allocate N variable worth of lists # Loop over the variables and store - print("bdt_data.keys()",bdt_data.keys()) + final_var_lst = [] + print(base_variables) + print(len(base_variables)) for cat in bdt_data.keys(): # Store the variables - for ivar, var in enumerate(base_variables): + for ivar, var in enumerate(VAR_LST_TMP_FULL): + final_var_lst.append(var) X_train_of[ivar] = np.concatenate((X_train_of[ivar], bdt_data[cat]["train"][var])) X_test_of[ivar] = np.concatenate((X_test_of[ivar], bdt_data[cat]["test"][var])) # Take the last variable and get length and create labels @@ -286,9 +311,6 @@ def main(): y_test_of = np.concatenate((y_test_of, np.full(len(bdt_data[cat]["test"][var]), label_d[cat]))) w_train_of = np.concatenate((w_train_of, bdt_data[cat]["train"]["base_bdt_weight"])) w_test_of = np.concatenate((w_test_of, bdt_data[cat]["test"]["base_bdt_weight"])) - #for ivar, var in enumerate(sf_variables): - #X_train_sf[ivar] = np.concatenate((X_train_sf[ivar], bdt_data[cat]["train"][var])) - #X_test_sf[ivar] = np.concatenate((X_test_sf[ivar], bdt_data[cat]["test"][var])) # Turn them into numpy array X_train_of = np.array(X_train_of) @@ -307,37 +329,6 @@ def main(): print(y_test_of.shape) print(w_test_of.shape) - ##### TMVA - - ## Building TMVA output numpy array (which needs to be in a flat matrix) - #y_train_tmva_of = [] - #y_test_tmva_of = [] - - ## Allocate some empty arrays where we will store things (there are 3 scores in tmva) - ##tmva_scores_categories = ["wwz", "zh", "bkg"] - #tmva_scores_categories = ["zzz", "wzz", "bkg"] - #for var in range(len(tmva_scores_categories)): y_train_tmva_of.append(np.array([])) # Allocate N variable worth of lists - #for var in range(len(tmva_scores_categories)): y_test_tmva_of.append(np.array([])) # Allocate N variable worth of lists - - ## Loop over the variables and store - #for cat in bdt_data.keys(): - # # Store the variables - # for ivar, var in enumerate(tmva_scores_categories): - # y_train_tmva_of[ivar] = np.concatenate((y_train_tmva_of[ivar], bdt_data[cat]["train"][f"tmva_bdt_base_{var}"])) - # y_test_tmva_of[ivar] = np.concatenate((y_test_tmva_of[ivar], bdt_data[cat]["test"][f"tmva_bdt_base_{var}"])) - # #y_train_tmva_sf[ivar] = np.concatenate((y_train_tmva_sf[ivar], bdt_data[cat]["train"][f"tmva_bdt_sf_{var}"])) - # #y_test_tmva_sf[ivar] = np.concatenate((y_test_tmva_sf[ivar], bdt_data[cat]["test"][f"tmva_bdt_sf_{var}"])) - - #y_train_tmva_of = np.array(y_train_tmva_of) - #y_test_tmva_of = np.array(y_test_tmva_of) - #y_train_tmva_sf = np.array(y_train_tmva_sf) - #y_test_tmva_sf = np.array(y_test_tmva_sf) - - #y_train_tmva_of = y_train_tmva_of.T - #y_test_tmva_of = y_test_tmva_of.T - #y_train_tmva_sf = y_train_tmva_sf.T - #y_test_tmva_sf = y_test_tmva_sf.T - dd = {} dd["X_train_of"] = X_train_of dd["y_train_of"] = y_train_of @@ -345,8 +336,8 @@ def main(): dd["X_test_of"] = X_test_of dd["y_test_of"] = y_test_of dd["w_test_of"] = w_test_of - #dd["y_train_tmva_of"] = y_train_tmva_of - #dd["y_test_tmva_of"] = y_test_tmva_of + + dd["var_name_lst"] = VAR_LST_TMP_FULL with gzip.open("bdt.pkl.gz", "wb") as fout: cloudpickle.dump(dd, fout) diff --git a/analysis/bdt_training/binary_train.py b/analysis/bdt_training/binary_train.py index c8a345a7..0ff9d38d 100644 --- a/analysis/bdt_training/binary_train.py +++ b/analysis/bdt_training/binary_train.py @@ -39,12 +39,8 @@ def main(): if not os.path.exists(out_dir_name): os.mkdir(out_dir_name) - #dd = pickle.load(gzip.open("bdt.pkl.gz")) - #dd = pickle.load(gzip.open("bdt_v0.pkl.gz")) - #dd = pickle.load(gzip.open("bdt_zzz_siphon02.pkl.gz")) - #dd = pickle.load(gzip.open("bdt_zzz_siphon02_18vars.pkl.gz")) - #dd = pickle.load(gzip.open("bdt_siphon04_18vars_j00.pkl.gz")) - dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars.pkl.gz")) + #dd = pickle.load(gzip.open("bdt_zzz_siphon05_16vars_k_ordered.pkl.gz")) + dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars_p_ordered.pkl.gz")) X_train_of = dd["X_train_of"] y_train_of = dd["y_train_of"] @@ -53,24 +49,27 @@ def main(): y_test_of = dd["y_test_of"] w_test_of = dd["w_test_of"] - print("X_train_of",X_train_of,len(X_train_of)) - print("y_train_of",y_train_of,len(y_train_of)) - print("w_train_of",w_train_of,len(w_train_of)) - print("X_test_of",X_train_of,len(X_test_of)) - print("y_test_of",y_train_of,len(y_test_of)) - print("w_test_of",w_train_of,len(w_test_of)) + var_name_lst = dd["var_name_lst"] + + print("X_train_of",X_train_of.shape) + print("y_train_of",y_train_of.shape) + print("w_train_of",w_train_of.shape) + print("X_test_of" ,X_test_of.shape) + print("y_test_of" ,y_test_of.shape) + print("w_test_of" ,w_test_of.shape) ############# Plot input vars ############# shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) - nvars = 16 + nvars = len(var_name_lst) sig_key = 0 bkg_key = 1 - for i in range(nvars): - print(i) + print("var_name_lst",var_name_lst) + for i,var_name in enumerate(var_name_lst): + print(i,var_name) ### Get the capped arrays ### # Train @@ -106,16 +105,14 @@ def main(): plt.hist(var_train_i_bkg,weights=wgt_train_bkg,bins=60,range=hrange,histtype="step",label="train bkg",density=True) plt.hist(var_test_i_bkg, weights=wgt_test_bkg, bins=60,range=hrange,histtype="step",label="test bkg",density=True) plt.legend() - plt.xlabel('mystery variable') - plt.title(f'Variable {i}, density=T') - plt.savefig(f"{out_dir_name}/var_{i}.png") - plt.savefig(f"{out_dir_name}/var_{i}.pdf") + plt.xlabel(f'{var_name}') + plt.title(f'Variable {i} {var_name}, density=T') + plt.savefig(f"{out_dir_name}/var_{i}_{var_name}.png") + plt.savefig(f"{out_dir_name}/var_{i}_{var_name}.pdf") plt.clf() - #break - + #exit() - exit() ############# Define model and fit ############# xgb_clf = xgb.XGBClassifier( @@ -125,13 +122,13 @@ def main(): grow_policy="depthwise", #learning_rate=2.0, learning_rate=0.5, - n_estimators=5, - #n_estimators=1000, + #n_estimators=100, + n_estimators=800, eval_metric=['error','logloss'], device="cuda", seed=42, - #n_jobs=32 - n_jobs=64 + n_jobs=32 + #n_jobs=64 ) params = xgb_clf.get_params() print("XGBoost Classifier Parameters:\n") @@ -153,7 +150,7 @@ def main(): ############# Make resutls plots ############# - ### Make hists ### + ### Make prob hist ### p_test_0 = xgb_clf.predict_proba(X_test_of)[:,0] p_test_1 = xgb_clf.predict_proba(X_test_of)[:,1] @@ -172,7 +169,7 @@ def main(): print("s",p_test_0[is_sig_test],max(p_test_0[is_sig_test])) print("b",p_test_0[is_bkg_test],max(p_test_0[is_bkg_test])) - fig, ax = plt.subplots(figsize=(9,5)) + fig, ax = plt.subplots(figsize=(5,5)) plt.hist(p_train_1[is_sig_train],bins=100,histtype="step",label="train is_sig",density=True) plt.hist(p_test_1[is_sig_test],bins=100,histtype="step",label="test is_sig",density=True) plt.hist(p_train_0[is_bkg_train],bins=100,histtype="step",label="train is_bkg",density=True) @@ -195,7 +192,7 @@ def main(): x_axis = range(0, epochs) # xgboost 'logloss' plot - fig, ax = plt.subplots(figsize=(9,5)) + fig, ax = plt.subplots(figsize=(5,5)) ax.plot(x_axis, results['validation_0']['logloss'], label='Train') ax.plot(x_axis, results['validation_1']['logloss'], label='Test') ax.legend() @@ -209,7 +206,7 @@ def main(): plt.clf() # xgboost 'error' plot - fig, ax = plt.subplots(figsize=(9,5)) + fig, ax = plt.subplots(figsize=(5,5)) ax.plot(x_axis, results['validation_0']['error'], label='Train') ax.plot(x_axis, results['validation_1']['error'], label='Test') ax.legend() @@ -219,7 +216,7 @@ def main(): plt.savefig(f"{out_dir_name}/error.pdf") plt.clf() - ### Make ROC plots ### + ### Make ROC plot ### y_prob_of = xgb_clf.predict_proba(X_test_of) n_classes = len(np.unique(y_test_of)) From 8e0495640696afe1c821242cb44e11ad84bc92c3 Mon Sep 17 00:00:00 2001 From: Kelci Mohrman Date: Fri, 11 Apr 2025 12:05:54 -0700 Subject: [PATCH 6/7] Bkg 0 sig 1 --- analysis/bdt_training/binary_ingress.py | 4 +- analysis/bdt_training/binary_train.py | 160 ++++++++++++------------ 2 files changed, 80 insertions(+), 84 deletions(-) diff --git a/analysis/bdt_training/binary_ingress.py b/analysis/bdt_training/binary_ingress.py index 7284aadf..399d86b9 100644 --- a/analysis/bdt_training/binary_ingress.py +++ b/analysis/bdt_training/binary_ingress.py @@ -276,8 +276,8 @@ def main(): #label_d["WZZ"] = 1 #label_d["Bkg"] = 2 - label_d["Sig"] = 0 - label_d["Bkg"] = 1 + label_d["Sig"] = 1 + label_d["Bkg"] = 0 # Remove duplicates and get only unique variable names # Since we looped over multiple categories same variable names were put in diff --git a/analysis/bdt_training/binary_train.py b/analysis/bdt_training/binary_train.py index 0ff9d38d..fa6a8154 100644 --- a/analysis/bdt_training/binary_train.py +++ b/analysis/bdt_training/binary_train.py @@ -39,8 +39,8 @@ def main(): if not os.path.exists(out_dir_name): os.mkdir(out_dir_name) - #dd = pickle.load(gzip.open("bdt_zzz_siphon05_16vars_k_ordered.pkl.gz")) - dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars_p_ordered.pkl.gz")) + #dd = pickle.load(gzip.open("bdt.pkl.gz")) + dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars_p_ordered_b0s1_2.pkl.gz")) X_train_of = dd["X_train_of"] y_train_of = dd["y_train_of"] @@ -62,73 +62,76 @@ def main(): ############# Plot input vars ############# shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) - - nvars = len(var_name_lst) - sig_key = 0 - bkg_key = 1 - print("var_name_lst",var_name_lst) - for i,var_name in enumerate(var_name_lst): - print(i,var_name) - - ### Get the capped arrays ### - # Train - var_train_i = X_train_of[:,i] - cap_val_train = get_cap_val(var_train_i) - # Test - var_test_i = X_test_of[:,i] - cap_val_test = get_cap_val(var_test_i) - # Cap - cap_val = max(cap_val_test,cap_val_train) - var_train_i = np.where(var_train_i Date: Sat, 12 Apr 2025 08:10:10 -0700 Subject: [PATCH 7/7] Better prints --- analysis/bdt_training/binary_train.py | 124 +++++++++++++++++++------- 1 file changed, 93 insertions(+), 31 deletions(-) diff --git a/analysis/bdt_training/binary_train.py b/analysis/bdt_training/binary_train.py index fa6a8154..f2e5f3a6 100644 --- a/analysis/bdt_training/binary_train.py +++ b/analysis/bdt_training/binary_train.py @@ -3,7 +3,6 @@ import pickle import gzip from sklearn.metrics import roc_curve, auc -from sklearn.preprocessing import label_binarize import xgboost as xgb import os import shutil @@ -40,7 +39,8 @@ def main(): os.mkdir(out_dir_name) #dd = pickle.load(gzip.open("bdt.pkl.gz")) - dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars_p_ordered_b0s1_2.pkl.gz")) + #dd = pickle.load(gzip.open("bdt_zzz_siphon05_18vars_p_ordered_b0s1_2.pkl.gz")) + dd = pickle.load(gzip.open("bdt_zzz_siphon06_18vars_p_ordered.pkl.gz")) X_train_of = dd["X_train_of"] y_train_of = dd["y_train_of"] @@ -51,20 +51,31 @@ def main(): var_name_lst = dd["var_name_lst"] - print("X_train_of",X_train_of.shape) - print("y_train_of",y_train_of.shape) - print("w_train_of",w_train_of.shape) - print("X_test_of" ,X_test_of.shape) - print("y_test_of" ,y_test_of.shape) - print("w_test_of" ,w_test_of.shape) + # Print out shape info and write to file + print("Shapes:") + print(f"\tX_train_of: {X_train_of.shape}") + print(f"\ty_train_of: {y_train_of.shape}") + print(f"\tw_train_of: {w_train_of.shape}") + print(f"\tX_test_of: {X_test_of.shape}") + print(f"\ty_test_of: {y_test_of.shape}") + print(f"\tw_test_of: {w_test_of.shape}") + print(f"\tvar_name_lst: {var_name_lst}") + fout = open(f"{out_dir_name}/info.txt", "w") + fout.write(f"\nShapes:\n") + fout.write(f"\tX_train_of: {X_train_of.shape}\n") + fout.write(f"\ty_train_of: {y_train_of.shape}\n") + fout.write(f"\tw_train_of: {w_train_of.shape}\n") + fout.write(f"\tX_test_of: {X_test_of.shape}\n") + fout.write(f"\ty_test_of: {y_test_of.shape}\n") + fout.write(f"\tw_test_of: {w_test_of.shape}\n") + fout.write(f"var_name_lst: {var_name_lst}\n") + fout.close() + shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) ############# Plot input vars ############# - shutil.copyfile("/home/users/phchang/public_html/dump/forKelci/index.php.txt", os.path.join(out_dir_name,"index.php")) - print("var_name_lst",var_name_lst) - - make_input_var_plots = 1 + make_input_var_plots = 0 if make_input_var_plots: nvars = len(var_name_lst) sig_key = 1 @@ -119,24 +130,54 @@ def main(): ############# Define model and fit ############# xgb_clf = xgb.XGBClassifier( + #objective='multi:softprob', + #objective='multi:softmax', + #num_class=2, + #objective='reg:squarederror', + #objective='reg:logistic', + objective='binary:logistic', - #missing=1, - booster="gbtree", - grow_policy="depthwise", - learning_rate=0.05, - n_estimators=800, - #n_estimators=1500, + n_estimators=500, eval_metric=['error','logloss'], - device="cuda", + #eval_metric=['merror','mlogloss'], seed=42, - #n_jobs=32, n_jobs=64, - max_depth=6, + scale_pos_weight=400, + + #device="cuda", + #grow_policy="depthwise", + #booster="gbtree", + + learning_rate=0.2, + reg_lambda=1.5, + reg_alpha=1.5, + max_depth=4, + min_child_weight=2, + #min_split_loss=100, # Does not seem to help + #max_delta_step=10, # Seems to not do anything + + subsample=0.5, ) + + # Print the params and write out for ref + fout = open(f"{out_dir_name}/info.txt", "a") params = xgb_clf.get_params() - print("XGBoost Classifier Parameters:\n") + print("\nXGBoost Classifier Parameters:\n") + fout.write("\nXGBoost Classifier Parameters:\n") for param, value in params.items(): - print(f" {param}: {value}") + print(f"\t{param}: {value}") + fout.write(f" {param}: {value}\n") + fout.close() + + + #integral_train_of_sig = np.sum(w_train_of[y_train_of==1]) + #integral_train_of_bkg = np.sum(w_train_of[y_train_of==0]) + #integral_test_of_sig = np.sum(w_test_of[y_test_of==1]) + #integral_test_of_bkg = np.sum(w_test_of[y_test_of==0]) + #w_train_of[y_train_of==1] = w_train_of[y_train_of==1] / integral_train_of_sig + #w_train_of[y_train_of==0] = w_train_of[y_train_of==0] / integral_train_of_bkg + #w_test_of[y_test_of==1] = w_test_of[y_test_of==1] / integral_test_of_sig + #w_test_of[y_test_of==0] = w_test_of[y_test_of==0] / integral_test_of_bkg # Train xgb_clf.fit( @@ -147,9 +188,11 @@ def main(): eval_set=[(X_train_of, y_train_of), (X_test_of, y_test_of)] ) + # Save the model xgb_clf.save_model(f"{out_dir_name}/bdt.json") + ############# Make resutls plots ############# @@ -163,16 +206,30 @@ def main(): p_test_sig = p_test[y_test_of==1] p_test_bkg = p_test[y_test_of==0] - print("test s, len, avg",p_test_sig,len(p_test_sig),sum(p_test_sig)/len(p_test_sig)) - print("test b, len, avg",p_test_bkg,len(p_test_bkg),sum(p_test_bkg)/len(p_test_bkg)) - print("train s, len, avg",p_train_sig,len(p_train_sig),sum(p_train_sig)/len(p_train_sig)) - print("train b, len, avg",p_train_bkg,len(p_train_bkg),sum(p_train_bkg)/len(p_train_bkg)) + w_train_sig = w_train_of[y_train_of==1] + w_train_bkg = w_train_of[y_train_of==0] + w_test_sig = w_test_of[y_test_of==1] + w_test_bkg = w_test_of[y_test_of==0] + + # Print test and train sig and bkg info + print(f"Test and train sig and bkg numbers:") + print(f"Test sig len, avg -> {len(p_test_sig)}, {sum(p_test_sig)/len(p_test_sig)}") + print(f"Test bkg: len, avg -> {len(p_test_bkg)}, {sum(p_test_bkg)/len(p_test_bkg)}") + print(f"Train sig: len, avg -> {len(p_train_sig)}, {sum(p_train_sig)/len(p_train_sig)}") + print(f"Train bkg: len, avg -> {len(p_train_bkg)} , {sum(p_train_bkg)/len(p_train_bkg)}") + fout = open(f"{out_dir_name}/info.txt", "a") + fout.write(f"\nTest and train sig and bkg numbers:\n") + fout.write(f"\tTest sig len, avg -> {len(p_test_sig)}, {sum(p_test_sig)/len(p_test_sig)}\n") + fout.write(f"\tTest bkg: len, avg -> {len(p_test_bkg)}, {sum(p_test_bkg)/len(p_test_bkg)}\n") + fout.write(f"\tTrain sig: len, avg -> {len(p_train_sig)}, {sum(p_train_sig)/len(p_train_sig)}\n") + fout.write(f"\tTrain bkg: len, avg -> {len(p_train_bkg)} , {sum(p_train_bkg)/len(p_train_bkg)}\n") + fout.close() fig, ax = plt.subplots(figsize=(5,5)) - plt.hist(p_train_sig,bins=100,histtype="step",label="train is_sig",density=True) - plt.hist(p_test_sig,bins=100,histtype="step",label="test is_sig",density=True) - plt.hist(p_train_bkg,bins=100,histtype="step",label="train is_bkg",density=True) - plt.hist(p_test_bkg,bins=100,histtype="step",label="test is_bkg",density=True) + plt.hist(p_train_sig,weights=w_train_sig,bins=100,histtype="step",label="train is_sig",density=True) + plt.hist(p_test_sig, weights=w_test_sig, bins=100,histtype="step",label="test is_sig",density=True) + plt.hist(p_train_bkg,weights=w_train_bkg,bins=100,histtype="step",label="train is_bkg",density=True) + plt.hist(p_test_bkg, weights=w_test_bkg, bins=100,histtype="step",label="test is_bkg",density=True) plt.legend() ax.set_xlim(0,1) @@ -189,12 +246,15 @@ def main(): # preparing evaluation metric plots results = xgb_clf.evals_result() epochs = len(results['validation_0']['logloss']) + #epochs = len(results['validation_0']['mlogloss']) x_axis = range(0, epochs) # xgboost 'logloss' plot fig, ax = plt.subplots(figsize=(5,5)) ax.plot(x_axis, results['validation_0']['logloss'], label='Train') ax.plot(x_axis, results['validation_1']['logloss'], label='Test') + #ax.plot(x_axis, results['validation_0']['mlogloss'], label='Train') + #ax.plot(x_axis, results['validation_1']['mlogloss'], label='Test') ax.legend() plt.ylabel('logloss') plt.title('GridSearchCV XGBoost logloss') @@ -209,6 +269,8 @@ def main(): fig, ax = plt.subplots(figsize=(5,5)) ax.plot(x_axis, results['validation_0']['error'], label='Train') ax.plot(x_axis, results['validation_1']['error'], label='Test') + #ax.plot(x_axis, results['validation_0']['merror'], label='Train') + #ax.plot(x_axis, results['validation_1']['merror'], label='Test') ax.legend() plt.ylabel('error') plt.title('GridSearchCV XGBoost error')