diff --git a/src/BTVNanoCommissioning/utils/correction.py b/src/BTVNanoCommissioning/utils/correction.py index 594e66ca..78b9f99f 100644 --- a/src/BTVNanoCommissioning/utils/correction.py +++ b/src/BTVNanoCommissioning/utils/correction.py @@ -992,12 +992,14 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): if "correctionlib" in str(type(correct_map["EGM"])): ## Reco SF -splitted pT if "Reco" in sf: - if "Summer22" not in correct_map["campaign"]: - ele_pt = np.where(ele.pt < 20.0, 20.0, ele.pt) - ele_pt_low = np.where(ele.pt >= 20.0, 19.9, ele.pt) + ## phi is used in Summer23 + ele_pt = np.clip(ele.pt, 20.1, 74.9) + ele_pt_low = np.where(ele.pt >= 20.0, 19.9, ele.pt) + ele_pt_high = np.clip(ele.pt, 75.0, 500.0) + if "Summer23" in correct_map["campaign"]: sfs_low = np.where( - (~mask) & (~masknone), - correct_map["EGM"][correct_map["EGM_cfg"][sf]].evaluate( + (ele.pt <= 20.0) & (~masknone), + correct_map["EGM"][sf.split(" ")[2]].evaluate( sf.split(" ")[1], "sf", "RecoBelow20", @@ -1007,18 +1009,36 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): ), 1.0, ) - sfs = np.where( - mask & (~masknone), + sfs_high = np.where( + (ele.pt > 75.0) & (~masknone), correct_map["EGM"][sf.split(" ")[2]].evaluate( - sf.split(" ")[1], "sf", "RecoAbove20", ele_eta, ele_pt + sf.split(" ")[1], + "sf", + "RecoAbove75", + ele_eta, + ele_pt_high, + ele.phi, ), sfs_low, ) + sfs = np.where( + (ele.pt > 20.0) & (ele.pt <= 75.0) & (~masknone), + correct_map["EGM"][sf.split(" ")[2]].evaluate( + sf.split(" ")[1], + "sf", + "Reco20to75", + ele_eta, + ele_pt, + ele.phi, + ), + sfs_high, + ) + sfs = np.where(masknone, 1.0, sfs) if syst: sfs_up_low = np.where( - ~mask & ~masknone, + (ele.pt <= 20.0) & ~masknone, correct_map["EGM"][sf.split(" ")[2]].evaluate( sf.split(" ")[1], "sfup", @@ -1030,7 +1050,7 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): 0.0, ) sfs_down_low = np.where( - ~mask & ~masknone, + (ele.pt <= 20.0) & ~masknone, correct_map["EGM"][sf.split(" ")[2]].evaluate( sf.split(" ")[1], "sfdown", @@ -1041,37 +1061,59 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): ), 0.0, ) + sfs_up_high = np.where( + (ele.pt > 20.0) & (ele.pt <= 75.0) & ~masknone, + correct_map["EGM"][sf.split(" ")[2]].evaluate( + sf.split(" ")[1], + "sfup", + "RecoAbove75", + ele_eta, + ele_pt_high, + ele.phi, + ), + sfs_up_low, + ) + sfs_down_high = np.where( + (ele.pt > 20.0) & (ele.pt <= 75.0) & ~masknone, + correct_map["EGM"][sf.split(" ")[2]].evaluate( + sf.split(" ")[1], + "sfdown", + "RecoAbove75", + ele_eta, + ele_pt_high, + ele.phi, + ), + sfs_down_low, + ) sfs_up = np.where( - mask & ~masknone, + (ele.pt > 20.0) & (ele.pt <= 75.0) & ~masknone, correct_map["EGM"][sf.split(" ")[2]].evaluate( sf.split(" ")[1], "sfup", - "RecoAbove20", + "Reco20to75", ele_eta, ele_pt, ele.phi, ), - sfs_up_low, + sfs_up_high, ) sfs_down = np.where( - mask & ~masknone, + (ele.pt > 20.0) & (ele.pt <= 75.0) & ~masknone, correct_map["EGM"][sf.split(" ")[2]].evaluate( sf.split(" ")[1], "sfdown", - "RecoAbove20", + "Reco20to75", ele_eta, ele_pt, ele.phi, ), - sfs_down_low, + sfs_down_high, ) sfs_up, sfs_down = np.where( masknone, 1.0, sfs_up ), np.where(masknone, 1.0, sfs_down) + else: - ele_pt = np.clip(ele.pt, 20.1, 74.9) - ele_pt_low = np.where(ele.pt >= 20.0, 19.9, ele.pt) - ele_pt_high = np.clip(ele.pt, 75.0, 500.0) sfs_low = np.where( (ele.pt <= 20.0) & (~masknone), @@ -1081,7 +1123,6 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): "RecoBelow20", ele_eta, ele_pt_low, - ele.phi, ), 1.0, ) @@ -1093,7 +1134,6 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): "RecoAbove75", ele_eta, ele_pt_high, - ele.phi, ), sfs_low, ) @@ -1116,7 +1156,6 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): "RecoBelow20", ele_eta, ele_pt_low, - ele.phi, ), 0.0, ) @@ -1128,7 +1167,6 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): "RecoBelow20", ele_eta, ele_pt_low, - ele.phi, ), 0.0, ) @@ -1151,7 +1189,6 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): "RecoAbove75", ele_eta, ele_pt_high, - ele.phi, ), sfs_down_low, ) @@ -1163,7 +1200,6 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): "Reco20to75", ele_eta, ele_pt, - ele.phi, ), sfs_up_high, ) @@ -1175,7 +1211,6 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): "Reco20to75", ele_eta, ele_pt, - ele.phi, ), sfs_down_high, ) @@ -1185,43 +1220,82 @@ def eleSFs(ele, correct_map, weights, syst=True, isHLT=False): ), np.where(masknone, 1.0, sfs_down) ## Other files else: - sfs = np.where( - masknone, - 1.0, - correct_map["EGM"][sf.split(" ")[2]].evaluate( - sf.split(" ")[1], - "sf", - correct_map["EGM_cfg"][sf], - ele_eta, - ele_pt, - ele.phi, - ), - ) - - if syst: - sfs_up = np.where( + if "Summer23" in correct_map["campaign"]: + sfs = np.where( masknone, 1.0, correct_map["EGM"][sf.split(" ")[2]].evaluate( sf.split(" ")[1], - "sfup", + "sf", correct_map["EGM_cfg"][sf], ele_eta, ele_pt, ele.phi, ), ) - sfs_down = np.where( + + if syst: + sfs_up = np.where( + masknone, + 1.0, + correct_map["EGM"][sf.split(" ")[2]].evaluate( + sf.split(" ")[1], + "sfup", + correct_map["EGM_cfg"][sf], + ele_eta, + ele_pt, + ele.phi, + ), + ) + sfs_down = np.where( + masknone, + 1.0, + correct_map["EGM"][sf.split(" ")[2]].evaluate( + sf.split(" ")[1], + "sfdown", + correct_map["EGM_cfg"][sf], + ele_eta, + ele_pt, + ele.phi, + ), + ) + else: + sfs = np.where( masknone, 1.0, correct_map["EGM"][sf.split(" ")[2]].evaluate( sf.split(" ")[1], - "sfdown", + "sf", correct_map["EGM_cfg"][sf], ele_eta, ele_pt, ), ) + + if syst: + sfs_up = np.where( + masknone, + 1.0, + correct_map["EGM"][sf.split(" ")[2]].evaluate( + sf.split(" ")[1], + "sfup", + correct_map["EGM_cfg"][sf], + ele_eta, + ele_pt, + ), + ) + sfs_down = np.where( + masknone, + 1.0, + correct_map["EGM"][sf.split(" ")[2]].evaluate( + sf.split(" ")[1], + "sfdown", + correct_map["EGM_cfg"][sf], + ele_eta, + ele_pt, + ), + ) + else: if "ele_Trig" in sf: sfs = np.where( @@ -1681,9 +1755,9 @@ def weight_manager(pruned_ev, SF_map, isSyst): syst_wei, ) if "MUO" in SF_map.keys() and "SelMuon" in pruned_ev.fields: - muSFs(pruned_ev.Muon, SF_map, weights, syst_wei, False) + muSFs(pruned_ev.SelMuon, SF_map, weights, syst_wei, False) if "EGM" in SF_map.keys() and "SelElectron" in pruned_ev.fields: - eleSFs(pruned_ev.Electron, SF_map, weights, syst_wei, False) + eleSFs(pruned_ev.SelElectron, SF_map, weights, syst_wei, False) if "BTV" in SF_map.keys() and "SelJet" in pruned_ev.fields: btagSFs(pruned_ev.SelJet, SF_map, weights, "DeepJetC", syst_wei) btagSFs(pruned_ev.SelJet, SF_map, weights, "DeepJetB", syst_wei) diff --git a/src/BTVNanoCommissioning/utils/histogrammer.py b/src/BTVNanoCommissioning/utils/histogrammer.py index 4600b4e9..eb3b679f 100644 --- a/src/BTVNanoCommissioning/utils/histogrammer.py +++ b/src/BTVNanoCommissioning/utils/histogrammer.py @@ -11,7 +11,7 @@ def histogrammer(events, workflow, year="2022", campaign="Summer22"): _hist_dict = {} - ## Common variables + ## Common axis flav_axis = Hist.axis.IntCategory([0, 1, 4, 5, 6], name="flav", label="Genflavour") syst_axis = Hist.axis.StrCategory([], name="syst", growth=True) pt_axis = Hist.axis.Regular(60, 0, 300, name="pt", label=" $p_{T}$ [GeV]") @@ -36,6 +36,7 @@ def histogrammer(events, workflow, year="2022", campaign="Summer22"): ptratio_axis = Hist.axis.Regular(50, 0, 1, name="ratio", label="ratio") n_axis = Hist.axis.Integer(0, 10, name="n", label="N obj") osss_axis = Hist.axis.IntCategory([1, -1], name="osss", label="OS(+)/SS(-)") + ## create histograms for each workflow ### Workflow specific if "example" == workflow: obj_list = [ @@ -650,6 +651,7 @@ def histogrammer(events, workflow, year="2022", campaign="Summer22"): return _hist_dict +# Filled common histogram def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): exclude_btv = [ "DeepCSVC", @@ -678,16 +680,19 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): genflavor = ak.zeros_like(pruned_ev.SelJet.pt, dtype=int) if "MuonJet" in pruned_ev.fields: smflav = ak.zeros_like(pruned_ev.MuonJet.pt, dtype=int) - + # Loop over the systematic variations for syst in systematics: if isSyst == False and syst != "nominal": break + # weight modifications for systematics weight = ( weights.weight() if syst == "nominal" or syst not in list(weights.variations) else weights.weight(modifier=syst) ) + # Loop over the histograms for histname, h in output.items(): + # tagger score histograms if ( "Deep" in histname and "btag" not in histname @@ -705,6 +710,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): )[0] ), ) + # PFcands histograms elif ( "PFCands" in pruned_ev.fields and "PFCands" in histname @@ -740,7 +746,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): )[0] ), ) - + # leading lepton histograms elif ( "hl_" in histname and "hl" in pruned_ev.fields @@ -752,6 +758,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): flatten(pruned_ev.hl[histname.replace("hl_", "")]), weight=weight, ) + # subleading lepton histograms elif ( "sl_" in histname and "sl" in pruned_ev.fields @@ -762,6 +769,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): flatten(pruned_ev.sl[histname.replace("sl_", "")]), weight=weight, ) + # Selected electron histograms elif ( "ele_" in histname and histname.replace("ele_", "") in pruned_ev.SelElectron.fields @@ -773,6 +781,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): flatten(pruned_ev.SelElectron[histname.replace("ele_", "")]), weight=weight, ) + # Selected muon histograms elif ( "mu_" in histname and histname.replace("mu_", "") in pruned_ev.SelMuon.fields @@ -783,6 +792,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): flatten(pruned_ev.SelMuon[histname.replace("mu_", "")]), weight=weight, ) + # Negatively charged lepton histograms-in DY workflow elif ( "negl_" in histname and histname.replace("negl_", "") in pruned_ev.negl.fields @@ -792,6 +802,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): flatten(pruned_ev.negl[histname.replace("negl_", "")]), weight=weight, ) + # Posively charged lepton histograms-in DY workflow elif ( "posl_" in histname and histname.replace("posl_", "") in pruned_ev.posl.fields @@ -801,6 +812,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): flatten(pruned_ev.posl[histname.replace("posl_", "")]), weight=weight, ) + # Soft muon histograms elif "soft_l" in histname and not "ptratio" in histname: h.fill( syst, @@ -810,7 +822,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): ) elif "njet" == histname: output["njet"].fill(syst, pruned_ev.njet, weight=weight) - + # Jet kinmeatics & deltaR between jet and lepton elif ( "jet" in histname and "posl" not in histname and "negl" not in histname ): @@ -847,6 +859,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): ) # filled discriminants elif "btag" in histname or "PNet" in histname: + # Events with muon jet if "MuonJet" in pruned_ev.fields: flavs, seljets = smflav, pruned_ev.MuonJet nj = 1 @@ -873,7 +886,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): discr=seljet[histname.replace(f"_{i}", "")], weight=weight, ) - + # pT ratio if "hl" in pruned_ev.fields: output["hl_ptratio"].fill( syst, @@ -898,7 +911,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): output["dr_negljet"].fill( syst, genflavor, pruned_ev.negl.delta_r(pruned_ev.SelJet), weight=weight ) - + # Muon enriched jet histograms if "MuonJet" in pruned_ev.fields: if ( "hl" not in pruned_ev.fields @@ -952,7 +965,7 @@ def histo_writter(pruned_ev, output, weights, systematics, isSyst, SF_map): output["dr_hmusmu"].fill( syst, pruned_ev.SelMuon.delta_r(pruned_ev.SoftMuon), weight=weight ) - + # dilepton system histograms: DY workflow if "dilep" in pruned_ev.fields: output["dilep_pt"].fill(syst, flatten(pruned_ev.dilep.pt), weight=weight) output["dilep_pt"].fill(syst, flatten(pruned_ev.dilep.eta), weight=weight)