Merge pull request #445 from TopEFT/run3id

Object selection organization and run 3D ID

analysis/topeft_run2/analysis_processor.py

@@ -229,27 +229,19 @@ def process(self, events):
 
 
         if is_run3:
-            elePresId = ele.mvaNoIso_WP90
-            ele_pres_id_tag = "wp90noiso"
-            eleFOId = ele.mvaNoIso_WP80
-            eleMVATTH = ele.mvaTTH
-            muMVATTH = mu.mvaTTH
+            leptonSelection = te_os.run3leptonselection()
             jetsRho = events.Rho["fixedGridRhoFastjetAll"]
         elif is_run2:
-            elePresId = ele.mvaFall17V2noIso_WPL
-            ele_pres_id_tag = "wpLnoiso"
-            eleFOId = ele.mvaFall17V2noIso_WP90
-            eleMVATTH = ele.mvaTTHUL
-            muMVATTH = mu.mvaTTHUL
+            leptonSelection = te_os.run2leptonselection()
             jetsRho = events.fixedGridRhoFastjetAll
 
         # An array of lenght events that is just 1 for each event
         # Probably there's a better way to do this, but we use this method elsewhere so I guess why not..
         events.nom = ak.ones_like(events.MET.pt)
 
         ele["idEmu"] = te_os.ttH_idEmu_cuts_E3(ele.hoe, ele.eta, ele.deltaEtaSC, ele.eInvMinusPInv, ele.sieie)
-        ele["conept"] = te_os.coneptElec(ele.pt, eleMVATTH, ele.jetRelIso)
-        mu["conept"] = te_os.coneptMuon(mu.pt, muMVATTH, mu.jetRelIso, mu.mediumId)
+        ele["conept"] = leptonSelection.coneptElec(ele)
+        mu["conept"] = leptonSelection.coneptMuon(mu)
         ele["btagDeepFlavB"] = ak.fill_none(ele.matched_jet.btagDeepFlavB, -99)
         mu["btagDeepFlavB"] = ak.fill_none(mu.matched_jet.btagDeepFlavB, -99)
         if not isData:
@@ -286,18 +278,18 @@ def process(self, events):
 
         ################### Electron selection ####################
 
-        ele["isPres"] = te_os.isPresElec(ele.pt, ele.eta, ele.dxy, ele.dz, ele.miniPFRelIso_all, ele.sip3d, elePresId)
-        ele["isLooseE"] = te_os.isLooseElec(ele.miniPFRelIso_all,ele.sip3d,ele.lostHits)
-        ele["isFO"] = te_os.isFOElec(ele.pt, ele.conept, ele.btagDeepFlavB, ele.idEmu, ele.convVeto, ele.lostHits, eleMVATTH, ele.jetRelIso, eleFOId, year)
-        ele["isTightLep"] = te_os.tightSelElec(ele.isFO, eleMVATTH)
+        ele["isPres"] = leptonSelection.isPresElec(ele)
+        ele["isLooseE"] = leptonSelection.isLooseElec(ele)
+        ele["isFO"] = leptonSelection.isFOElec(ele, year)
+        ele["isTightLep"] = leptonSelection.tightSelElec(ele)
 
         ################### Muon selection ####################
 
         mu["pt"] = ApplyRochesterCorrections(year, mu, isData) # Run3 ready
-        mu["isPres"] = te_os.isPresMuon(mu.dxy, mu.dz, mu.sip3d, mu.eta, mu.pt, mu.miniPFRelIso_all)
-        mu["isLooseM"] = te_os.isLooseMuon(mu.miniPFRelIso_all,mu.sip3d,mu.looseId)
-        mu["isFO"] = te_os.isFOMuon(mu.pt, mu.conept, mu.btagDeepFlavB, muMVATTH, mu.jetRelIso, year)
-        mu["isTightLep"]= te_os.tightSelMuon(mu.isFO, mu.mediumId, muMVATTH)
+        mu["isPres"] = leptonSelection.isPresMuon(mu)
+        mu["isLooseM"] = leptonSelection.isLooseMuon(mu)
+        mu["isFO"] = leptonSelection.isFOMuon(mu, year)
+        mu["isTightLep"]= leptonSelection.tightSelMuon(mu)
 
         ################### Loose selection ####################
 
@@ -322,7 +314,7 @@ def process(self, events):
         e_fo = ele[ele.isPres & ele.isLooseE & ele.isFO]
 
         # Attach the lepton SFs to the electron and muons collections
-        AttachElectronSF(e_fo, year=year, looseWP=ele_pres_id_tag) #Run3 ready
+        AttachElectronSF(e_fo, year=year, looseWP="wpLnoiso" if is_run3 else "none") #Run3 ready
         AttachMuonSF(m_fo,year=year) #Run3 ready
 
         # Attach per lepton fake rates

topeft/modules/corrections.py

@@ -1168,33 +1168,7 @@ def AttachElectronSF(electrons, year, looseWP=None):
     reco_up = ak.unflatten(reco_up_flat, ak.num(pt))
     reco_do = ak.unflatten(reco_do_flat, ak.num(pt))
 
-    if is_run3:
-        loose_sf_flat = None
-        loose_up_flat = None
-        loose_do_flat = None
-        pt_mask = ak.flatten((pt >= 10))
-        egm_args = [eta_flat, pt_bin_flat]
-        if "2023" in year:
-            egm_args.append(phi_flat)
-        loose_sf_flat = ak.where(
-            ~pt_mask,
-            1,
-            ceval[egm_tag].evaluate(egm_year, "sf", looseWP, *egm_args)
-        )
-        loose_up_flat = ak.where(
-            ~pt_mask,
-            1,
-            ceval[egm_tag].evaluate(egm_year, "sfup", looseWP, *egm_args)
-        )
-        loose_do_flat = ak.where(
-            ~pt_mask,
-            1,
-            ceval[egm_tag].evaluate(egm_year, "sfdown", looseWP, *egm_args)
-        )
-        loose_sf = ak.unflatten(loose_sf_flat, ak.num(pt))
-        loose_up = ak.unflatten(loose_up_flat, ak.num(pt))
-        loose_do = ak.unflatten(loose_do_flat, ak.num(pt))
-    else:
+    if not is_run3: # run 3 dont need loose for the id because we dont have ID working point (to be checked!!)
         loose_sf  = SFevaluator['ElecLooseSF_{year}'.format(year=year)](np.abs(eta),pt)
         loose_err = SFevaluator['ElecLooseSF_{year}_er'.format(year=year)](np.abs(eta),pt)
         loose_up = loose_sf + loose_err

topeft/modules/object_selection.py

@@ -72,100 +72,147 @@ def smoothBFlav(jetpt,ptmin,ptmax,year,scale_loose=1.0):
     x = np.minimum(np.maximum(0, jetpt - ptmin)/(ptmax-ptmin), 1.0)
     return x*wploose*scale_loose + (1-x)*wpmedium
 
-def coneptElec(pt, mvaTTHUL, jetRelIso):
-    conePt = (0.90 * pt * (1 + jetRelIso))
-    return ak.where((mvaTTHUL>get_tc_param("mva_TTH_e_cut")),pt,conePt)
-
-def coneptMuon(pt, mvaTTHUL, jetRelIso, mediumId):
-    conePt = (0.90 * pt * (1 + jetRelIso))
-    return ak.where(((mvaTTHUL>get_tc_param("mva_TTH_m_cut"))&(mediumId>0)),pt,conePt)
-
-def isPresElec(pt, eta, dxy, dz, miniIso, sip3D, eleId):
-    pt_mask    = (pt       > get_te_param("pres_e_pt_cut"))
-    eta_mask   = (abs(eta) < get_te_param("eta_e_cut"))
-    dxy_mask   = (abs(dxy) < get_te_param("dxy_cut"))
-    dz_mask    = (abs(dz)  < get_te_param("dz_cut"))
-    iso_mask   = (miniIso  < get_te_param("iso_cut"))
-    sip3d_mask = (sip3D    < get_te_param("sip3d_cut"))
-    return (pt_mask & eta_mask & dxy_mask & dz_mask & iso_mask & sip3d_mask & eleId)
-
-def isPresMuon(dxy, dz, sip3D, eta, pt, miniRelIso):
-    pt_mask    = (pt         > get_te_param("pres_m_pt_cut"))
-    eta_mask   = (abs(eta)   < get_te_param("eta_m_cut"))
-    dxy_mask   = (abs(dxy)   < get_te_param("dxy_cut"))
-    dz_mask    = (abs(dz)    < get_te_param("dz_cut"))
-    iso_mask   = (miniRelIso < get_te_param("iso_cut"))
-    sip3d_mask = (sip3D      < get_te_param("sip3d_cut"))
-    return (pt_mask & eta_mask & dxy_mask & dz_mask & iso_mask & sip3d_mask)
-
-def isLooseElec(miniPFRelIso_all,sip3d,lostHits):
-    return (miniPFRelIso_all<get_te_param("iso_cut")) & (sip3d<get_te_param("sip3d_cut")) & (lostHits<=1)
-
-def isLooseMuon(miniPFRelIso_all,sip3d,looseId):
-    return (miniPFRelIso_all<get_te_param("iso_cut")) & (sip3d<get_te_param("sip3d_cut")) & (looseId)
-
-def isFOElec(pt, conept, jetBTagDeepFlav, ttH_idEmu_cuts_E3, convVeto, lostHits, mvaTTHUL, jetRelIso, mvaFall17V2noIso_WP90, year):
-
-    # Get the btag cut for the year
-    if (year == "2016"):
-        bTagCut = get_tc_param("btag_wp_medium_UL16")
-    elif (year == "2016APV"):
-        bTagCut = get_tc_param("btag_wp_medium_UL16APV")
-    elif (year == "2017"):
-        bTagCut = get_tc_param("btag_wp_medium_UL17")
-    elif (year == "2018"):
-        bTagCut = get_tc_param("btag_wp_medium_UL18")
-    elif (year == "2022"):
-        bTagCut = get_tc_param("btag_wp_medium_2022")
-    elif (year == "2022EE"):
-        bTagCut = get_tc_param("btag_wp_medium_2022EE")
-    elif (year == "2023"):
-        bTagCut = get_tc_param("btag_wp_medium_2023")
-    elif (year == "2023BPix"):
-        bTagCut = get_tc_param("btag_wp_medium_2023BPix")
-
-    else:
-        raise Exception(f"Error: Unknown year \"{year}\". Exiting...")
-
-    btabReq    = (jetBTagDeepFlav<bTagCut)
-    ptReq      = (conept>get_te_param("fo_pt_cut"))
-    qualityReq = (ttH_idEmu_cuts_E3 & convVeto & (lostHits==0))
-    mvaReq     = ((mvaTTHUL>get_tc_param("mva_TTH_e_cut")) | ((mvaFall17V2noIso_WP90) & (jetBTagDeepFlav<smoothBFlav(0.9*pt*(1+jetRelIso),20,45,year)) & (jetRelIso < get_te_param("fo_e_jetRelIso_cut"))))
-
-    return ptReq & btabReq & qualityReq & mvaReq
-
-def isFOMuon(pt, conept, jetBTagDeepFlav, mvaTTHUL, jetRelIso, year):
+def smoothSip3D(jetpt, sipmin, sipmax, ptmin, ptmax):
+    x = np.minimum(np.maximum(0, jetpt - ptmin)/(ptmax-ptmin), 1.0)
+    return x*sipmax+(1-x)*sipmin
 
+def get_medium_btag_foryear(year):
     # Get the btag cut for the year
     if (year == "2016"):
-        bTagCut = get_tc_param("btag_wp_medium_UL16")
+        return get_tc_param("btag_wp_medium_UL16")
     elif (year == "2016APV"):
-        bTagCut = get_tc_param("btag_wp_medium_UL16APV")
+        return get_tc_param("btag_wp_medium_UL16APV")
     elif (year == "2017"):
-        bTagCut = get_tc_param("btag_wp_medium_UL17")
+        return get_tc_param("btag_wp_medium_UL17")
     elif (year == "2018"):
-        bTagCut = get_tc_param("btag_wp_medium_UL18")
+        return get_tc_param("btag_wp_medium_UL18")
     elif (year == "2022"):
-        bTagCut = get_tc_param("btag_wp_medium_2022")
+        return get_tc_param("btag_wp_medium_2022")
     elif (year == "2022EE"):
-        bTagCut = get_tc_param("btag_wp_medium_2022EE")
+        return get_tc_param("btag_wp_medium_2022EE")
     elif (year == "2023"):
-        bTagCut = get_tc_param("btag_wp_medium_2023")
+        return get_tc_param("btag_wp_medium_2023")
     elif (year == "2023BPix"):
-        bTagCut = get_tc_param("btag_wp_medium_2023BPix")
+        return get_tc_param("btag_wp_medium_2023BPix")
     else:
         raise Exception(f"Error: Unknown year \"{year}\". Exiting...")
 
-    btagReq = (jetBTagDeepFlav<bTagCut)
-    ptReq   = (conept>get_te_param("fo_pt_cut"))
-    mvaReq  = ((mvaTTHUL>get_tc_param("mva_TTH_m_cut")) | ((jetBTagDeepFlav<smoothBFlav(0.9*pt*(1+jetRelIso),20,45,year)) & (jetRelIso < get_te_param("fo_m_jetRelIso_cut"))))
-    return ptReq & btagReq & mvaReq
-
-def tightSelElec(clean_and_FO_selection_TTH, mvaTTHUL):
-    return (clean_and_FO_selection_TTH) & (mvaTTHUL > get_tc_param("mva_TTH_e_cut"))
-
-def tightSelMuon(clean_and_FO_selection_TTH, mediumId, mvaTTHUL):
-    return (clean_and_FO_selection_TTH) & (mediumId>0) & (mvaTTHUL > get_tc_param("mva_TTH_m_cut"))
+class run2leptonselection:
+
+    def __init__(self):
+        pass
+
+    def coneptElec(self, ele):
+        conePt = (0.90 * ele.pt * (1 + ele.jetRelIso))
+        return ak.where((ele.mvaTTHUL>get_te_param("mva_TTH_e_cut")),ele.pt,conePt)
+
+    def coneptMuon(self, muo):
+        conePt = (0.90 * muo.pt * (1 + muo.jetRelIso))
+        return ak.where(((muo.mvaTTHUL>get_te_param("mva_TTH_m_cut"))&(muo.mediumId>0)),muo.pt,conePt)
+
+    def isPresElec(self, ele):
+        pt_mask    = (ele.pt       > get_te_param("pres_e_pt_cut"))
+        eta_mask   = (abs(ele.eta) < get_te_param("eta_e_cut"))
+        dxy_mask   = (abs(ele.dxy) < get_te_param("dxy_cut"))
+        dz_mask    = (abs(ele.dz)  < get_te_param("dz_cut"))
+        iso_mask   = (ele.miniPFRelIso_all  < get_te_param("iso_cut"))
+        sip3d_mask = (ele.sip3d    < get_te_param("sip3d_cut"))
+        return (pt_mask & eta_mask & dxy_mask & dz_mask & iso_mask & sip3d_mask & ele.mvaFall17V2noIso_WPL)
+
+    def isPresMuon(self, muon):
+        pt_mask    = (muon.pt         > get_te_param("pres_m_pt_cut"))
+        eta_mask   = (abs(muon.eta)   < get_te_param("eta_m_cut"))
+        dxy_mask   = (abs(muon.dxy)   < get_te_param("dxy_cut"))
+        dz_mask    = (abs(muon.dz)    < get_te_param("dz_cut"))
+        iso_mask   = (muon.miniPFRelIso_all < get_te_param("iso_cut"))
+        sip3d_mask = (muon.sip3d      < get_te_param("sip3d_cut"))
+        return (pt_mask & eta_mask & dxy_mask & dz_mask & iso_mask & sip3d_mask)
+
+    def isLooseElec(self, ele):
+        return (ele.miniPFRelIso_all<get_te_param("iso_cut")) & (ele.sip3d<get_te_param("sip3d_cut")) & (ele.lostHits<=1)
+
+    def isLooseMuon(self, muon):
+        return (muon.miniPFRelIso_all<get_te_param("iso_cut")) & (muon.sip3d<get_te_param("sip3d_cut")) & (muon.looseId)
+
+    def isFOElec(self, ele, year):
+        bTagCut=get_medium_btag_foryear(year)
+        btabReq    = (ele.jetBTagDeepFlav<bTagCut)
+        ptReq      = (ele.conept>get_te_param("fo_pt_cut"))
+        qualityReq = (ele.idEmu & ele.convVeto & (ele.lostHits==0))
+        mvaReq     = ((ele.mvaTTHUL>get_te_param("mva_TTH_e_cut")) | ((ele.mvaFall17V2noIso_WP90) & (ele.jetBTagDeepFlav<smoothBFlav(0.9*ele.pt*(1+ele.jetRelIso),20,45,year)) & (ele.jetRelIso < get_te_param("fo_e_jetRelIso_cut"))))
+
+        return ptReq & btabReq & qualityReq & mvaReq
+
+    def isFOMuon(self, muo, year):
+        bTagCut=get_medium_btag_foryear(year)
+        btagReq = (muo.jetBTagDeepFlav<bTagCut)
+        ptReq   = (muo.conept>get_te_param("fo_pt_cut"))
+        mvaReq  = ((muo.mvaTTHUL>get_te_param("mva_TTH_m_cut")) | ((muo.jetBTagDeepFlav<smoothBFlav(0.9*muo.pt*(1+muo.jetRelIso),20,45,year)) & (muo.jetRelIso < get_te_param("fo_m_jetRelIso_cut"))))
+        return ptReq & btagReq & mvaReq
+
+    def tightSelElec(self, ele):
+        return (ele.isFO) & (ele.mvaTTHUL > get_te_param("mva_TTH_e_cut"))
+
+    def tightSelMuon(self, muo):
+        return (muo.isFO) & (muo.mediumId>0) & (muo.mvaTTHUL > get_te_param("mva_TTH_m_cut"))
+
+class run3leptonselection:
+
+    def __init__(self):
+        pass
+
+    def coneptElec(self, ele):
+        conePt = (0.90 * ele.pt * (1 + ele.jetRelIso))
+        return ak.where((ele.mvaTTH_Run3>get_te_param("mva_TTH_e_cut_run3")),ele.pt,conePt)
+
+    def coneptMuon(self, muo):
+        conePt = (0.90 * muo.pt * (1 + muo.jetRelIso))
+        return ak.where(((muo.mvaTTH_Run3>get_te_param("mva_TTH_m_cut_run3"))&(muo.mediumId>0)),muo.pt,conePt)
+
+    def isPresElec(self, ele):
+        pt_mask    = (ele.pt       > get_te_param("pres_e_pt_cut"))
+        eta_mask   = (abs(ele.eta) < get_te_param("eta_e_cut"))
+        dxy_mask   = (abs(ele.dxy) < get_te_param("dxy_cut"))
+        dz_mask    = (abs(ele.dz)  < get_te_param("dz_cut"))
+        iso_mask   = (ele.miniPFRelIso_all  < get_te_param("iso_cut"))
+        sip3d_mask = (ele.sip3d    < get_te_param("sip3d_cut"))
+        ecal_crack_mask = (((abs(ele.etaSC) < 1.4442) | (abs(ele.etaSC) > 1.566)))
+        return (pt_mask & eta_mask & dxy_mask & dz_mask & iso_mask & sip3d_mask & ecal_crack_mask)
+
+    def isPresMuon(self, muon):
+        pt_mask    = (muon.pt         > get_te_param("pres_m_pt_cut"))
+        eta_mask   = (abs(muon.eta)   < get_te_param("eta_m_cut"))
+        dxy_mask   = (abs(muon.dxy)   < get_te_param("dxy_cut"))
+        dz_mask    = (abs(muon.dz)    < get_te_param("dz_cut"))
+        iso_mask   = (muon.miniPFRelIso_all < get_te_param("iso_cut"))
+        sip3d_mask = (muon.sip3d      < get_te_param("sip3d_cut"))
+        return (pt_mask & eta_mask & dxy_mask & dz_mask & iso_mask & sip3d_mask)
+
+    def isLooseElec(self, ele):
+        return (ele.miniPFRelIso_all<get_te_param("iso_cut")) & (ele.sip3d<get_te_param("sip3d_cut")) & (ele.lostHits<=1)
+
+    def isLooseMuon(self, muon):
+        return (muon.miniPFRelIso_all<get_te_param("iso_cut")) & (muon.sip3d<get_te_param("sip3d_cut")) & (muon.mediumId)
+
+    def isFOElec(self, ele, year):
+        bTagCut=get_medium_btag_foryear(year)
+        btabReq    = (ele.jetBTagDeepFlav<bTagCut)
+        ptReq      = (ele.conept>get_te_param("fo_pt_cut"))
+        qualityReq = (ele.idEmu & ele.convVeto & (ele.lostHits==0))
+        mvaReq     = ((ele.mvaTTH_Run3>get_te_param("mva_TTH_e_cut_run3")) | ((ele.mvaIso > get_te_param("fo_e_mvaiso_cut_run3"))  & (ele.jetRelIso < get_te_param("fo_e_jetRelIso_cut"))))
+
+    def isFOMuon(self, muo, year):
+        bTagCut=get_medium_btag_foryear(year)
+        btagReq = (muo.jetBTagDeepFlav<bTagCut)
+        ptReq   = (muo.conept>get_te_param("fo_pt_cut"))
+        mvaReq  = ((muo.mvaTTH_Run3>get_te_param("mva_TTH_m_cut_run3")) | ((muo.jetBTagDeepFlav<smoothBFlav(0.9*muo.pt*(1+muo.jetRelIso),20,45,year)) & (muo.jetRelIso < get_te_param("fo_m_jetRelIso_cut")) & (muo.sip3d < smoothSip3D(0.9*muo.pt*(1+muo.jetRelIso),2.5,8.,15,45))))
+        return ptReq & btagReq & mvaReq
+
+    def tightSelElec(self, ele):
+        return (ele.isFO) & (ele.mvaTTH_Run3 > get_te_param("mva_TTH_e_cut_run3"))
+
+    def tightSelMuon(self, muo):
+        return (muo.isFO) & (muo.mediumId>0) & (muo.mvaTTH_Run3 > get_te_param("mva_TTH_m_cut_run3"))
 
 def isClean(obj_A, obj_B, drmin=0.4):
     objB_near, objB_DR = obj_A.nearest(obj_B, return_metric=True)

topeft/params/params.json

@@ -8,11 +8,18 @@
     "pres_e_pt_cut" : 7.0,
     "pres_m_pt_cut" : 5.0,
 
+    "mva_TTH_e_cut" : 0.90,
+    "mva_TTH_m_cut" : 0.85,
+    "mva_TTH_e_cut_run3" : 0.90,
+    "mva_TTH_m_cut_run3" : 0.64,
+
     "iso_cut" : 0.4,
 
     "fo_m_jetRelIso_cut" : 0.5,
     "fo_e_jetRelIso_cut" : 1.0,
 
+    "fo_e_mvaiso_cut_run3" : 0.6,
+
     "jet_id_cut" : 0,
 
     "dxy_cut" : 0.05,