diff --git a/DataFormats/EgammaCandidates/interface/GsfElectron.h b/DataFormats/EgammaCandidates/interface/GsfElectron.h index 7e98ad1834645..e1b5c1a7c367a 100644 --- a/DataFormats/EgammaCandidates/interface/GsfElectron.h +++ b/DataFormats/EgammaCandidates/interface/GsfElectron.h @@ -787,6 +787,7 @@ namespace reco { // setters void setCorrectedEcalEnergyError(float newEnergyError); void setCorrectedEcalEnergy(float newEnergy); + void setCorrectedEcalEnergy(float newEnergy, bool rescaleDependentValues); void setTrackMomentumError(float trackMomentumError); void setP4(P4Kind kind, const LorentzVector &p4, float p4Error, bool setCandidate); using RecoCandidate::setP4; @@ -812,8 +813,8 @@ namespace reco { //bool isMomentumCorrected() const { return corrections_.isMomentumCorrected ; } float caloEnergy() const { return correctedEcalEnergy(); } bool isEnergyScaleCorrected() const { return isEcalEnergyCorrected(); } - void correctEcalEnergy(float newEnergy, float newEnergyError) { - setCorrectedEcalEnergy(newEnergy); + void correctEcalEnergy(float newEnergy, float newEnergyError, bool corrEovP = true) { + setCorrectedEcalEnergy(newEnergy, corrEovP); setEcalEnergyError(newEnergyError); } void correctMomentum(const LorentzVector &p4, float trackMomentumError, float p4Error) { diff --git a/DataFormats/EgammaCandidates/src/GsfElectron.cc b/DataFormats/EgammaCandidates/src/GsfElectron.cc index d2344432ebe77..3197b3cbc837c 100644 --- a/DataFormats/EgammaCandidates/src/GsfElectron.cc +++ b/DataFormats/EgammaCandidates/src/GsfElectron.cc @@ -171,14 +171,18 @@ void GsfElectron::setCorrectedEcalEnergyError(float energyError) { corrections_.correctedEcalEnergyError = energyError; } -void GsfElectron::setCorrectedEcalEnergy(float newEnergy) { +void GsfElectron::setCorrectedEcalEnergy(float newEnergy) { setCorrectedEcalEnergy(newEnergy, true); } + +void GsfElectron::setCorrectedEcalEnergy(float newEnergy, bool rescaleDependentValues) { math::XYZTLorentzVectorD momentum = p4(); momentum *= newEnergy / momentum.e(); setP4(momentum); - showerShape_.hcalDepth1OverEcal *= corrections_.correctedEcalEnergy / newEnergy; - showerShape_.hcalDepth2OverEcal *= corrections_.correctedEcalEnergy / newEnergy; - trackClusterMatching_.eSuperClusterOverP *= newEnergy / corrections_.correctedEcalEnergy; - corrections_.correctedEcalEnergyError *= newEnergy / corrections_.correctedEcalEnergy; + if (corrections_.correctedEcalEnergy > 0. && rescaleDependentValues) { + showerShape_.hcalDepth1OverEcal *= corrections_.correctedEcalEnergy / newEnergy; + showerShape_.hcalDepth2OverEcal *= corrections_.correctedEcalEnergy / newEnergy; + trackClusterMatching_.eSuperClusterOverP *= newEnergy / corrections_.correctedEcalEnergy; + corrections_.correctedEcalEnergyError *= newEnergy / corrections_.correctedEcalEnergy; + } corrections_.correctedEcalEnergy = newEnergy; corrections_.isEcalEnergyCorrected = true; } diff --git a/PhysicsTools/NanoAOD/python/lowPtElectrons_cff.py b/PhysicsTools/NanoAOD/python/lowPtElectrons_cff.py new file mode 100644 index 0000000000000..571bc72fcfd0a --- /dev/null +++ b/PhysicsTools/NanoAOD/python/lowPtElectrons_cff.py @@ -0,0 +1,203 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.nano_eras_cff import * +from PhysicsTools.NanoAOD.common_cff import * + +################################################################################ +# Modules +################################################################################ + +from RecoEgamma.EgammaTools.lowPtElectronModifier_cfi import lowPtElectronModifier +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectrons_cfi import lowPtRegressionModifier +modifiedLowPtElectrons = cms.EDProducer( + "ModifiedElectronProducer", + src = cms.InputTag("slimmedLowPtElectrons"), + modifierConfig = cms.PSet( + modifications = cms.VPSet(lowPtElectronModifier,lowPtRegressionModifier) + ) +) + +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronID_cfi import lowPtGsfElectronID +lowPtPATElectronID = lowPtGsfElectronID.clone( + usePAT = True, + electrons = "modifiedLowPtElectrons", + unbiased = "", + ModelWeights = [ + 'RecoEgamma/ElectronIdentification/data/LowPtElectrons/LowPtElectrons_ID_2020Nov28.root', + ], +) + +isoForLowPtEle = cms.EDProducer( + "EleIsoValueMapProducer", + src = cms.InputTag("modifiedLowPtElectrons"), + relative = cms.bool(False), + rho_MiniIso = cms.InputTag("fixedGridRhoFastjetAll"), + rho_PFIso = cms.InputTag("fixedGridRhoFastjetAll"), + EAFile_MiniIso = cms.FileInPath("RecoEgamma/ElectronIdentification/data/Fall17/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_94X.txt"), + EAFile_PFIso = cms.FileInPath("RecoEgamma/ElectronIdentification/data/Fall17/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_94X.txt"), +) + +updatedLowPtElectronsWithUserData = cms.EDProducer( + "PATElectronUserDataEmbedder", + src = cms.InputTag("modifiedLowPtElectrons"), + userFloats = cms.PSet( + ID = cms.InputTag("lowPtPATElectronID"), + miniIsoChg = cms.InputTag("isoForLowPtEle:miniIsoChg"), + miniIsoAll = cms.InputTag("isoForLowPtEle:miniIsoAll"), + ), + userIntFromBools = cms.PSet(), + userInts = cms.PSet(), + userCands = cms.PSet(), +) + +finalLowPtElectrons = cms.EDFilter( + "PATElectronRefSelector", + src = cms.InputTag("updatedLowPtElectronsWithUserData"), + cut = cms.string("pt > 1. && userFloat('ID') > -0.25"), +) + +################################################################################ +# electronTable +################################################################################ + +lowPtElectronTable = cms.EDProducer( + "SimpleCandidateFlatTableProducer", + src = cms.InputTag("finalLowPtElectrons"), + cut = cms.string(""), + name= cms.string("LowPtElectron"), + doc = cms.string("slimmedLowPtElectrons after basic selection (" + finalLowPtElectrons.cut.value()+")"), + singleton = cms.bool(False), # the number of entries is variable + extension = cms.bool(False), # this is the main table for the electrons + variables = cms.PSet( + # Basic variables + CandVars, + # BDT scores and WPs + embeddedID = Var("electronID('ID')",float,doc="ID, BDT (raw) score"), + ID = Var("userFloat('ID')",float,doc="New ID, BDT (raw) score"), + unbiased = Var("electronID('unbiased')",float,doc="ElectronSeed, pT- and dxy- agnostic BDT (raw) score"), + ptbiased = Var("electronID('ptbiased')",float,doc="ElectronSeed, pT- and dxy- dependent BDT (raw) score"), + # Isolation + miniPFRelIso_chg = Var("userFloat('miniIsoChg')/pt",float, + doc="mini PF relative isolation, charged component"), + miniPFRelIso_all = Var("userFloat('miniIsoAll')/pt",float, + doc="mini PF relative isolation, total (with scaled rho*EA PU corrections)"), + # Conversions + convVeto = Var("passConversionVeto()",bool,doc="pass conversion veto"), + convWP = Var("userInt('convOpen')*1 + userInt('convLoose')*2 + userInt('convTight')*4", + int,doc="conversion flag bit map: 1=Veto, 2=Loose, 3=Tight"), + convVtxRadius = Var("userFloat('convVtxRadius')",float,doc="conversion vertex radius (cm)",precision=7), + # Tracking + lostHits = Var("gsfTrack.hitPattern.numberOfLostHits('MISSING_INNER_HITS')","uint8",doc="number of missing inner hits"), + # Cluster-related + energyErr = Var("p4Error('P4_COMBINATION')",float,doc="energy error of the cluster-track combination",precision=6), + deltaEtaSC = Var("superCluster().eta()-eta()",float,doc="delta eta (SC,ele) with sign",precision=10), + r9 = Var("full5x5_r9()",float,doc="R9 of the SC, calculated with full 5x5 region",precision=10), + sieie = Var("full5x5_sigmaIetaIeta()",float,doc="sigma_IetaIeta of the SC, calculated with full 5x5 region",precision=10), + eInvMinusPInv = Var("(1-eSuperClusterOverP())/ecalEnergy()",float,doc="1/E_SC - 1/p_trk",precision=10), + scEtOverPt = Var("(superCluster().energy()/(pt*cosh(superCluster().eta())))-1",float,doc="(SC energy)/pt-1",precision=8), + hoe = Var("hadronicOverEm()",float,doc="H over E",precision=8), + # Displacement + dxy = Var("dB('PV2D')",float,doc="dxy (with sign) wrt first PV, in cm",precision=10), + dxyErr = Var("edB('PV2D')",float,doc="dxy uncertainty, in cm",precision=6), + dz = Var("dB('PVDZ')",float,doc="dz (with sign) wrt first PV, in cm",precision=10), + dzErr = Var("abs(edB('PVDZ'))",float,doc="dz uncertainty, in cm",precision=6), + ip3d = Var("abs(dB('PV3D'))",float,doc="3D impact parameter wrt first PV, in cm",precision=10), + sip3d = Var("abs(dB('PV3D')/edB('PV3D'))",float,doc="3D impact parameter significance wrt first PV, in cm",precision=10), + # Cross-referencing + #jetIdx + #photonIdx + ), +) + +################################################################################ +# electronTable (MC) +################################################################################ + +from PhysicsTools.NanoAOD.particlelevel_cff import particleLevel +particleLevelForMatchingLowPt = particleLevel.clone( + lepMinPt = cms.double(1.), + phoMinPt = cms.double(1), +) + +tautaggerForMatchingLowPt = cms.EDProducer( + "GenJetTauTaggerProducer", + src = cms.InputTag('particleLevelForMatchingLowPt:leptons') +) + +matchingLowPtElecPhoton = cms.EDProducer( + "GenJetGenPartMerger", + srcJet =cms.InputTag("particleLevelForMatchingLowPt:leptons"), + srcPart=cms.InputTag("particleLevelForMatchingLowPt:photons"), + hasTauAnc=cms.InputTag("tautaggerForMatchingLowPt"), +) + +lowPtElectronsMCMatchForTableAlt = cms.EDProducer( + "GenJetMatcherDRPtByDR", # cut on deltaR, deltaPt/Pt; pick best by deltaR + src = lowPtElectronTable.src, # final reco collection + matched = cms.InputTag("matchingLowPtElecPhoton:merged"), # final mc-truth particle collection + mcPdgId = cms.vint32(11,22), # one or more PDG ID (11 = el, 22 = pho); absolute values (see below) + checkCharge = cms.bool(False), # True = require RECO and MC objects to have the same charge + mcStatus = cms.vint32(), + maxDeltaR = cms.double(0.3), # Minimum deltaR for the match + maxDPtRel = cms.double(0.5), # Minimum deltaPt/Pt for the match + resolveAmbiguities = cms.bool(True), # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality = cms.bool(True), # False = just match input in order; True = pick lowest deltaR pair first +) + +lowPtElectronsMCMatchForTable = cms.EDProducer( + "MCMatcher", # cut on deltaR, deltaPt/Pt; pick best by deltaR + src = lowPtElectronTable.src, # final reco collection + matched = cms.InputTag("finalGenParticles"), # final mc-truth particle collection + mcPdgId = cms.vint32(11), # one or more PDG ID (11 = ele); absolute values (see below) + checkCharge = cms.bool(False), # True = require RECO and MC objects to have the same charge + mcStatus = cms.vint32(1), # PYTHIA status code (1 = stable, 2 = shower, 3 = hard scattering) + maxDeltaR = cms.double(0.3), # Minimum deltaR for the match + maxDPtRel = cms.double(0.5), # Minimum deltaPt/Pt for the match + resolveAmbiguities = cms.bool(True), # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality = cms.bool(True), # False = just match input in order; True = pick lowest deltaR pair first +) + +from PhysicsTools.NanoAOD.electrons_cff import electronMCTable +lowPtElectronMCTable = cms.EDProducer( + "CandMCMatchTableProducer", + src = lowPtElectronTable.src, + mcMapDressedLep = cms.InputTag("lowPtElectronsMCMatchForTableAlt"), + mcMap = cms.InputTag("lowPtElectronsMCMatchForTable"), + mapTauAnc = cms.InputTag("matchingLowPtElecPhoton:hasTauAnc"), + objName = lowPtElectronTable.name, + objType = electronMCTable.objType, + branchName = cms.string("genPart"), + docString = cms.string("MC matching to status==1 electrons or photons"), + genparticles = cms.InputTag("finalGenParticles"), +) + +################################################################################ +# Sequences +################################################################################ + +lowPtElectronSequence = cms.Sequence(modifiedLowPtElectrons + +lowPtPATElectronID + +isoForLowPtEle + +updatedLowPtElectronsWithUserData + +finalLowPtElectrons) +lowPtElectronTables = cms.Sequence(lowPtElectronTable) +lowPtElectronMC = cms.Sequence( + particleLevelForMatchingLowPt + +tautaggerForMatchingLowPt + +matchingLowPtElecPhoton + +lowPtElectronsMCMatchForTable + +lowPtElectronsMCMatchForTableAlt + +lowPtElectronMCTable) + +################################################################################ +# Modifiers +################################################################################ + +_modifiers = ( run2_miniAOD_80XLegacy | + run2_nanoAOD_94XMiniAODv1 | + run2_nanoAOD_94XMiniAODv2 | + run2_nanoAOD_94X2016 | + run2_nanoAOD_102Xv1 | + run2_nanoAOD_106Xv1 ) +(_modifiers).toReplaceWith(lowPtElectronSequence,cms.Sequence()) +(_modifiers).toReplaceWith(lowPtElectronTables,cms.Sequence()) +(_modifiers).toReplaceWith(lowPtElectronMC,cms.Sequence()) diff --git a/PhysicsTools/NanoAOD/python/nanoDQM_cff.py b/PhysicsTools/NanoAOD/python/nanoDQM_cff.py index 1772e9d76f736..1d74b21d089e2 100644 --- a/PhysicsTools/NanoAOD/python/nanoDQM_cff.py +++ b/PhysicsTools/NanoAOD/python/nanoDQM_cff.py @@ -108,6 +108,7 @@ ## MC nanoDQMMC = nanoDQM.clone() nanoDQMMC.vplots.Electron.sels.Prompt = cms.string("genPartFlav == 1") +nanoDQMMC.vplots.LowPtElectron.sels.Prompt = cms.string("genPartFlav == 1") nanoDQMMC.vplots.Muon.sels.Prompt = cms.string("genPartFlav == 1") nanoDQMMC.vplots.Photon.sels.Prompt = cms.string("genPartFlav == 1") nanoDQMMC.vplots.Tau.sels.Prompt = cms.string("genPartFlav == 5") diff --git a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py index ee2867025e2d6..19d6f305b16f2 100644 --- a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py +++ b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py @@ -112,6 +112,53 @@ Plot1D('dEsigmaDown', 'dEsigmaDown', 100, -0.1, 0.1, '#Delta E sigmaDown'), ) ), + + LowPtElectron = cms.PSet( + sels = cms.PSet( + Good = cms.string('pt > 1. && ID > 5.') + ), + plots = cms.VPSet( + # + Count1D('_size', 8, -0.5, 7.5, 'slimmedLowPtElectrons after basic selection'), + # CandVars + Plot1D('charge', 'charge', 3, -1.5, 1.5, 'electric charge'), + Plot1D('eta', 'eta', 20, -3., 3., 'eta'), + NoPlot('mass'), + Plot1D('pdgId', 'pdgId', 101, -50.5, 50.5, 'PDG code assigned by the event reconstruction (not by MC truth)'), + Plot1D('phi', 'phi', 20, -3.14159, 3.14159, 'phi'), + Plot1D('pt', 'pt', 40, 0., 20., 'pt (corrected)'), + # BDT scores and WPs + Plot1D('embeddedID', 'embeddedID', 40, -10., 10., 'Embedded ID, BDT (raw) score'), + Plot1D('ID', 'ID', 40, -10., 10., 'ID, BDT (raw) score'), + Plot1D('unbiased', 'unbiased', 40, -10., 10., 'ElectronSeed, pT- and dxy- agnostic BDT (raw) score'), + Plot1D('ptbiased', 'ptbiased', 40, -10., 10., 'ElectronSeed, pT- and dxy- dependent BDT (raw) score'), + # Isolation + Plot1D('miniPFRelIso_chg', 'miniPFRelIso_chg', 20, 0, 2, 'mini PF relative isolation, charged component'), + Plot1D('miniPFRelIso_all', 'miniPFRelIso_all', 20, 0, 2, 'mini PF relative isolation, total (with scaled rho*EA PU corrections)'), + # Conversions + Plot1D('convVeto', 'convVeto', 2, -0.5, 1.5, 'pass conversion veto'), + Plot1D('convWP', 'convWP', 8, -0.5, 7.5, 'conversion flag bit map: 1=Veto, 2=Loose, 3=Tight'), + Plot1D('convVtxRadius', 'convVtxRadius', 40, 0., 20.0, 'conversion vertex radius (cm)'), + # Tracking + Plot1D('lostHits', 'lostHits', 4, -0.5, 3.5, 'number of missing inner hits'), + # Cluster-related + Plot1D('energyErr', 'energyErr', 40, 0., 20., 'energy error of the cluster from regression'), + Plot1D('deltaEtaSC', 'deltaEtaSC', 20, -0.2, 0.2, 'delta eta (SC,ele) with sign'), + Plot1D('r9', 'r9', 20, 0, 1.1, 'R9 of the supercluster, calculated with full 5x5 region'), + Plot1D('sieie', 'sieie', 20, 0, 0.05, 'sigma_IetaIeta of the supercluster, calculated with full 5x5 region'), + Plot1D('eInvMinusPInv', 'eInvMinusPInv', 20, -0.1, 0.1, '1/E_SC - 1/p_trk'), + Plot1D('scEtOverPt', 'scEtOverPt', 20, -0.5, 0.5, '(supercluster transverse energy)/pt - 1'), + Plot1D('hoe', 'hoe', 20, 0, 0.6, 'H over E'), + # Displacement + Plot1D('dxy', 'dxy', 20, -0.1, 0.1, 'dxy (with sign) wrt first PV, in cm'), + Plot1D('dz', 'dz', 20, -0.3, 0.3, 'dz (with sign) wrt first PV, in cm'), + Plot1D('dxyErr', 'dxyErr', 20, 0., 0.2, 'dxy uncertainty, in cm'), + Plot1D('dzErr', 'dzErr', 20, 0., 0.2, 'dz uncertainty, in cm'), + Plot1D('ip3d', 'ip3d', 20, 0., 0.2, '3D impact parameter wrt first PV, in cm'), + Plot1D('sip3d', 'sip3d', 20, 0., 20., '3D impact parameter significance wrt first PV, in cm'), + ), + ), + FatJet = cms.PSet( sels = cms.PSet(), plots = cms.VPSet( diff --git a/PhysicsTools/NanoAOD/python/nano_cff.py b/PhysicsTools/NanoAOD/python/nano_cff.py index 54d614355fa4d..b26bff15d8ffc 100644 --- a/PhysicsTools/NanoAOD/python/nano_cff.py +++ b/PhysicsTools/NanoAOD/python/nano_cff.py @@ -7,6 +7,7 @@ from PhysicsTools.NanoAOD.taus_cff import * from PhysicsTools.NanoAOD.boostedTaus_cff import * from PhysicsTools.NanoAOD.electrons_cff import * +from PhysicsTools.NanoAOD.lowPtElectrons_cff import * from PhysicsTools.NanoAOD.photons_cff import * from PhysicsTools.NanoAOD.globals_cff import * from PhysicsTools.NanoAOD.extraflags_cff import * @@ -106,10 +107,10 @@ (run2_miniAOD_80XLegacy | run2_nanoAOD_94X2016 | run2_nanoAOD_94XMiniAODv1 | run2_nanoAOD_94XMiniAODv2 | run2_nanoAOD_102Xv1).toModify(l1bits, storeUnprefireableBit=False) nanoSequenceCommon = cms.Sequence( - nanoMetadata + jetSequence + muonSequence + tauSequence + boostedTauSequence + electronSequence+photonSequence+vertexSequence+ + nanoMetadata + jetSequence + muonSequence + tauSequence + boostedTauSequence + electronSequence + lowPtElectronSequence + photonSequence+vertexSequence+ isoTrackSequence + jetLepSequence + # must be after all the leptons linkedObjects + - jetTables + muonTables + tauTables + boostedTauTables + electronTables + photonTables + globalTables +vertexTables+ metTables+simpleCleanerTable + isoTrackTables + jetTables + muonTables + tauTables + boostedTauTables + electronTables + lowPtElectronTables + photonTables + globalTables +vertexTables+ metTables+simpleCleanerTable + isoTrackTables ) #remove boosted tau from previous eras (run2_miniAOD_80XLegacy | run2_nanoAOD_92X | run2_nanoAOD_94XMiniAODv1 | run2_nanoAOD_94X2016 | run2_nanoAOD_94XMiniAODv2 | run2_nanoAOD_102Xv1 | run2_nanoAOD_106Xv1).toReplaceWith(nanoSequenceCommon, nanoSequenceCommon.copyAndExclude([boostedTauSequence, boostedTauTables])) @@ -119,7 +120,7 @@ nanoSequence = cms.Sequence(nanoSequenceCommon + nanoSequenceOnlyData + nanoSequenceOnlyFullSim) -nanoSequenceFS = cms.Sequence(genParticleSequence + genVertexTables + particleLevelSequence + nanoSequenceCommon + jetMC + muonMC + electronMC + photonMC + tauMC + boostedTauMC + metMC + ttbarCatMCProducers + globalTablesMC + btagWeightTable + genWeightsTable + genVertexTable + genParticleTables + particleLevelTables + lheInfoTable + ttbarCategoryTable ) +nanoSequenceFS = cms.Sequence(genParticleSequence + genVertexTables + particleLevelSequence + nanoSequenceCommon + jetMC + muonMC + electronMC + lowPtElectronMC + photonMC + tauMC + boostedTauMC + metMC + ttbarCatMCProducers + globalTablesMC + btagWeightTable + genWeightsTable + genVertexTable + genParticleTables + particleLevelTables + lheInfoTable + ttbarCategoryTable ) (run2_nanoAOD_92X | run2_miniAOD_80XLegacy | run2_nanoAOD_94X2016 | run2_nanoAOD_94X2016 | \ run2_nanoAOD_94XMiniAODv1 | run2_nanoAOD_94XMiniAODv2 | \ diff --git a/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py b/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py index 684512a4ce0fd..df6d1c4180143 100644 --- a/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py +++ b/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py @@ -1,5 +1,7 @@ import FWCore.ParameterSet.Config as cms +from RecoEgamma.EgammaTools.lowPtElectronModifier_cfi import lowPtElectronModifier + slimmedLowPtElectrons = cms.EDProducer("PATElectronSlimmer", src = cms.InputTag("selectedPatLowPtElectrons"), dropSuperCluster = cms.string("0"), # you can put a cut to slim selectively, e.g. pt < 10 @@ -39,6 +41,7 @@ modifierName = cms.string('EGExtraInfoModifierFromPackedCandPtrValueMaps'), photon_config = cms.PSet() ), + lowPtElectronModifier, ) ) ) diff --git a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc index bba7853345515..ffa52a3468caa 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc +++ b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc @@ -37,9 +37,11 @@ class LowPtGsfElectronIDProducer final : public edm::global::EDProducer<> { double eval( const std::string& name, const edm::Ptr&, double rho, float unbiased, float field_z) const; - const edm::EDGetTokenT > electrons_; + const bool usePAT_; + edm::EDGetTokenT electrons_; + edm::EDGetTokenT patElectrons_; const edm::EDGetTokenT rho_; - const edm::EDGetTokenT > unbiased_; + edm::EDGetTokenT > unbiased_; const std::vector names_; const bool passThrough_; const double minPtThreshold_; @@ -52,15 +54,23 @@ class LowPtGsfElectronIDProducer final : public edm::global::EDProducer<> { //////////////////////////////////////////////////////////////////////////////// // LowPtGsfElectronIDProducer::LowPtGsfElectronIDProducer(const edm::ParameterSet& conf) - : electrons_(consumes >(conf.getParameter("electrons"))), + : usePAT_(conf.getParameter("usePAT")), + electrons_(), + patElectrons_(), rho_(consumes(conf.getParameter("rho"))), - unbiased_(consumes >(conf.getParameter("unbiased"))), + unbiased_(), names_(conf.getParameter >("ModelNames")), passThrough_(conf.getParameter("PassThrough")), minPtThreshold_(conf.getParameter("MinPtThreshold")), maxPtThreshold_(conf.getParameter("MaxPtThreshold")), thresholds_(conf.getParameter >("ModelThresholds")), version_(conf.getParameter("Version")) { + if (usePAT_) { + patElectrons_ = consumes(conf.getParameter("electrons")); + } else { + electrons_ = consumes(conf.getParameter("electrons")); + unbiased_ = consumes >(conf.getParameter("unbiased")); + } for (auto& weights : conf.getParameter >("ModelWeights")) { models_.push_back(createGBRForest(edm::FileInPath(weights))); } @@ -97,39 +107,52 @@ void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const throw cms::Exception("InvalidHandle", os.str()); } - // Retrieve GsfElectrons from Event - edm::Handle > electrons; - event.getByToken(electrons_, electrons); - if (!electrons.isValid()) { - std::ostringstream os; - os << "Problem accessing low-pT electrons collection" << std::endl; - throw cms::Exception("InvalidHandle", os.str()); + // Retrieve pat::Electrons or reco::GsfElectrons from Event + edm::Handle patElectrons; + edm::Handle electrons; + if (usePAT_) { + event.getByToken(patElectrons_, patElectrons); + } else { + event.getByToken(electrons_, electrons); } // ElectronSeed unbiased BDT edm::Handle > unbiasedH; - event.getByToken(unbiased_, unbiasedH); + if (!unbiased_.isUninitialized()) { + event.getByToken(unbiased_, unbiasedH); + } // Iterate through Electrons, evaluate BDT, and store result std::vector > output; + unsigned int nElectrons = usePAT_ ? patElectrons->size() : electrons->size(); for (unsigned int iname = 0; iname < names_.size(); ++iname) { - output.emplace_back(electrons->size(), -999.); + output.emplace_back(nElectrons, -999.); } - for (unsigned int iele = 0; iele < electrons->size(); iele++) { - edm::Ptr ele(electrons, iele); - if (ele->core().isNull()) { - continue; - } - const auto& gsf = ele->core()->gsfTrack(); // reco::GsfTrackRef - if (gsf.isNull()) { - continue; + if (usePAT_) { + for (unsigned int iele = 0; iele < nElectrons; iele++) { + edm::Ptr ele(patElectrons, iele); + if (!ele->isElectronIDAvailable("unbiased")) { + continue; + } + for (unsigned int iname = 0; iname < names_.size(); ++iname) { + output[iname][iele] = eval(names_[iname], ele, *rho, ele->electronID("unbiased"), zfield.z()); + } } - float unbiased = (*unbiasedH)[gsf]; - - //if ( !passThrough_ && ( ele->pt() < minPtThreshold_ ) ) { continue; } - for (unsigned int iname = 0; iname < names_.size(); ++iname) { - output[iname][iele] = eval(names_[iname], ele, *rho, unbiased, zfield.z()); + } else { + for (unsigned int iele = 0; iele < nElectrons; iele++) { + edm::Ptr ele(electrons, iele); + if (ele->core().isNull()) { + continue; + } + const auto& gsf = ele->core()->gsfTrack(); // reco::GsfTrackRef + if (gsf.isNull()) { + continue; + } + float unbiased = (*unbiasedH)[gsf]; + for (unsigned int iname = 0; iname < names_.size(); ++iname) { + output[iname][iele] = eval(names_[iname], ele, *rho, unbiased, zfield.z()); + } } } @@ -137,7 +160,11 @@ void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const for (unsigned int iname = 0; iname < names_.size(); ++iname) { auto ptr = std::make_unique >(edm::ValueMap()); edm::ValueMap::Filler filler(*ptr); - filler.insert(electrons, output[iname].begin(), output[iname].end()); + if (usePAT_) { + filler.insert(patElectrons, output[iname].begin(), output[iname].end()); + } else { + filler.insert(electrons, output[iname].begin(), output[iname].end()); + } filler.fill(); event.put(std::move(ptr), names_[iname]); } @@ -167,8 +194,9 @@ double LowPtGsfElectronIDProducer::eval( // void LowPtGsfElectronIDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; + desc.add("usePAT", false); desc.add("electrons", edm::InputTag("lowPtGsfElectrons")); - desc.add("unbiased", edm::InputTag("lowPtGsfElectronSeedValueMaps:unbiased")); + desc.addOptional("unbiased", edm::InputTag("lowPtGsfElectronSeedValueMaps:unbiased")); desc.add("rho", edm::InputTag("fixedGridRhoFastjetAll")); desc.add >("ModelNames", {""}); desc.add >( diff --git a/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectrons_cfi.py b/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectrons_cfi.py index 272a416aa6cf2..ed420bdc39921 100644 --- a/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectrons_cfi.py +++ b/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectrons_cfi.py @@ -2,7 +2,7 @@ from RecoEgamma.EgammaTools.regressionModifier_cfi import regressionModifier106XUL -_lowPtRegressionModifier = regressionModifier106XUL.clone( +lowPtRegressionModifier = regressionModifier106XUL.clone( modifierName = 'EGRegressionModifierV3', rhoTag = 'fixedGridRhoFastjetAll', eleRegs = dict( @@ -42,7 +42,7 @@ from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronFinalizer_cfi import lowPtGsfElectronFinalizer lowPtGsfElectrons = lowPtGsfElectronFinalizer.clone( previousGsfElectronsTag = "lowPtGsfElectronsPreRegression", - regressionConfig = _lowPtRegressionModifier, + regressionConfig = lowPtRegressionModifier, ) from Configuration.ProcessModifiers.run2_miniAOD_UL_cff import run2_miniAOD_UL diff --git a/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc b/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc index b09f7037b16cf..dacbcffc197f1 100644 --- a/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc +++ b/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc @@ -3,6 +3,7 @@ #include "DataFormats/GsfTrackReco/interface/GsfTrack.h" #include "DataFormats/ParticleFlowReco/interface/PFCluster.h" #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h" +#include "DataFormats/PatCandidates/interface/Electron.h" #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "TVector3.h" @@ -183,7 +184,7 @@ namespace lowptgsfeleid { // GSF tracks if (ele.core().isNonnull()) { - reco::GsfTrackRef gsf = ele.core()->gsfTrack(); + reco::GsfTrackRef gsf = ele.gsfTrack(); if (gsf.isNonnull()) { gsf_mode_p = gsf->pMode(); eid_gsf_nhits = (float)gsf->found(); @@ -198,7 +199,7 @@ namespace lowptgsfeleid { // Super clusters if (ele.core().isNonnull()) { - reco::SuperClusterRef sc = ele.core()->superCluster(); + reco::SuperClusterRef sc = ele.superCluster(); if (sc.isNonnull()) { eid_sc_E = sc->energy(); eid_sc_eta = sc->eta(); @@ -230,9 +231,9 @@ namespace lowptgsfeleid { // Clusters if (ele.core().isNonnull()) { - reco::GsfTrackRef gsf = ele.core()->gsfTrack(); + reco::GsfTrackRef gsf = ele.gsfTrack(); if (gsf.isNonnull()) { - reco::SuperClusterRef sc = ele.core()->superCluster(); + reco::SuperClusterRef sc = ele.superCluster(); if (sc.isNonnull()) { // Propagate electron track to ECAL surface double mass2 = 0.000511 * 0.000511; diff --git a/RecoEgamma/EgammaTools/BuildFile.xml b/RecoEgamma/EgammaTools/BuildFile.xml index 4557f9cbd48e8..6a5f6b7815913 100644 --- a/RecoEgamma/EgammaTools/BuildFile.xml +++ b/RecoEgamma/EgammaTools/BuildFile.xml @@ -1,9 +1,13 @@ + + + + diff --git a/RecoEgamma/EgammaTools/interface/LowPtConversion.h b/RecoEgamma/EgammaTools/interface/LowPtConversion.h new file mode 100644 index 0000000000000..744bce5caff8f --- /dev/null +++ b/RecoEgamma/EgammaTools/interface/LowPtConversion.h @@ -0,0 +1,68 @@ +#ifndef RecoEgamma_EgammaTools_LowPtConversion_h +#define RecoEgamma_EgammaTools_LowPtConversion_h + +#include "CommonTools/Statistics/interface/ChiSquaredProbability.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/EgammaCandidates/interface/Conversion.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Common/interface/RefToBase.h" +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/TrackReco/interface/Track.h" + +class LowPtConversion { +public: + LowPtConversion() = default; + ~LowPtConversion() = default; + + bool wpOpen() const; // Matched to any conversion (without selections) + bool wpLoose() const; // Nancy's baseline selections for conversions + bool wpTight() const; // Nancy's selection for analysis of conversions + + void addUserVars(pat::Electron& ele) const; // adds minimal set of flags to electron userData + void addExtraUserVars(pat::Electron& ele) const; // adds all variables to electron userData + + bool match(const reco::BeamSpot& beamSpot, const reco::ConversionCollection& conversions, const pat::Electron& ele); + + static float mee(float ipx1, float ipy1, float ipz1, float ipx2, float ipy2, float ipz2); + +private: + // quality + bool valid_ = false; + float chi2prob_ = -1.; + bool quality_high_purity_ = false; + bool quality_high_efficiency_ = false; + + // tracks + uint ntracks_ = 0; + float min_trk_pt_ = -1.; + int ilead_ = -1; + int itrail_ = -1; + + // displacement + float l_xy_ = -1.; + float vtx_radius_ = -1.; + + // invariant mass + float mass_from_conv_ = -1.; + float mass_from_Pin_ = -1.; + float mass_before_fit_ = -1.; + float mass_after_fit_ = -1.; + + // hits before vertex + uint lead_nhits_before_vtx_ = 0; + uint trail_nhits_before_vtx_ = 0; + uint max_nhits_before_vtx_ = 0; + uint sum_nhits_before_vtx_ = 0; + int delta_expected_nhits_inner_ = 0; + + // opening angle + float delta_cot_from_Pin_ = -1.; + + // match? + bool matched_ = false; + edm::RefToBase matched_lead_; + edm::RefToBase matched_trail_; +}; + +#endif // RecoEgamma_EgammaTools_LowPtConversion_h diff --git a/RecoEgamma/EgammaTools/plugins/EGRegressionModifierV3.cc b/RecoEgamma/EgammaTools/plugins/EGRegressionModifierV3.cc index a7d4f91fec372..0c24f9f4298d2 100644 --- a/RecoEgamma/EgammaTools/plugins/EGRegressionModifierV3.cc +++ b/RecoEgamma/EgammaTools/plugins/EGRegressionModifierV3.cc @@ -109,8 +109,7 @@ void EGRegressionModifierV3::modifyObject(reco::GsfElectron& ele) const { return; // do not apply corrections in case of missing info (slimmed MiniAOD electrons) - if (!superClus->clusters().isAvailable()) - return; + bool rescaleDependentValues = superClus->clusters().isAvailable(); //check if fbrem is filled as its needed for E/p combination so abort if its set to the default value //this will be the case for <5 (or current cuts) for miniAOD electrons @@ -141,7 +140,7 @@ void EGRegressionModifierV3::modifyObject(reco::GsfElectron& ele) const { const float corrEnergy = (rawEnergy + rawESEnergy) * ecalMeanCorr; const float corrEnergyErr = corrEnergy * ecalSigma; - ele.setCorrectedEcalEnergy(corrEnergy); + ele.setCorrectedEcalEnergy(corrEnergy, rescaleDependentValues); ele.setCorrectedEcalEnergyError(corrEnergyErr); std::pair combEnergyAndErr = eleRegs_->epComb.combine(ele); diff --git a/RecoEgamma/EgammaTools/plugins/LowPtElectronsModifier.cc b/RecoEgamma/EgammaTools/plugins/LowPtElectronsModifier.cc new file mode 100644 index 0000000000000..22b55a300e21c --- /dev/null +++ b/RecoEgamma/EgammaTools/plugins/LowPtElectronsModifier.cc @@ -0,0 +1,82 @@ +#include "CommonTools/CandAlgos/interface/ModifyObjectValueBase.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "RecoEgamma/EgammaTools/interface/LowPtConversion.h" + +//////////////////////////////////////////////////////////////////////////////// +// +class LowPtElectronModifier : public ModifyObjectValueBase { +public: + LowPtElectronModifier(const edm::ParameterSet& conf, edm::ConsumesCollector&); + ~LowPtElectronModifier() override = default; + + void setEvent(const edm::Event&) final; + void setEventContent(const edm::EventSetup&) final; + + void modifyObject(pat::Electron& ele) const final; + +private: + const edm::EDGetTokenT convT_; + reco::ConversionCollection const* conv_ = nullptr; + const edm::EDGetTokenT beamSpotT_; + reco::BeamSpot const* beamSpot_ = nullptr; + const edm::EDGetTokenT verticesT_; + reco::VertexCollection const* vertices_ = nullptr; + bool extra_; +}; + +//////////////////////////////////////////////////////////////////////////////// +// +LowPtElectronModifier::LowPtElectronModifier(const edm::ParameterSet& conf, edm::ConsumesCollector& cc) + : ModifyObjectValueBase(conf), + convT_(cc.consumes(conf.getParameter("conversions"))), + conv_(), + beamSpotT_(cc.consumes(conf.getParameter("beamSpot"))), + beamSpot_(), + verticesT_(cc.consumes(conf.getParameter("vertices"))), + vertices_(), + extra_(conf.getParameter("addExtraUserVars")) { + ; +} + +//////////////////////////////////////////////////////////////////////////////// +// +void LowPtElectronModifier::setEvent(const edm::Event& iEvent) { + conv_ = &iEvent.get(convT_); + beamSpot_ = &iEvent.get(beamSpotT_); + vertices_ = &iEvent.get(verticesT_); +} + +//////////////////////////////////////////////////////////////////////////////// +// +void LowPtElectronModifier::setEventContent(const edm::EventSetup& iSetup) {} + +//////////////////////////////////////////////////////////////////////////////// +// +void LowPtElectronModifier::modifyObject(pat::Electron& ele) const { + // Embed Conversion info + LowPtConversion conv; + conv.match(*beamSpot_, *conv_, ele); + conv.addUserVars(ele); + if (extra_) { + conv.addExtraUserVars(ele); + } + // Set impact parameters + auto const& gsfTrack = *ele.gsfTrack(); + if (!vertices_->empty()) { + const reco::Vertex& pv = vertices_->front(); + ele.setDB(gsfTrack.dxy(pv.position()), + gsfTrack.dxyError(pv.position(), pv.covariance()), + pat::Electron::PV2D); // PV2D + ele.setDB(gsfTrack.dz(pv.position()), std::hypot(gsfTrack.dzError(), pv.zError()), + pat::Electron::PVDZ); // PVDZ + } + ele.setDB(gsfTrack.dxy(*beamSpot_), gsfTrack.dxyError(*beamSpot_), + pat::Electron::BS2D); // BS2D +} + +//////////////////////////////////////////////////////////////////////////////// +// +DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, LowPtElectronModifier, "LowPtElectronModifier"); diff --git a/RecoEgamma/EgammaTools/python/lowPtElectronModifier_cfi.py b/RecoEgamma/EgammaTools/python/lowPtElectronModifier_cfi.py new file mode 100644 index 0000000000000..72636d5f00c3a --- /dev/null +++ b/RecoEgamma/EgammaTools/python/lowPtElectronModifier_cfi.py @@ -0,0 +1,9 @@ +import FWCore.ParameterSet.Config as cms + +lowPtElectronModifier = cms.PSet( + modifierName = cms.string('LowPtElectronModifier'), + beamSpot = cms.InputTag('offlineBeamSpot'), + conversions = cms.InputTag('gsfTracksOpenConversions:gsfTracksOpenConversions'), + addExtraUserVars = cms.bool(True), + vertices = cms.InputTag("offlineSlimmedPrimaryVertices"), +) diff --git a/RecoEgamma/EgammaTools/src/LowPtConversion.cc b/RecoEgamma/EgammaTools/src/LowPtConversion.cc new file mode 100644 index 0000000000000..1862490d4a180 --- /dev/null +++ b/RecoEgamma/EgammaTools/src/LowPtConversion.cc @@ -0,0 +1,212 @@ +#include "RecoEgamma/EgammaTools/interface/LowPtConversion.h" + +//////////////////////////////////////////////////////////////////////////////// +// Matched to any conversion (without selections) +// +bool LowPtConversion::wpOpen() const { return matched_; } + +//////////////////////////////////////////////////////////////////////////////// +// Nancy's baseline selections for conversions +// Based on: https://github.com/CMSBParking/BParkingNANO/blob/b2664ed/BParkingNano/plugins/ConversionSelector.cc#L253-L300 +bool LowPtConversion::wpLoose() const { + return (wpOpen() && ntracks_ == 2 && valid_ && quality_high_purity_ && chi2prob_ > 0.0005); +} + +//////////////////////////////////////////////////////////////////////////////// +// Nancy's selection for analysis of conversions +// Based on: slide 20 of https://indico.cern.ch/event/814402/contributions/3401312/ +bool LowPtConversion::wpTight() const { + return (wpLoose() && sum_nhits_before_vtx_ <= 1 && l_xy_ > 0. && mass_from_conv_ > 0. && // sanity check + mass_from_conv_ < 0.05); +} + +//////////////////////////////////////////////////////////////////////////////// +// adds minimal set of flags to electron userData +void LowPtConversion::addUserVars(pat::Electron& ele) const { + ele.addUserInt("convOpen", matched_ ? 1 : 0); + ele.addUserInt("convLoose", wpLoose() ? 1 : 0); + ele.addUserInt("convTight", wpTight() ? 1 : 0); + ele.addUserInt("convLead", matched_lead_.isNonnull() ? 1 : 0); + ele.addUserInt("convTrail", matched_trail_.isNonnull() ? 1 : 0); + if (ele.hasUserInt("convExtra") == false) { + ele.addUserInt("convExtra", 0); + } +} + +//////////////////////////////////////////////////////////////////////////////// +// adds all variables to electron userData +void LowPtConversion::addExtraUserVars(pat::Electron& ele) const { + // Flag that indicates if extra variables are added to electron userData + ele.addUserInt("convExtra", 1, true); // overwrite + + // quality + ele.addUserInt("convValid", valid_ ? 1 : 0); + ele.addUserFloat("convChi2Prob", chi2prob_); + ele.addUserInt("convQualityHighPurity", quality_high_purity_ ? 1 : 0); + ele.addUserInt("convQualityHighEff", quality_high_efficiency_ ? 1 : 0); + + // tracks + ele.addUserInt("convTracksN", ntracks_); + ele.addUserFloat("convMinTrkPt", min_trk_pt_); + ele.addUserInt("convLeadIdx", ilead_); + ele.addUserInt("convTrailIdx", itrail_); + + // displacement + ele.addUserFloat("convLxy", l_xy_); + ele.addUserFloat("convVtxRadius", vtx_radius_); + + // invariant mass + ele.addUserFloat("convMass", mass_from_conv_); + ele.addUserFloat("convMassFromPin", mass_from_Pin_); + ele.addUserFloat("convMassBeforeFit", mass_before_fit_); + ele.addUserFloat("convMassAfterFit", mass_after_fit_); + + // hits before vertex + ele.addUserInt("convLeadNHitsBeforeVtx", lead_nhits_before_vtx_); + ele.addUserInt("convTrailNHitsBeforeVtx", trail_nhits_before_vtx_); + ele.addUserInt("convMaxNHitsBeforeVtx", max_nhits_before_vtx_); + ele.addUserInt("convSumNHitsBeforeVtx", sum_nhits_before_vtx_); + ele.addUserInt("convDeltaExpectedNHitsInner", delta_expected_nhits_inner_); + + // opening angle + ele.addUserFloat("convDeltaCotFromPin", delta_cot_from_Pin_); +} + +//////////////////////////////////////////////////////////////////////////////// +// +bool LowPtConversion::match(const reco::BeamSpot& beamSpot, + const reco::ConversionCollection& conversions, + const pat::Electron& ele) { + // Iterate through conversions and calculate quantities (requirement from Nancy) + for (const auto& conv : conversions) { + // Filter + if (conv.tracks().size() != 2) { + continue; + } + + // Quality + valid_ = conv.conversionVertex().isValid(); // (=true) + chi2prob_ = ChiSquaredProbability(conv.conversionVertex().chi2(), conv.conversionVertex().ndof()); // (<0.005) + quality_high_purity_ = conv.quality(reco::Conversion::highPurity); // (=true) + quality_high_efficiency_ = conv.quality(reco::Conversion::highEfficiency); // (none) + + // Tracks + ntracks_ = conv.tracks().size(); // (=2) + min_trk_pt_ = -1.; // (>0.5) + for (const auto& trk : conv.tracks()) { + if (trk.isNonnull() && trk.isAvailable() && (min_trk_pt_ < 0. || trk->pt() < min_trk_pt_)) { + min_trk_pt_ = trk->pt(); + } + } + ilead_ = -1; + itrail_ = -1; + if (conv.tracks().size() == 2) { + const edm::RefToBase& trk1 = conv.tracks().front(); + const edm::RefToBase& trk2 = conv.tracks().back(); + if (trk1.isNonnull() && trk1.isAvailable() && trk2.isNonnull() && trk2.isAvailable()) { + if (trk1->pt() > trk2->pt()) { + ilead_ = 0; + itrail_ = 1; + } else { + ilead_ = 1; + itrail_ = 0; + } + } + } + + // Transverse displacement (with respect to beamspot) and vertex radius + math::XYZVectorF p_refitted = conv.refittedPairMomentum(); + float dx = conv.conversionVertex().x() - beamSpot.x0(); + float dy = conv.conversionVertex().y() - beamSpot.y0(); + l_xy_ = (p_refitted.x() * dx + p_refitted.y() * dy) / p_refitted.rho(); + vtx_radius_ = sqrt(conv.conversionVertex().position().perp2()); // (1.5= 2 && ilead_ > -1 && itrail_ > -1) { + math::XYZVectorF lead_Pin = conv.tracksPin().at(ilead_); + math::XYZVectorF trail_Pin = conv.tracksPin().at(itrail_); + mass_from_Pin_ = mee(lead_Pin.x(), lead_Pin.y(), lead_Pin.z(), trail_Pin.x(), trail_Pin.y(), trail_Pin.z()); + // Opening angle + delta_cot_from_Pin_ = 1. / tan(trail_Pin.theta()) - 1. / tan(lead_Pin.theta()); + } + + // Invariant mass before fit to common vertex + if (conv.tracks().size() >= 2 && ilead_ > -1 && itrail_ > -1) { + auto lead_before_vtx_fit = conv.tracks().at(ilead_)->momentum(); + auto trail_before_vtx_fit = conv.tracks().at(itrail_)->momentum(); + mass_before_fit_ = mee(lead_before_vtx_fit.x(), + lead_before_vtx_fit.y(), + lead_before_vtx_fit.z(), + trail_before_vtx_fit.x(), + trail_before_vtx_fit.y(), + trail_before_vtx_fit.z()); + } + + // Invariant mass after the fit to common vertex + if (conv.conversionVertex().refittedTracks().size() >= 2 && ilead_ > -1 && itrail_ > -1) { + auto const& lead_after_vtx_fit = conv.conversionVertex().refittedTracks().at(ilead_); + auto const& trail_after_vtx_fit = conv.conversionVertex().refittedTracks().at(itrail_); + mass_after_fit_ = mee(lead_after_vtx_fit.px(), + lead_after_vtx_fit.py(), + lead_after_vtx_fit.pz(), + trail_after_vtx_fit.px(), + trail_after_vtx_fit.py(), + trail_after_vtx_fit.pz()); + // Difference in expeted hits + delta_expected_nhits_inner_ = + lead_after_vtx_fit.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) - + trail_after_vtx_fit.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS); + } + + // Hits prior to vertex + if (ilead_ > -1 && itrail_ > -1) { + auto const& nHits = conv.nHitsBeforeVtx(); + bool enoughTracks = nHits.size() > 1; + lead_nhits_before_vtx_ = enoughTracks ? nHits.at(ilead_) : 0; + trail_nhits_before_vtx_ = enoughTracks ? nHits.at(itrail_) : 0; + max_nhits_before_vtx_ = enoughTracks ? std::max(nHits[0], nHits[1]) : 0; + sum_nhits_before_vtx_ = enoughTracks ? nHits[0] + nHits[1] : 0; + } + + // Attempt to match conversion track to electron + for (uint itrk = 0; itrk < conv.tracks().size(); ++itrk) { + const edm::RefToBase trk = conv.tracks()[itrk]; + if (trk.isNull()) { + continue; + } + reco::GsfTrackRef ref = ele.core()->gsfTrack(); + reco::GsfTrackRef gsf = ele.gsfTrack(); + if (gsf.isNull()) { + continue; + } + if (ref.id() == trk.id() && ref.key() == trk.key()) { + matched_ = true; + if (static_cast(itrk) == ilead_) { + matched_lead_ = trk; + } + if (static_cast(itrk) == itrail_) { + matched_trail_ = trk; + } + } + } // track loop + } // conversions loop + + return matched_; +} + +//////////////////////////////////////////////////////////////////////////////// +// +float LowPtConversion::mee(float px1, float py1, float pz1, float px2, float py2, float pz2) { + const float m = 0.000511; + const float px = px1 + px2; + const float py = py1 + py2; + const float pz = pz1 + pz2; + const float p1 = px1 * px1 + py1 * py1 + pz1 * pz1; + const float p2 = px2 * px2 + py2 * py2 + pz2 * pz2; + const float e = sqrt(p1 + m * m) + sqrt(p2 + m * m); + const float mass = (e * e - px * px - py * py - pz * pz); + return mass > 0. ? sqrt(mass) : -1.; +}