diff --git a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py index f126d1da55eff..ac80309bfeb94 100644 --- a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py +++ b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py @@ -789,6 +789,7 @@ Plot1D('charge', 'charge', 3, -1.5, 1.5, 'electric charge'), Plot1D('chargedIso', 'chargedIso', 20, 0, 200, 'charged isolation'), Plot1D('decayMode', 'decayMode', 12, -0.5, 11.5, 'decayMode()'), + Plot1D('decayModePNet', 'decayModePNet', 13, -1.5, 11.5, 'decay mode of the highest PNet tau score'), Plot1D('dxy', 'dxy', 20, -1000, 1000, 'd_{xy} of lead track with respect to PV, in cm (with sign)'), Plot1D('dz', 'dz', 20, -20, 20, 'd_{z} of lead track with respect to PV, in cm (with sign)'), Plot1D('eta', 'eta', 20, -3, 3, 'eta'), @@ -797,6 +798,7 @@ Plot1D('idAntiEleDeadECal', 'idAntiEleDeadECal', 2, -0.5, 1.5, "tauID('againstElectronDeadECAL')"), Plot1D('idAntiMu', 'idAntiMu', 11, -0.5, 10.5, 'Anti-muon discriminator V3: : int 1 = Loose, 2 = Tight'), Plot1D('idDecayModeOldDMs', 'idDecayModeOldDMs', 2, -0.5, 1.5, "tauID('decayModeFinding')"), + Plot1D('idDecayModeNewDMs', 'idDecayModeNewDMs', 2, -0.5, 1.5, "tauID('decayModeFindingNewDMs')"), Plot1D('idDeepTau2017v2p1VSe', 'idDeepTau2017v2p1VSe', 11, -0.5, 10.5, 'byDeepTau2017v2p1VSe ID working points (deepTau2017v2p1): int 1 = VVVLoose, 2 = VVLoose, 3 = VLoose, 4 = Loose, 5 = Medium, 6 = Tight, 7 = VTight, 8 = VVTight'), Plot1D('idDeepTau2017v2p1VSjet', 'idDeepTau2017v2p1VSjet', 11, -0.5, 10.5, 'byDeepTau2017v2p1VSjet ID working points (deepTau2017v2p1): int 1 = VVVLoose, 2 = VVLoose, 3 = VLoose, 4 = Loose, 5 = Medium, 6 = Tight, 7 = VTight, 8 = VVTight'), Plot1D('idDeepTau2017v2p1VSmu', 'idDeepTau2017v2p1VSmu', 11, -0.5, 10.5, 'byDeepTau2017v2p1VSmu ID working points (deepTau2017v2p1): int 1 = VLoose, 2 = Loose, 3 = Medium, 4 = Tight'), @@ -819,8 +821,18 @@ Plot1D('rawDeepTau2018v2p5VSe', 'rawDeepTau2018v2p5VSe', 20, 0, 1, 'byDeepTau2018v2p5VSe raw output discriminator (deepTau2018v2p5)'), Plot1D('rawDeepTau2018v2p5VSjet', 'rawDeepTau2018v2p5VSjet', 20, 0, 1, 'byDeepTau2018v2p5VSjet raw output discriminator (deepTau2018v2p5)'), Plot1D('rawDeepTau2018v2p5VSmu', 'rawDeepTau2018v2p5VSmu', 20, 0, 1, 'byDeepTau2018v2p5VSmu raw output discriminator (deepTau2018v2p5)'), + Plot1D('rawPNetVSe', 'rawPNetVSe', 20, 0, 1, 'byPNetVSe raw output discriminator (PNet 2023)'), + Plot1D('rawPNetVSjet', 'rawPNetVSjet', 20, 0, 1, 'byPNetVSjet raw output discriminator (PNet 2023)'), + Plot1D('rawPNetVSmu', 'rawPNetVSmu', 20, 0, 1, 'byPNetVSmu raw output discriminator (PNet 2023)'), Plot1D('rawIso', 'rawIso', 20, 0, 200, 'combined isolation (deltaBeta corrections)'), Plot1D('rawIsodR03', 'rawIsodR03', 20, 0, 200, 'combined isolation (deltaBeta corrections, dR=0.3)'), + Plot1D('ptCorrPNet', 'ptCorrPNet', 20, 0, 2, 'pt correction (PNet 2023)'), + Plot1D('qConfPNet', 'qConfPNet', 20, -0.5, 0.5, 'signed charge confidence (PNet 2023)'), + Plot1D('probDM0PNet', 'probDM0PNet', 20, 0, 1, 'normalised probablity of decayMode 0, 1h+0pi0 (PNet 2023)'), + Plot1D('probDM1PNet', 'probDM1PNet', 20, 0, 1, 'normalised probablity of decayMode 1, 1h+1pi0 (PNet 2023)'), + Plot1D('probDM2PNet', 'probDM2PNet', 20, 0, 1, 'normalised probablity of decayMode 2, 1h+2pi0 (PNet 2023)'), + Plot1D('probDM10PNet', 'probDM10PNet', 20, 0, 1, 'normalised probablity of decayMode 10, 3h+0pi0 (PNet 2023)'), + Plot1D('probDM11PNet', 'probDM11PNet', 20, 0, 1, 'normalised probablity of decayMode 11, 3h+1pi0 (PNet 2023)'), ) ), TkMET = cms.PSet( diff --git a/PhysicsTools/NanoAOD/python/nano_cff.py b/PhysicsTools/NanoAOD/python/nano_cff.py index 1786b0016ef87..f2da46e9fb7b4 100644 --- a/PhysicsTools/NanoAOD/python/nano_cff.py +++ b/PhysicsTools/NanoAOD/python/nano_cff.py @@ -135,6 +135,51 @@ def nanoAOD_addBoostedTauIds(process, idsToRun=[]): return process +def nanoAOD_addPNetToTaus(process, addPNetInfo=False, runPNetCHSAK4=False): + if addPNetInfo: + originalTauName = process.finalTaus.src.value() + updatedTauName = originalTauName+'WithPNet' + jetCollection = "updatedJets" + process.load('RecoBTag.ONNXRuntime.pfParticleNetFromMiniAODAK4_cff') + pnetTagName = "pfParticleNetFromMiniAODAK4CHSCentralJetTag" + pnetDiscriminators = []; + for tag in getattr(process,pnetTagName+"s").flav_names.value(): + pnetDiscriminators.append(pnetTagName+"s:"+tag) + + # Define "hybridTau" producer + from PhysicsTools.PatAlgos.patTauHybridProducer_cfi import patTauHybridProducer + setattr(process, updatedTauName, patTauHybridProducer.clone( + src = originalTauName, + jetSource = jetCollection, + dRMax = 0.4, + jetPtMin = 15, + jetEtaMax = 2.5, + pnetLabel = pnetTagName+"s", + pnetScoreNames = pnetDiscriminators, + tauScoreMin = -1, + vsJetMin = 0.05, + checkTauScoreIsBest = False, + chargeAssignmentProbMin = 0.2, + addGenJetMatch = False, + genJetMatch = "" + )) + process.finalTaus.src = updatedTauName + + # run PNet for CHS AK4 jets if requested + if runPNetCHSAK4: + from PhysicsTools.NanoAOD.jetsAK4_CHS_cff import nanoAOD_addDeepInfoAK4CHS + process = nanoAOD_addDeepInfoAK4CHS(process, + addDeepBTag = False, + addDeepFlavour = False, + addParticleNet = True + ) + + #remember to adjust the selection and tables with added IDs + + process.tauTask.add(process.jetTask, getattr(process, updatedTauName)) + + return process + from PhysicsTools.SelectorUtils.tools.vid_id_tools import * def nanoAOD_activateVID(process): @@ -184,13 +229,26 @@ def nanoAOD_customizeCommon(process): ) nanoAOD_tau_switch = cms.PSet( - idsToAdd = cms.vstring() + idsToAdd = cms.vstring(), + runPNetAK4 = cms.bool(False), + addPNet = cms.bool(True) ) (run2_nanoAOD_106Xv2 | run3_nanoAOD_122).toModify( nanoAOD_tau_switch, idsToAdd = ["deepTau2018v2p5"] ).toModify( process, lambda p : nanoAOD_addTauIds(p, nanoAOD_tau_switch.idsToAdd.value()) ) + # Add PNet info to taus + # enable rerun of PNet for CHS jets for early run3 eras + # (it is rerun for run2 within jet tasks while is not needed for newer + # run3 eras as it is present in miniAOD) + (run3_nanoAOD_122 | run3_nanoAOD_124).toModify( + nanoAOD_tau_switch, runPNetAK4 = True + ) + nanoAOD_addPNetToTaus(process, + addPNetInfo = nanoAOD_tau_switch.addPNet.value(), + runPNetCHSAK4 = nanoAOD_tau_switch.runPNetAK4.value() + ) nanoAOD_boostedTau_switch = cms.PSet( idsToAdd = cms.vstring() ) diff --git a/PhysicsTools/NanoAOD/python/taus_cff.py b/PhysicsTools/NanoAOD/python/taus_cff.py index 8aa3f8077d7ba..f78f5b957511e 100644 --- a/PhysicsTools/NanoAOD/python/taus_cff.py +++ b/PhysicsTools/NanoAOD/python/taus_cff.py @@ -18,12 +18,12 @@ finalTaus = cms.EDFilter("PATTauRefSelector", src = cms.InputTag("slimmedTaus"), - cut = cms.string("pt > 18 && tauID('decayModeFindingNewDMs') && (tauID('byLooseCombinedIsolationDeltaBetaCorr3Hits') || (tauID('chargedIsoPtSumdR03')+max(0.,tauID('neutralIsoPtSumdR03')-0.072*tauID('puCorrPtSum'))<2.5) || tauID('byVVVLooseDeepTau2017v2p1VSjet') || tauID('byVVVLooseDeepTau2018v2p5VSjet'))") + cut = cms.string("pt > 18 && ((tauID('decayModeFindingNewDMs') > 0.5 && (tauID('byLooseCombinedIsolationDeltaBetaCorr3Hits') || (tauID('chargedIsoPtSumdR03')+max(0.,tauID('neutralIsoPtSumdR03')-0.072*tauID('puCorrPtSum'))<2.5) || tauID('byVVVLooseDeepTau2017v2p1VSjet') || tauID('byVVVLooseDeepTau2018v2p5VSjet'))) || (?isTauIDAvailable('byPNetVSjetraw')?tauID('byPNetVSjetraw'):-1) > {})".format(0.05)) ) run3_nanoAOD_124.toModify( finalTaus, - cut = cms.string("pt > 18 && tauID('decayModeFindingNewDMs') && (tauID('byLooseCombinedIsolationDeltaBetaCorr3Hits') || (tauID('chargedIsoPtSumdR03')+max(0.,tauID('neutralIsoPtSumdR03')-0.072*tauID('puCorrPtSum'))<2.5) || tauID('byVVVLooseDeepTau2017v2p1VSjet') || (tauID('byDeepTau2018v2p5VSjetraw') > {}))".format(WORKING_POINTS_v2p5["jet"]["VVVLoose"])) + cut = cms.string("pt > 18 && ((tauID('decayModeFindingNewDMs') > 0.5 && (tauID('byLooseCombinedIsolationDeltaBetaCorr3Hits') || (tauID('chargedIsoPtSumdR03')+max(0.,tauID('neutralIsoPtSumdR03')-0.072*tauID('puCorrPtSum'))<2.5) || tauID('byVVVLooseDeepTau2017v2p1VSjet') || (tauID('byDeepTau2018v2p5VSjetraw') > {}))) || (?isTauIDAvailable('byPNetVSjetraw')?tauID('byPNetVSjetraw'):-1) > {})".format(WORKING_POINTS_v2p5["jet"]["VVVLoose"],0.05)) ) ##################### Tables for final output and docs ########################## @@ -38,8 +38,10 @@ def _tauIdWPMask(pattern, choices, doc="", from_raw=False, wp_thrs=None): wp_definition = "test_bit(tauID('{}')-{}+1,0)".format(pattern, wp_thrs[wp_name]) var_definition.append(wp_definition) var_definition = " + ".join(var_definition) + var_definition = ("?isTauIDAvailable('%s')?(" % pattern) + var_definition + "):0" else: var_definition = " + ".join(["tauID('%s')" % (pattern % c) for c in choices]) + var_definition = ("?isTauIDAvailable('%s')?(" % (pattern % choices[0])) + var_definition + "):0" doc = doc + ": "+", ".join(["%d = %s" % (i,c) for (i,c) in enumerate(choices, start=1)]) return Var(var_definition, "uint8", doc=doc) @@ -60,32 +62,32 @@ def _tauIdWPMask(pattern, choices, doc="", from_raw=False, wp_thrs=None): svIdx2 = Var("?overlaps('vertices').size()>1?overlaps('vertices')[1].key():-1", "int16", doc="index of second matching secondary vertex"), nSVs = Var("?hasOverlaps('vertices')?overlaps('vertices').size():0", "uint8", doc="number of secondary vertices in the tau"), decayMode = Var("decayMode()", "uint8"), - idDecayModeOldDMs = Var("tauID('decayModeFinding')", bool), + idDecayModeOldDMs = Var("(?isTauIDAvailable('decayModeFinding')?tauID('decayModeFinding'):-1) > 0", bool), + idDecayModeNewDMs = Var("(?isTauIDAvailable('decayModeFindingNewDMs')?tauID('decayModeFindingNewDMs'):-1) > 0", bool), + leadTkPtOverTauPt = Var("?leadChargedHadrCand.isNonnull()?leadChargedHadrCand.pt/pt:1",float, doc="pt of the leading track divided by tau pt",precision=10), + leadTkDeltaEta = Var("?leadChargedHadrCand.isNonnull()?(leadChargedHadrCand.eta - eta):0",float, doc="eta of the leading track, minus tau eta",precision=8), + leadTkDeltaPhi = Var("?leadChargedHadrCand.isNonnull()?deltaPhi(leadChargedHadrCand.phi, phi):0",float, doc="phi of the leading track, minus tau phi",precision=8), - leadTkPtOverTauPt = Var("leadChargedHadrCand.pt/pt ",float, doc="pt of the leading track divided by tau pt",precision=10), - leadTkDeltaEta = Var("leadChargedHadrCand.eta - eta ",float, doc="eta of the leading track, minus tau eta",precision=8), - leadTkDeltaPhi = Var("deltaPhi(leadChargedHadrCand.phi, phi) ",float, doc="phi of the leading track, minus tau phi",precision=8), - - dxy = Var("leadChargedHadrCand().dxy()",float, doc="d_{xy} of lead track with respect to PV, in cm (with sign)",precision=10), - dz = Var("leadChargedHadrCand().dz()",float, doc="d_{z} of lead track with respect to PV, in cm (with sign)",precision=14), + dxy = Var("?leadChargedHadrCand.isNonnull()?leadChargedHadrCand().dxy():0",float, doc="d_{xy} of lead track with respect to PV, in cm (with sign)",precision=10), + dz = Var("?leadChargedHadrCand.isNonnull()?leadChargedHadrCand().dz():0",float, doc="d_{z} of lead track with respect to PV, in cm (with sign)",precision=14), # these are too many, we may have to suppress some - rawIso = Var( "tauID('byCombinedIsolationDeltaBetaCorrRaw3Hits')", float, doc = "combined isolation (deltaBeta corrections)", precision=10), - rawIsodR03 = Var( "(tauID('chargedIsoPtSumdR03')+max(0.,tauID('neutralIsoPtSumdR03')-0.072*tauID('puCorrPtSum')))", float, doc = "combined isolation (deltaBeta corrections, dR=0.3)", precision=10), - chargedIso = Var( "tauID('chargedIsoPtSum')", float, doc = "charged isolation", precision=10), - neutralIso = Var( "tauID('neutralIsoPtSum')", float, doc = "neutral (photon) isolation", precision=10), - puCorr = Var( "tauID('puCorrPtSum')", float, doc = "pileup correction", precision=10), - photonsOutsideSignalCone = Var( "tauID('photonPtSumOutsideSignalCone')", float, doc = "sum of photons outside signal cone", precision=10), + rawIso = Var("?isTauIDAvailable('byCombinedIsolationDeltaBetaCorrRaw3Hits')?tauID('byCombinedIsolationDeltaBetaCorrRaw3Hits'):-1", float, doc = "combined isolation (deltaBeta corrections)", precision=10), + rawIsodR03 = Var("?isTauIDAvailable('chargedIsoPtSumdR03')?(tauID('chargedIsoPtSumdR03')+max(0.,tauID('neutralIsoPtSumdR03')-0.072*tauID('puCorrPtSum'))):-1", float, doc = "combined isolation (deltaBeta corrections, dR=0.3)", precision=10), + chargedIso = Var("?isTauIDAvailable('chargedIsoPtSum')?tauID('chargedIsoPtSum'):-1", float, doc = "charged isolation", precision=10), + neutralIso = Var("?isTauIDAvailable('neutralIsoPtSum')?tauID('neutralIsoPtSum'):-1", float, doc = "neutral (photon) isolation", precision=10), + puCorr = Var("?isTauIDAvailable('puCorrPtSum')?tauID('puCorrPtSum'):-1", float, doc = "pileup correction", precision=10), + photonsOutsideSignalCone = Var("?isTauIDAvailable('photonPtSumOutsideSignalCone')?tauID('photonPtSumOutsideSignalCone'):-1", float, doc = "sum of photons outside signal cone", precision=10), idAntiMu = _tauIdWPMask("againstMuon%s3", choices=("Loose","Tight"), doc= "Anti-muon discriminator V3: "), - idAntiEleDeadECal = Var("tauID('againstElectronDeadECAL')", bool, doc = "Anti-electron dead-ECal discriminator"), + idAntiEleDeadECal = Var("(?isTauIDAvailable('againstElectronDeadECAL')?tauID('againstElectronDeadECAL'):-1) > 0", bool, doc = "Anti-electron dead-ECal discriminator"), ) _deepTauVars2017v2p1 = cms.PSet( - rawDeepTau2017v2p1VSe = Var("tauID('byDeepTau2017v2p1VSeraw')", float, doc="byDeepTau2017v2p1VSe raw output discriminator (deepTau2017v2p1)", precision=10), - rawDeepTau2017v2p1VSmu = Var("tauID('byDeepTau2017v2p1VSmuraw')", float, doc="byDeepTau2017v2p1VSmu raw output discriminator (deepTau2017v2p1)", precision=10), - rawDeepTau2017v2p1VSjet = Var("tauID('byDeepTau2017v2p1VSjetraw')", float, doc="byDeepTau2017v2p1VSjet raw output discriminator (deepTau2017v2p1)", precision=10), + rawDeepTau2017v2p1VSe = Var("?isTauIDAvailable('byDeepTau2017v2p1VSeraw')?tauID('byDeepTau2017v2p1VSeraw'):-1", float, doc="byDeepTau2017v2p1VSe raw output discriminator (deepTau2017v2p1)", precision=10), + rawDeepTau2017v2p1VSmu = Var("?isTauIDAvailable('byDeepTau2017v2p1VSmuraw')?tauID('byDeepTau2017v2p1VSmuraw'):-1", float, doc="byDeepTau2017v2p1VSmu raw output discriminator (deepTau2017v2p1)", precision=10), + rawDeepTau2017v2p1VSjet = Var("?isTauIDAvailable('byDeepTau2017v2p1VSjetraw')?tauID('byDeepTau2017v2p1VSjetraw'):-1", float, doc="byDeepTau2017v2p1VSjet raw output discriminator (deepTau2017v2p1)", precision=10), idDeepTau2017v2p1VSe = _tauIdWPMask("by%sDeepTau2017v2p1VSe", choices=("VVVLoose","VVLoose","VLoose","Loose","Medium","Tight","VTight","VVTight"), doc="byDeepTau2017v2p1VSe ID working points (deepTau2017v2p1)"), @@ -97,9 +99,9 @@ def _tauIdWPMask(pattern, choices, doc="", from_raw=False, wp_thrs=None): doc="byDeepTau2017v2p1VSjet ID working points (deepTau2017v2p1)"), ) _deepTauVars2018v2p5 = cms.PSet( - rawDeepTau2018v2p5VSe = Var("tauID('byDeepTau2018v2p5VSeraw')", float, doc="byDeepTau2018v2p5VSe raw output discriminator (deepTau2018v2p5)", precision=10), - rawDeepTau2018v2p5VSmu = Var("tauID('byDeepTau2018v2p5VSmuraw')", float, doc="byDeepTau2018v2p5VSmu raw output discriminator (deepTau2018v2p5)", precision=10), - rawDeepTau2018v2p5VSjet = Var("tauID('byDeepTau2018v2p5VSjetraw')", float, doc="byDeepTau2018v2p5VSjet raw output discriminator (deepTau2018v2p5)", precision=10), + rawDeepTau2018v2p5VSe = Var("?isTauIDAvailable('byDeepTau2018v2p5VSeraw')?tauID('byDeepTau2018v2p5VSeraw'):-1", float, doc="byDeepTau2018v2p5VSe raw output discriminator (deepTau2018v2p5)", precision=10), + rawDeepTau2018v2p5VSmu = Var("?isTauIDAvailable('byDeepTau2018v2p5VSmuraw')?tauID('byDeepTau2018v2p5VSmuraw'):-1", float, doc="byDeepTau2018v2p5VSmu raw output discriminator (deepTau2018v2p5)", precision=10), + rawDeepTau2018v2p5VSjet = Var("?isTauIDAvailable('byDeepTau2018v2p5VSjetraw')?tauID('byDeepTau2018v2p5VSjetraw'):-1", float, doc="byDeepTau2018v2p5VSjet raw output discriminator (deepTau2018v2p5)", precision=10), idDeepTau2018v2p5VSe = _tauIdWPMask("by%sDeepTau2018v2p5VSe", choices=("VVVLoose","VVLoose","VLoose","Loose","Medium","Tight","VTight","VVTight"), doc="byDeepTau2018v2p5VSe ID working points (deepTau2018v2p5)"), @@ -111,10 +113,25 @@ def _tauIdWPMask(pattern, choices, doc="", from_raw=False, wp_thrs=None): doc="byDeepTau2018v2p5VSjet ID working points (deepTau2018v2p5)"), ) +_pnet2023 = cms.PSet( + decayModePNet = Var("?isTauIDAvailable('byPNetDecayMode')?tauID('byPNetDecayMode'):-1", "int16",doc="decay mode of the highest tau score of ParticleNet"), + rawPNetVSe = Var("?isTauIDAvailable('byPNetVSeraw')?tauID('byPNetVSeraw'):-1", float, doc="raw output of ParticleNetVsE discriminator (PNet 2023)", precision=10), + rawPNetVSmu = Var("?isTauIDAvailable('byPNetVSmuraw')?tauID('byPNetVSmuraw'):-1", float, doc="raw output of ParticleNetVsMu discriminator (PNet 2023)", precision=10), + rawPNetVSjet = Var("?isTauIDAvailable('byPNetVSjetraw')?tauID('byPNetVSjetraw'):-1", float, doc="raw output of ParticleNetVsJet discriminator (PNet 2023)", precision=10), + ptCorrPNet = Var("?isTauIDAvailable('byPNetPtCorr')?tauID('byPNetPtCorr'):1", float, doc="pt correction (PNet 2023)", precision=10), + qConfPNet = Var("?isTauIDAvailable('byPNetQConf')?tauID('byPNetQConf'):0", float, doc="signed charge confidence (PNet 2023)", precision=10), + probDM0PNet = Var("?isTauIDAvailable('byPNetProb1h0pi0')?tauID('byPNetProb1h0pi0'):-1", float, doc="normalised probablity of decayMode 0, 1h+0pi0 (PNet 2023)", precision=10), + probDM1PNet = Var("?isTauIDAvailable('byPNetProb1h1pi0')?tauID('byPNetProb1h1pi0'):-1", float, doc="normalised probablity of decayMode 1, 1h+1pi0 (PNet 2023)", precision=10), + probDM2PNet = Var("?isTauIDAvailable('byPNetProb1h2pi0')?tauID('byPNetProb1h2pi0'):-1", float, doc="normalised probablity of decayMode 2, 1h+2pi0 (PNet 2023)", precision=10), + probDM10PNet = Var("?isTauIDAvailable('byPNetProb3h0pi0')?tauID('byPNetProb3h0pi0'):-1", float, doc="normalised probablity of decayMode 10, 3h+0pi0 (PNet 2023)", precision=10), + probDM11PNet = Var("?isTauIDAvailable('byPNetProb3h1pi0')?tauID('byPNetProb3h1pi0'):-1", float, doc="normalised probablity of decayMode 11, 3h+1pi0 (PNet 2023)", precision=10), +) + _variablesMiniV2 = cms.PSet( _tauVarsBase, _deepTauVars2017v2p1, - _deepTauVars2018v2p5 + _deepTauVars2018v2p5, + _pnet2023 ) tauTable.variables = _variablesMiniV2 diff --git a/PhysicsTools/PatAlgos/plugins/PATTauHybridProducer.cc b/PhysicsTools/PatAlgos/plugins/PATTauHybridProducer.cc new file mode 100644 index 0000000000000..e1d8d7c783c19 --- /dev/null +++ b/PhysicsTools/PatAlgos/plugins/PATTauHybridProducer.cc @@ -0,0 +1,450 @@ +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/PatCandidates/interface/Tau.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "FWCore/Utilities/interface/transform.h" +#include "DataFormats/JetReco/interface/GenJetCollection.h" +#include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" + +class PATTauHybridProducer : public edm::stream::EDProducer<> { +public: + explicit PATTauHybridProducer(const edm::ParameterSet&); + ~PATTauHybridProducer() override{}; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + void produce(edm::Event&, const edm::EventSetup&) override; + +private: + void fillTauFromJet(reco::PFTau& pfTau, const reco::JetBaseRef& jet); + //--- configuration parameters + edm::EDGetTokenT tausToken_; + edm::EDGetTokenT jetsToken_; + bool addGenJetMatch_; + edm::EDGetTokenT> genJetMatchToken_; + const float dR2Max_, jetPtMin_, jetEtaMax_; + const std::string pnetLabel_; + std::vector pnetTauScoreNames_; + std::vector pnetJetScoreNames_; + std::vector pnetLepScoreNames_; + std::string pnetPtCorrName_; + const float tauScoreMin_, vsJetMin_, chargeAssignmentProbMin_; + const bool checkTauScoreIsBest_; + const bool usePFLeptonsAsChargedHadrons_; + + const std::map tagToDM_; + enum class tauId_pnet_idx : size_t { dm = 0, vsjet, vse, vsmu, ptcorr, qconf, pdm0, pdm1, pdm2, pdm10, pdm11, last }; + enum class tauId_min_idx : size_t { hpsnew = 0, last }; +}; + +PATTauHybridProducer::PATTauHybridProducer(const edm::ParameterSet& cfg) + : tausToken_(consumes(cfg.getParameter("src"))), + jetsToken_(consumes(cfg.getParameter("jetSource"))), + addGenJetMatch_(cfg.getParameter("addGenJetMatch")), + dR2Max_(std::pow(cfg.getParameter("dRMax"), 2)), + jetPtMin_(cfg.getParameter("jetPtMin")), + jetEtaMax_(cfg.getParameter("jetEtaMax")), + pnetLabel_(cfg.getParameter("pnetLabel")), + pnetPtCorrName_(cfg.getParameter("pnetPtCorrName")), + tauScoreMin_(cfg.getParameter("tauScoreMin")), + vsJetMin_(cfg.getParameter("vsJetMin")), + chargeAssignmentProbMin_(cfg.getParameter("chargeAssignmentProbMin")), + checkTauScoreIsBest_(cfg.getParameter("checkTauScoreIsBest")), + usePFLeptonsAsChargedHadrons_(cfg.getParameter("usePFLeptonsAsChargedHadrons")), + tagToDM_({{"1h0p", 0}, {"1h1or2p", 1}, {"1h1p", 1}, {"1h2p", 2}, {"3h0p", 10}, {"3h1p", 11}}) { + // Read the different PNet score names + std::vector pnetScoreNames = cfg.getParameter>("pnetScoreNames"); + for (const auto& scoreName : pnetScoreNames) { + size_t labelLenght = scoreName.find(':') == std::string::npos ? 0 : scoreName.find(':') + 1; + std::string name = scoreName.substr(labelLenght); + if (name.find("prob") == std::string::npos) + continue; + if (name.find("tau") != std::string::npos) + pnetTauScoreNames_.push_back(name); + else if (name.find("ele") != std::string::npos || name.find("mu") != std::string::npos) + pnetLepScoreNames_.push_back(name); + else + pnetJetScoreNames_.push_back(name); + if (pnetPtCorrName_.find(':') != std::string::npos) + pnetPtCorrName_ = pnetPtCorrName_.substr(pnetPtCorrName_.find(':') + 1); + // GenJet matching + if (addGenJetMatch_) { + genJetMatchToken_ = + consumes>(cfg.getParameter("genJetMatch")); + } + } + + produces>(); + //FIXME: produce a separate collection for PNet-recovered taus? +} + +void PATTauHybridProducer::produce(edm::Event& evt, const edm::EventSetup& es) { + // Get the vector of taus + edm::Handle inputTaus; + evt.getByToken(tausToken_, inputTaus); + + auto outputTaus = std::make_unique>(); + outputTaus->reserve(inputTaus->size()); + + // Get the vector of jets + edm::Handle jets; + evt.getByToken(jetsToken_, jets); + + // Switch off gen-matching for real data + if (evt.isRealData()) { + addGenJetMatch_ = false; + } + edm::Handle> genJetMatch; + if (addGenJetMatch_) + evt.getByToken(genJetMatchToken_, genJetMatch); + + // Minimal HPS-like tauID list + std::vector tauIds_minimal((size_t)tauId_min_idx::last); + tauIds_minimal[(size_t)tauId_min_idx::hpsnew] = std::make_pair("decayModeFindingNewDMs", -1); + + // PNet tauID list + std::vector tauIds_pnet((size_t)tauId_pnet_idx::last); + tauIds_pnet[(size_t)tauId_pnet_idx::dm] = std::make_pair("byPNetDecayMode", reco::PFTau::kNull); + tauIds_pnet[(size_t)tauId_pnet_idx::vsjet] = std::make_pair("byPNetVSjetraw", -1); + tauIds_pnet[(size_t)tauId_pnet_idx::vse] = std::make_pair("byPNetVSeraw", -1); + tauIds_pnet[(size_t)tauId_pnet_idx::vsmu] = std::make_pair("byPNetVSmuraw", -1); + tauIds_pnet[(size_t)tauId_pnet_idx::ptcorr] = std::make_pair("byPNetPtCorr", 1); + tauIds_pnet[(size_t)tauId_pnet_idx::qconf] = std::make_pair("byPNetQConf", 0); + tauIds_pnet[(size_t)tauId_pnet_idx::pdm0] = std::make_pair("byPNetProb1h0pi0", -1); + tauIds_pnet[(size_t)tauId_pnet_idx::pdm1] = std::make_pair("byPNetProb1h1pi0", -1); + tauIds_pnet[(size_t)tauId_pnet_idx::pdm2] = std::make_pair("byPNetProb1h2pi0", -1); + tauIds_pnet[(size_t)tauId_pnet_idx::pdm10] = std::make_pair("byPNetProb3h0pi0", -1); + tauIds_pnet[(size_t)tauId_pnet_idx::pdm11] = std::make_pair("byPNetProb3h1pi0", -1); + + std::set matched_taus; + size_t jet_idx = 0; + for (const auto& jet : *jets) { + jet_idx++; + if (jet.pt() < jetPtMin_) + continue; + if (std::abs(jet.eta()) > jetEtaMax_) + continue; + size_t tau_idx = 0; + bool matched = false; + + // Analyse PNet scores + std::pair bestPnetTauScore("probtauundef", -1); + float sumOfPnetTauScores = 0; + std::vector tauPerDMScores(5); + float plusChargeProb = 0; + for (const auto& scoreName : pnetTauScoreNames_) { + float score = jet.bDiscriminator(pnetLabel_ + ":" + scoreName); + sumOfPnetTauScores += score; + if (scoreName.find("taup") != std::string::npos) + plusChargeProb += score; + if (score > bestPnetTauScore.second) { + bestPnetTauScore.second = score; + bestPnetTauScore.first = scoreName; + } + if (scoreName.find("1h0p") != std::string::npos) + tauPerDMScores[0] += score; + else if (scoreName.find("1h1") != + std::string:: + npos) //Note: final "p" in "1p" ommited to enble matching also with "1h1or2p" from early trainings + tauPerDMScores[1] += score; + else if (scoreName.find("1h2p") != std::string::npos) + tauPerDMScores[2] += score; + else if (scoreName.find("3h0p") != std::string::npos) + tauPerDMScores[3] += score; + else if (scoreName.find("3h1p") != std::string::npos) + tauPerDMScores[4] += score; + } + if (sumOfPnetTauScores > 0) + plusChargeProb /= sumOfPnetTauScores; + + float sumOfPnetEleScores = 0, sumOfPnetMuScores = 0; + bool isTauScoreBest = (sumOfPnetTauScores > 0); + for (const auto& scoreName : pnetLepScoreNames_) { + float score = jet.bDiscriminator(pnetLabel_ + ":" + scoreName); + if (scoreName.find("ele") != std::string::npos) + sumOfPnetEleScores += score; + else if (scoreName.find("mu") != std::string::npos) + sumOfPnetMuScores += score; + if (score > bestPnetTauScore.second) + isTauScoreBest = false; + } + if (checkTauScoreIsBest_ && isTauScoreBest) { //if needed iterate over jet scores + for (const auto& scoreName : pnetJetScoreNames_) + if (jet.bDiscriminator(pnetLabel_ + ":" + scoreName) > bestPnetTauScore.second) + isTauScoreBest = false; + } + + // PNet discriminants vs jets, electrons and muons + tauIds_pnet[(size_t)tauId_pnet_idx::vsjet].second = + sumOfPnetTauScores / + (1. - sumOfPnetEleScores - + sumOfPnetMuScores); //vsJet: tau scores by sum of tau and jet scores or equally by 1 - sum of lepton scores + tauIds_pnet[(size_t)tauId_pnet_idx::vse].second = + sumOfPnetTauScores / + (sumOfPnetTauScores + sumOfPnetEleScores); //vsEle: tau scores by sum of tau and ele scores + tauIds_pnet[(size_t)tauId_pnet_idx::vsmu].second = + sumOfPnetTauScores / (sumOfPnetTauScores + sumOfPnetMuScores); //vsMu: tau scores by sum of tau and mu scores + + // Decay mode and charge of the highest tau score + int bestCharge = 0; + size_t pos = + bestPnetTauScore.first.find("tau") + 3; //this is well defined by constraction as name is "probtauXXXX" + const char q = (pos < bestPnetTauScore.first.size()) ? bestPnetTauScore.first[pos] : 'u'; + if (q == 'm') { //minus + pos++; + bestCharge = -1; + } else if (q == 'p') { //plus + pos++; + bestCharge = 1; + } + auto pNetDM = tagToDM_.find(bestPnetTauScore.first.substr(pos)); + if (pNetDM != tagToDM_.end()) + tauIds_pnet[(size_t)tauId_pnet_idx::dm].second = pNetDM->second; + + // PNet Pt correction + float ptcorr = jet.bDiscriminator(pnetLabel_ + ":" + pnetPtCorrName_); + if (ptcorr > -1000.) // -1000. is default for not found discriminantor + tauIds_pnet[(size_t)tauId_pnet_idx::ptcorr].second = ptcorr; + + // PNet charge confidence + tauIds_pnet[(size_t)tauId_pnet_idx::qconf].second = (plusChargeProb - 0.5); + + // PNet per decay mode normalised score + tauIds_pnet[(size_t)tauId_pnet_idx::pdm0].second = tauPerDMScores[0] / sumOfPnetTauScores; + tauIds_pnet[(size_t)tauId_pnet_idx::pdm1].second = tauPerDMScores[1] / sumOfPnetTauScores; + tauIds_pnet[(size_t)tauId_pnet_idx::pdm2].second = tauPerDMScores[2] / sumOfPnetTauScores; + tauIds_pnet[(size_t)tauId_pnet_idx::pdm10].second = tauPerDMScores[3] / sumOfPnetTauScores; + tauIds_pnet[(size_t)tauId_pnet_idx::pdm11].second = tauPerDMScores[4] / sumOfPnetTauScores; + + // Search for matching tau + for (const auto& inputTau : *inputTaus) { + tau_idx++; + if (matched_taus.count(tau_idx - 1) > 0) + continue; + float dR2 = deltaR2(jet, inputTau); + // select 1st found match rather than best match (both should be equivalent for reasonable dRMax) + if (dR2 < dR2Max_) { + matched_taus.insert(tau_idx - 1); + pat::Tau outputTau(inputTau); + const size_t nTauIds = inputTau.tauIDs().size(); + std::vector tauIds(nTauIds + tauIds_pnet.size()); + for (size_t i = 0; i < nTauIds; ++i) + tauIds[i] = inputTau.tauIDs()[i]; + for (size_t i = 0; i < tauIds_pnet.size(); ++i) + tauIds[nTauIds + i] = tauIds_pnet[i]; + outputTau.setTauIDs(tauIds); + matched = true; + outputTaus->push_back(outputTau); + + break; + } + } // end of tau loop + if (matched) + continue; + + // Accept only jets passing minimal tau-like selection, i.e. with one of the tau score being globally the best and above some threshold, and with good quality of charge assignment + if ((checkTauScoreIsBest_ && !isTauScoreBest) || bestPnetTauScore.second < tauScoreMin_ || + tauIds_pnet[(size_t)tauId_pnet_idx::vsjet].second < vsJetMin_ || + std::abs(0.5 - plusChargeProb) < chargeAssignmentProbMin_) + continue; + + // Build taus from non-matched jets + // "Null" pftau with raw (uncorrected) jet kinematics + reco::PFTau pfTauFromJet(bestCharge, jet.correctedP4("Uncorrected")); + // Set PDGid + pfTauFromJet.setPdgId(bestCharge < 0 ? 15 : -15); + // and decay mode predicted by PNet + pfTauFromJet.setDecayMode( + static_cast(int(tauIds_pnet[(size_t)tauId_pnet_idx::dm].second))); + // Fill tau content using only jet consistunets within cone around leading + // charged particle + // FIXME: more sophisticated finding of tau constituents will be considered later + pfTauFromJet.setSignalConeSize( + std::clamp(3.6 / jet.correctedP4("Uncorrected").pt(), 0.08, 0.12)); // shrinking cone in function of jet-Pt + const edm::Ref jetRef(jets, jet_idx - 1); + fillTauFromJet(pfTauFromJet, reco::JetBaseRef(jetRef)); + + // PATTau + pat::Tau outputTauFromJet(pfTauFromJet); + // Add tauIDs + std::vector newtauIds(tauIds_minimal.size() + tauIds_pnet.size()); + for (size_t i = 0; i < tauIds_minimal.size(); ++i) + newtauIds[i] = tauIds_minimal[i]; + for (size_t i = 0; i < tauIds_pnet.size(); ++i) + newtauIds[tauIds_minimal.size() + i] = tauIds_pnet[i]; + outputTauFromJet.setTauIDs(newtauIds); + // Add genTauJet match + if (addGenJetMatch_) { + reco::GenJetRef genJetTau = (*genJetMatch)[jetRef]; + if (genJetTau.isNonnull() && genJetTau.isAvailable()) { + outputTauFromJet.setGenJet(genJetTau); + } + } + outputTaus->push_back(outputTauFromJet); + + } // end of jet loop + + // Taus non-matched to jets (usually at pt-threshold or/and eta boundaries) + if (matched_taus.size() < inputTaus->size()) { + for (size_t iTau = 0; iTau < inputTaus->size(); ++iTau) { + if (matched_taus.count(iTau) > 0) + continue; + const pat::Tau& inputTau = inputTaus->at(iTau); + pat::Tau outputTau(inputTau); + const size_t nTauIds = inputTau.tauIDs().size(); + std::vector tauIds(nTauIds + tauIds_pnet.size()); + for (size_t i = 0; i < nTauIds; ++i) + tauIds[i] = inputTau.tauIDs()[i]; + for (size_t i = 0; i < tauIds_pnet.size(); ++i) { + tauIds[nTauIds + i] = tauIds_pnet[i]; + tauIds[nTauIds + i].second = + (i != (size_t)tauId_pnet_idx::ptcorr ? (i != (size_t)tauId_pnet_idx::qconf ? -1 : 0) : 1); + } + outputTau.setTauIDs(tauIds); + outputTaus->push_back(outputTau); + } + } //non-matched taus + + evt.put(std::move(outputTaus)); +} + +void PATTauHybridProducer::fillTauFromJet(reco::PFTau& pfTau, const reco::JetBaseRef& jet) { + // Use tools as in PFTau builders to select tau decay products and isolation candidates + typedef std::vector CandPtrs; + + // Get the charged hadron candidates + CandPtrs pfChs, pfChsSig; + // Check if we want to include electrons and muons in "charged hadron" + // collection. This is the preferred behavior, as the PF lepton selections + // are very loose. + if (!usePFLeptonsAsChargedHadrons_) { + pfChs = reco::tau::pfCandidatesByPdgId(*jet, 211); + } else { + pfChs = reco::tau::pfChargedCands(*jet); + } + // take 1st charged candidate with charge as of tau (collection is pt-sorted) + if (pfTau.charge() == 0 || pfChs.size() == 1) { + pfTau.setleadChargedHadrCand(pfChs[0]); + pfTau.setleadCand(pfChs[0]); + pfChsSig.push_back(pfChs[0]); + pfChs.erase(pfChs.begin()); + } else { + for (CandPtrs::iterator it = pfChs.begin(); it != pfChs.end();) { + if ((*it)->charge() == pfTau.charge()) { + pfTau.setleadChargedHadrCand(*it); + pfTau.setleadCand(*it); + pfChsSig.push_back(*it); + it = pfChs.erase(it); + break; + } else { + ++it; + } + } + // In case of lack of candidate with charge same as of tau use leading charged candidate + if (pfTau.leadChargedHadrCand().isNull() && !pfChs.empty()) { + pfTau.setleadChargedHadrCand(pfChs[0]); + pfTau.setleadCand(pfChs[0]); + pfChsSig.push_back(pfChs[0]); + pfChs.erase(pfChs.begin()); + } + } + // if more than one charged decay product is expected add all inside signal + // cone around the leading track + const double dR2Max = std::pow(pfTau.signalConeSize(), 2); + if (pfTau.decayMode() >= reco::PFTau::kThreeProng0PiZero && pfTau.leadChargedHadrCand().isNonnull()) { + for (CandPtrs::iterator it = pfChs.begin(); it != pfChs.end();) { + if (deltaR2((*it)->p4(), pfTau.leadChargedHadrCand()->p4()) < dR2Max) { + pfChsSig.push_back(*it); + it = pfChs.erase(it); + } else { + ++it; + } + } + } + // Clean isolation candidates from low-pt and leptonic ones + for (CandPtrs::iterator it = pfChs.begin(); it != pfChs.end();) { + if ((*it)->pt() < 0.5 || std::abs((*it)->pdgId()) != 211) { + it = pfChs.erase(it); + } else { + ++it; + } + } + // Set charged candidates + pfTau.setsignalChargedHadrCands(pfChsSig); + pfTau.setisolationChargedHadrCands(pfChs); + + // Get the gamma candidates (pi0 decay products) + CandPtrs pfGammas, pfGammasSig; + pfGammas = reco::tau::pfCandidatesByPdgId(*jet, 22); + // In case of lack of leading charged candidate substiute it with leading gamma candidate + if (pfTau.leadChargedHadrCand().isNull() && !pfGammas.empty()) { + pfTau.setleadChargedHadrCand(pfGammas[0]); + pfTau.setleadCand(pfGammas[0]); + pfGammasSig.push_back(pfGammas[0]); + pfChs.erase(pfGammas.begin()); + } + // Clean gamma candidates from low-pt ones + for (CandPtrs::iterator it = pfGammas.begin(); it != pfGammas.end();) { + if ((*it)->pt() < 0.5) { + it = pfGammas.erase(it); + } else { + ++it; + } + } + // if decay mode with pi0s is expected look for signal gamma candidates + // within eta-phi strips around leading track + if (pfTau.decayMode() % 5 != 0 && pfTau.leadChargedHadrCand().isNonnull()) { + for (CandPtrs::iterator it = pfGammas.begin(); it != pfGammas.end();) { + if (std::abs((*it)->eta() - pfTau.leadChargedHadrCand()->eta()) < + std::clamp(0.2 * std::pow((*it)->pt(), -0.66), 0.05, 0.15) && + deltaPhi((*it)->phi(), pfTau.leadChargedHadrCand()->phi()) < + std::clamp(0.35 * std::pow((*it)->pt(), -0.71), 0.05, 0.3)) { + pfGammasSig.push_back(*it); + it = pfGammas.erase(it); + } else { + ++it; + } + } + } + // Set gamma candidates + pfTau.setsignalGammaCands(pfGammasSig); + pfTau.setisolationGammaCands(pfGammas); +} + +void PATTauHybridProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + // patTauHybridProducer + edm::ParameterSetDescription desc; + + desc.add("src", edm::InputTag("slimmedTaus")); + desc.add("jetSource", edm::InputTag("slimmedJetsUpadted")); + desc.add("dRMax", 0.4); + desc.add("jetPtMin", 20.0); + desc.add("jetEtaMax", 2.5); + desc.add("pnetLabel", "pfParticleNetAK4JetTags"); + desc.add>( + "pnetScoreNames", + {"probmu", "probele", "probtaup1h0p", "probtaup1h1p", "probtaup1h2p", "probtaup3h0p", "probtaup3h1p", + "probtaum1h0p", "probtaum1h1p", "probtaum1h2p", "probtaum3h0p", "probtaum3h1p", "probb", "probc", + "probuds", "probg", "ptcorr", "ptreshigh", "ptreslow", "ptnu"}); + desc.add("pnetPtCorrName", "ptcorr"); + desc.add("tauScoreMin", -1)->setComment("Minimal value of the best tau score to built recovery tau"); + desc.add("vsJetMin", -1)->setComment("Minimal value of PNet tau-vs-jet discriminant to built recovery tau"); + desc.add("checkTauScoreIsBest", false) + ->setComment("If true, recovery tau is built only if one of tau scores is the highest"); + desc.add("chargeAssignmentProbMin", 0.2) + ->setComment("Minimal value of charge assignment probability to built recovery tau (0,0.5)"); + desc.add("addGenJetMatch", true)->setComment("add MC genTauJet matching"); + desc.add("genJetMatch", edm::InputTag("tauGenJetMatch")); + desc.add("usePFLeptonsAsChargedHadrons", true) + ->setComment("If true, all charged particles are used as charged hadron candidates"); + + descriptions.addWithDefaultLabel(desc); +} + +#include "FWCore/Framework/interface/MakerMacros.h" + +DEFINE_FWK_MODULE(PATTauHybridProducer);