diff --git a/CommonTools/MVAUtils/BuildFile.xml b/CommonTools/MVAUtils/BuildFile.xml index 9c7589c0a58b1..c6e15e70055d6 100644 --- a/CommonTools/MVAUtils/BuildFile.xml +++ b/CommonTools/MVAUtils/BuildFile.xml @@ -1,3 +1,4 @@ + diff --git a/CommonTools/MVAUtils/bin/BuildFile.xml b/CommonTools/MVAUtils/bin/BuildFile.xml new file mode 100644 index 0000000000000..147f2d969912d --- /dev/null +++ b/CommonTools/MVAUtils/bin/BuildFile.xml @@ -0,0 +1,4 @@ + + + + diff --git a/CommonTools/MVAUtils/bin/convertXMLToGBRForestROOT.cc b/CommonTools/MVAUtils/bin/convertXMLToGBRForestROOT.cc new file mode 100644 index 0000000000000..3c7c9f5805fd7 --- /dev/null +++ b/CommonTools/MVAUtils/bin/convertXMLToGBRForestROOT.cc @@ -0,0 +1,34 @@ +#include "CommonTools/MVAUtils/interface/GBRForestTools.h" + +#include "TFile.h" + +#include +#include + +int main(int argc, char **argv) { + if (argc != 3) { + std::cout << "Please pass a (gzipped) BDT weight file and a name for the output ROOT file." << std::endl; + return 1; + } + + char *inputFileName = argv[1]; + char *outputFileName = argv[2]; + + if (!std::filesystem::exists(inputFileName)) { + std::cout << "Input file " << inputFileName << " does not exists." << std::endl; + return 1; + } + + if (std::filesystem::exists(outputFileName)) { + std::cout << "Output file " << outputFileName << " already exists." << std::endl; + return 1; + } + + auto gbrForest = createGBRForest(inputFileName); + std::cout << "Read GBRForest " << inputFileName << " successfully." << std::endl; + + TFile{outputFileName, "RECREATE"}.WriteObject(gbrForest.get(), "gbrForest"); + std::cout << "GBRForest written to " << outputFileName << " successfully." << std::endl; + + return 0; +} diff --git a/CommonTools/MVAUtils/src/GBRForestTools.cc b/CommonTools/MVAUtils/src/GBRForestTools.cc index a277c1f63022e..71e39cd60f23e 100644 --- a/CommonTools/MVAUtils/src/GBRForestTools.cc +++ b/CommonTools/MVAUtils/src/GBRForestTools.cc @@ -3,11 +3,14 @@ #include "FWCore/ParameterSet/interface/FileInPath.h" #include "FWCore/Utilities/interface/Exception.h" +#include "TFile.h" + #include #include #include #include #include +#include namespace { @@ -117,6 +120,16 @@ namespace { } std::unique_ptr init(const std::string& weightsFileFullPath, std::vector& varNames) { + // + // Load weights file, for ROOT file + // + if (reco::details::hasEnding(weightsFileFullPath, ".root")) { + TFile gbrForestFile(weightsFileFullPath.c_str()); + std::unique_ptr up(gbrForestFile.Get("gbrForest")); + gbrForestFile.Close("nodelete"); + return up; + } + // // Load weights file, for gzipped or raw xml file // diff --git a/DataFormats/PatCandidates/src/classes_def_objects.xml b/DataFormats/PatCandidates/src/classes_def_objects.xml index 304414153857a..5769f1247b8ca 100644 --- a/DataFormats/PatCandidates/src/classes_def_objects.xml +++ b/DataFormats/PatCandidates/src/classes_def_objects.xml @@ -484,6 +484,8 @@ + + diff --git a/PhysicsTools/PatAlgos/plugins/BuildFile.xml b/PhysicsTools/PatAlgos/plugins/BuildFile.xml index 957b468c99b10..d17868f1629ff 100644 --- a/PhysicsTools/PatAlgos/plugins/BuildFile.xml +++ b/PhysicsTools/PatAlgos/plugins/BuildFile.xml @@ -9,6 +9,7 @@ + diff --git a/PhysicsTools/PatAlgos/plugins/LowPtGSFToPackedCandidateLinker.cc b/PhysicsTools/PatAlgos/plugins/LowPtGSFToPackedCandidateLinker.cc index 3ca2e5c385b12..a82b2fd3ac561 100644 --- a/PhysicsTools/PatAlgos/plugins/LowPtGSFToPackedCandidateLinker.cc +++ b/PhysicsTools/PatAlgos/plugins/LowPtGSFToPackedCandidateLinker.cc @@ -1,10 +1,11 @@ #include #include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/Common/interface/RefToPtr.h" #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" -#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h" #include "DataFormats/GsfTrackReco/interface/GsfTrack.h" #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h" +#include "DataFormats/PatCandidates/interface/Electron.h" #include "DataFormats/PatCandidates/interface/PackedCandidate.h" #include "DataFormats/Common/interface/Association.h" #include "FWCore/Framework/interface/global/EDProducer.h" @@ -19,6 +20,9 @@ #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h" #include "FWCore/Framework/interface/MakerMacros.h" +typedef edm::Ptr PackedCandidatePtr; +typedef std::vector PackedCandidatePtrCollection; + class LowPtGSFToPackedCandidateLinker : public edm::global::EDProducer<> { public: explicit LowPtGSFToPackedCandidateLinker(const edm::ParameterSet&); @@ -36,6 +40,7 @@ class LowPtGSFToPackedCandidateLinker : public edm::global::EDProducer<> { const edm::EDGetTokenT > lost2trk_; const edm::EDGetTokenT > gsf2trk_; const edm::EDGetTokenT > gsftracks_; + const edm::EDGetTokenT > electrons_; }; LowPtGSFToPackedCandidateLinker::LowPtGSFToPackedCandidateLinker(const edm::ParameterSet& iConfig) @@ -48,37 +53,26 @@ LowPtGSFToPackedCandidateLinker::LowPtGSFToPackedCandidateLinker(const edm::Para lost2trk_{consumes >( iConfig.getParameter("lostTracks"))}, gsf2trk_{consumes >(iConfig.getParameter("gsfToTrack"))}, - gsftracks_{consumes >(iConfig.getParameter("gsfTracks"))} { - produces >("packedCandidates"); - produces >("lostTracks"); + gsftracks_{consumes >(iConfig.getParameter("gsfTracks"))}, + electrons_{consumes >(iConfig.getParameter("electrons"))} { + produces >("gsf2packed"); + produces >("gsf2lost"); + produces >("ele2packed"); + produces >("ele2lost"); } LowPtGSFToPackedCandidateLinker::~LowPtGSFToPackedCandidateLinker() {} void LowPtGSFToPackedCandidateLinker::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { - edm::Handle pfcands; - iEvent.getByToken(pfcands_, pfcands); - - edm::Handle packed; - iEvent.getByToken(packed_, packed); - - edm::Handle lost_tracks; - iEvent.getByToken(lost_tracks_, lost_tracks); - - edm::Handle > pf2packed; - iEvent.getByToken(pf2packed_, pf2packed); - - edm::Handle > lost2trk_assoc; - iEvent.getByToken(lost2trk_, lost2trk_assoc); - - edm::Handle > gsftracks; - iEvent.getByToken(gsftracks_, gsftracks); - - edm::Handle tracks; - iEvent.getByToken(tracks_, tracks); - - edm::Handle > gsf2trk; - iEvent.getByToken(gsf2trk_, gsf2trk); + auto pfcands = iEvent.getHandle(pfcands_); + auto packed = iEvent.getHandle(packed_); + auto lost_tracks = iEvent.getHandle(lost_tracks_); + auto pf2packed = iEvent.getHandle(pf2packed_); + auto lost2trk_assoc = iEvent.getHandle(lost2trk_); + auto gsftracks = iEvent.getHandle(gsftracks_); + auto tracks = iEvent.getHandle(tracks_); + auto gsf2trk = iEvent.getHandle(gsf2trk_); + auto electrons = iEvent.getHandle(electrons_); // collection sizes, for reference const size_t npf = pfcands->size(); @@ -86,6 +80,7 @@ void LowPtGSFToPackedCandidateLinker::produce(edm::StreamID, edm::Event& iEvent, const size_t nlost = lost_tracks->size(); const size_t ntracks = tracks->size(); const size_t ngsf = gsftracks->size(); + const size_t nele = electrons->size(); //store index mapping in vectors for easy and fast access std::vector trk2packed(ntracks, npacked); @@ -94,6 +89,8 @@ void LowPtGSFToPackedCandidateLinker::produce(edm::StreamID, edm::Event& iEvent, //store auxiliary mappings for association std::vector gsf2pack(ngsf, -1); std::vector gsf2lost(ngsf, -1); + PackedCandidatePtrCollection ele2packedptr(nele, PackedCandidatePtr(packed, -1)); + PackedCandidatePtrCollection ele2lostptr(nele, PackedCandidatePtr(lost_tracks, -1)); //electrons will never store their track (they store the Gsf track) //map PackedPF <--> Track @@ -103,8 +100,7 @@ void LowPtGSFToPackedCandidateLinker::produce(edm::StreamID, edm::Event& iEvent, auto packed_ref = (*pf2packed)[pf_ref]; if (cand.charge() && packed_ref.isNonnull() && cand.trackRef().isNonnull() && cand.trackRef().id() == tracks.id()) { size_t trkid = cand.trackRef().index(); - size_t packid = packed_ref.index(); - trk2packed[trkid] = packid; + trk2packed[trkid] = packed_ref.index(); } } @@ -113,22 +109,20 @@ void LowPtGSFToPackedCandidateLinker::produce(edm::StreamID, edm::Event& iEvent, reco::TrackRef key(tracks, itrk); pat::PackedCandidateRef lostTrack = (*lost2trk_assoc)[key]; if (lostTrack.isNonnull()) { - size_t ilost = lostTrack.index(); //assumes that LostTracks are all made from the same track collection - trk2lost[itrk] = ilost; + trk2lost[itrk] = lostTrack.index(); // assumes that LostTracks are all made from the same track collection } } //map Track --> GSF and fill GSF --> PackedCandidates and GSF --> Lost associations for (unsigned int igsf = 0; igsf < ngsf; ++igsf) { - reco::GsfTrackRef gref(gsftracks, igsf); - reco::TrackRef trk = (*gsf2trk)[gref]; - if (trk.id() != tracks.id()) { + reco::GsfTrackRef gsf_ref(gsftracks, igsf); + reco::TrackRef trk_ref = (*gsf2trk)[gsf_ref]; + if (trk_ref.id() != tracks.id()) { throw cms::Exception( "WrongCollection", "The reco::Track collection used to match against the GSF Tracks was not used to produce such tracks"); } - size_t trkid = trk.index(); - + size_t trkid = trk_ref.index(); if (trk2packed[trkid] != npacked) { gsf2pack[igsf] = trk2packed[trkid]; } @@ -137,18 +131,51 @@ void LowPtGSFToPackedCandidateLinker::produce(edm::StreamID, edm::Event& iEvent, } } + //map Electron-->pat::PFCandidatePtr via Electron-->GsfTrack-->Track and Track-->pat::PFCandidatePtr + for (unsigned int iele = 0; iele < nele; ++iele) { + auto const& ele = (*electrons)[iele]; + reco::GsfTrackRef gsf_ref = ele.core()->gsfTrack(); + reco::TrackRef trk_ref = (*gsf2trk)[gsf_ref]; + if (trk_ref.id() != tracks.id()) { + throw cms::Exception( + "WrongCollection", + "The reco::Track collection used to match against the GSF Tracks was not used to produce such tracks"); + } + size_t trkid = trk_ref.index(); + auto packedIdx = trk2packed[trkid]; + if (packedIdx != npacked) { + ele2packedptr[iele] = PackedCandidatePtr(packed, packedIdx); + } + auto lostIdx = trk2lost[trkid]; + if (lostIdx != nlost) { + ele2lostptr[iele] = PackedCandidatePtr(lost_tracks, lostIdx); + } + } + // create output collections from the mappings auto assoc_gsf2pack = std::make_unique >(packed); edm::Association::Filler gsf2pack_filler(*assoc_gsf2pack); gsf2pack_filler.insert(gsftracks, gsf2pack.begin(), gsf2pack.end()); gsf2pack_filler.fill(); - iEvent.put(std::move(assoc_gsf2pack), "packedCandidates"); + iEvent.put(std::move(assoc_gsf2pack), "gsf2packed"); auto assoc_gsf2lost = std::make_unique >(lost_tracks); edm::Association::Filler gsf2lost_filler(*assoc_gsf2lost); gsf2lost_filler.insert(gsftracks, gsf2lost.begin(), gsf2lost.end()); gsf2lost_filler.fill(); - iEvent.put(std::move(assoc_gsf2lost), "lostTracks"); + iEvent.put(std::move(assoc_gsf2lost), "gsf2lost"); + + auto map_ele2packedptr = std::make_unique >(); + edm::ValueMap::Filler ele2packedptr_filler(*map_ele2packedptr); + ele2packedptr_filler.insert(electrons, ele2packedptr.begin(), ele2packedptr.end()); + ele2packedptr_filler.fill(); + iEvent.put(std::move(map_ele2packedptr), "ele2packed"); + + auto map_ele2lostptr = std::make_unique >(); + edm::ValueMap::Filler ele2lostptr_filler(*map_ele2lostptr); + ele2lostptr_filler.insert(electrons, ele2lostptr.begin(), ele2lostptr.end()); + ele2lostptr_filler.fill(); + iEvent.put(std::move(map_ele2lostptr), "ele2lost"); } void LowPtGSFToPackedCandidateLinker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { @@ -159,6 +186,7 @@ void LowPtGSFToPackedCandidateLinker::fillDescriptions(edm::ConfigurationDescrip desc.add("tracks", edm::InputTag("generalTracks")); desc.add("gsfToTrack", edm::InputTag("lowPtGsfToTrackLinks")); desc.add("gsfTracks", edm::InputTag("lowPtGsfEleGsfTracks")); + desc.add("electrons", edm::InputTag("selectedPatLowPtElectrons")); descriptions.add("lowPtGsfLinksDefault", desc); } diff --git a/PhysicsTools/PatAlgos/python/producersLayer1/lowPtElectronProducer_cff.py b/PhysicsTools/PatAlgos/python/producersLayer1/lowPtElectronProducer_cff.py index b64a33f7a3e3f..f6f745d01b956 100644 --- a/PhysicsTools/PatAlgos/python/producersLayer1/lowPtElectronProducer_cff.py +++ b/PhysicsTools/PatAlgos/python/producersLayer1/lowPtElectronProducer_cff.py @@ -11,42 +11,49 @@ ) patLowPtElectrons = patElectrons.clone( - # input collections + + # Input collection electronSource = sourceElectrons, - genParticleMatch = cms.InputTag("lowPtElectronMatch"), - # overrides - addElectronID = cms.bool(False), - addGenMatch = cms.bool(True), - addMVAVariables = cms.bool(False), - addPFClusterIso = cms.bool(False), - electronIDSources = cms.PSet(), - computeMiniIso = cms.bool(False), - isoDeposits = cms.PSet(), - isolationValues = cms.PSet(), - isolationValuesNoPFId = cms.PSet(), - miniIsoParamsB = cms.vdouble(), - miniIsoParamsE = cms.vdouble(), - usePfCandidateMultiMap = cms.bool(False), - # embedding - embedBasicClusters = cms.bool(False), - embedGenMatch = cms.bool(False), - embedGsfElectronCore = cms.bool(False), - embedGsfTrack = cms.bool(False), - embedHighLevelSelection = cms.bool(False), - embedPFCandidate = cms.bool(False), - embedPflowBasicClusters = cms.bool(False), - embedPflowPreshowerClusters = cms.bool(False), - embedPflowSuperCluster = cms.bool(False), - embedPreshowerClusters = cms.bool(False), - embedRecHits = cms.bool(False), - embedSeedCluster = cms.bool(False), - embedSuperCluster = cms.bool(False), - embedTrack = cms.bool(True), - ) + + # MC matching + genParticleMatch = "lowPtElectronMatch", + + # Electron ID + addElectronID = True, + electronIDSources = dict( + unbiased = cms.InputTag("rekeyLowPtGsfElectronSeedValueMaps:unbiased"), + ptbiased = cms.InputTag("rekeyLowPtGsfElectronSeedValueMaps:ptbiased"), + ID = cms.InputTag("lowPtGsfElectronID"), + ), + + # Embedding of RECO/AOD items + + # Embedding of RECO/AOD items + embedTrack = True, + embedGsfElectronCore = True, + embedGsfTrack = True, + embedSuperCluster = True, + embedSeedCluster = True, + embedBasicClusters = True, + embedPreshowerClusters = False, + embedRecHits = False, + embedPflowSuperCluster = False, + embedPflowBasicClusters = False, + embedPflowPreshowerClusters = False, + embedPFCandidate = False, + + # Miscellaneous flags + addMVAVariables = False, + embedHighLevelSelection = False, + isoDeposits = cms.PSet(), + isolationValues = cms.PSet(), + isolationValuesNoPFId = cms.PSet(), + +) makePatLowPtElectronsTask = cms.Task( lowPtElectronMatch, - patLowPtElectrons + patLowPtElectrons, ) makePatLowPtElectrons = cms.Sequence(makePatLowPtElectronsTask) @@ -61,3 +68,12 @@ electronSource = "gedGsfElectrons", genParticleMatch = "electronMatch" ) + +# Schedule rekeying of seed BDT ValueMaps by reco::GsfElectron for run2_miniAOD_UL +from Configuration.ProcessModifiers.run2_miniAOD_UL_cff import run2_miniAOD_UL +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronSeedValueMaps_cff import rekeyLowPtGsfElectronSeedValueMaps +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronID_cfi import lowPtGsfElectronID +_makePatLowPtElectronsTask = makePatLowPtElectronsTask.copy() +_makePatLowPtElectronsTask.add(rekeyLowPtGsfElectronSeedValueMaps) +_makePatLowPtElectronsTask.add(lowPtGsfElectronID) +run2_miniAOD_UL.toReplaceWith(makePatLowPtElectronsTask,_makePatLowPtElectronsTask) diff --git a/PhysicsTools/PatAlgos/python/producersLayer1/patCandidates_cff.py b/PhysicsTools/PatAlgos/python/producersLayer1/patCandidates_cff.py index c8e9e1177acd3..6f40251875909 100644 --- a/PhysicsTools/PatAlgos/python/producersLayer1/patCandidates_cff.py +++ b/PhysicsTools/PatAlgos/python/producersLayer1/patCandidates_cff.py @@ -38,8 +38,16 @@ _patCandidatesTask = patCandidatesTask.copy() _patCandidatesTask.remove(makePatOOTPhotonsTask) from Configuration.ProcessModifiers.pp_on_AA_cff import pp_on_AA -pp_on_AA.toReplaceWith(patCandidatesTask, _patCandidatesTask) +pp_on_AA.toReplaceWith(patCandidatesTask, _patCandidatesTask) pp_on_AA.toModify(patCandidateSummary.candidates, func = lambda list: list.remove(cms.InputTag("patOOTPhotons")) ) +from Configuration.Eras.Modifier_run2_miniAOD_94XFall17_cff import run2_miniAOD_94XFall17 +from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy +_mAOD = (run2_miniAOD_94XFall17 | run2_miniAOD_80XLegacy) +(pp_on_AA | _mAOD).toReplaceWith(patCandidatesTask, + patCandidatesTask.copyAndExclude([makePatLowPtElectronsTask])) +(pp_on_AA | _mAOD).toModify(patCandidateSummary.candidates, + func = lambda list: list.remove(cms.InputTag("patLowPtElectrons")) ) + patCandidates = cms.Sequence(patCandidateSummary, patCandidatesTask) diff --git a/PhysicsTools/PatAlgos/python/selectionLayer1/lowPtElectronSelector_cfi.py b/PhysicsTools/PatAlgos/python/selectionLayer1/lowPtElectronSelector_cfi.py index 17b510bf0b51b..7ea6b50afd8b8 100644 --- a/PhysicsTools/PatAlgos/python/selectionLayer1/lowPtElectronSelector_cfi.py +++ b/PhysicsTools/PatAlgos/python/selectionLayer1/lowPtElectronSelector_cfi.py @@ -6,9 +6,13 @@ # selectedPatLowPtElectrons = cms.EDFilter("PATElectronSelector", src = cms.InputTag("patLowPtElectrons"), - cut = cms.string("") + cut = cms.string("pt>1. && electronID('ID')>1.5"), ) +# Modifier for bParking (fully open selection) +from Configuration.Eras.Modifier_bParking_cff import bParking +bParking.toModify(selectedPatLowPtElectrons,cut = "") + # Modifiers for legacy AOD from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy from Configuration.Eras.Modifier_run2_miniAOD_94XFall17_cff import run2_miniAOD_94XFall17 diff --git a/PhysicsTools/PatAlgos/python/selectionLayer1/selectedPatCandidates_cff.py b/PhysicsTools/PatAlgos/python/selectionLayer1/selectedPatCandidates_cff.py index 3431fa77972cd..50381daef2318 100644 --- a/PhysicsTools/PatAlgos/python/selectionLayer1/selectedPatCandidates_cff.py +++ b/PhysicsTools/PatAlgos/python/selectionLayer1/selectedPatCandidates_cff.py @@ -38,3 +38,11 @@ from Configuration.ProcessModifiers.pp_on_AA_cff import pp_on_AA pp_on_AA.toReplaceWith(selectedPatCandidatesTask, selectedPatCandidatesTask.copyAndExclude([selectedPatOOTPhotons])) pp_on_AA.toModify(selectedPatCandidateSummary.candidates, func = lambda list: list.remove(cms.InputTag("selectedPatOOTPhotons")) ) + +from Configuration.Eras.Modifier_run2_miniAOD_94XFall17_cff import run2_miniAOD_94XFall17 +from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy +_mAOD = (run2_miniAOD_94XFall17 | run2_miniAOD_80XLegacy) +(pp_on_AA | _mAOD).toReplaceWith(selectedPatCandidatesTask, + selectedPatCandidatesTask.copyAndExclude([selectedPatLowPtElectrons])) +(pp_on_AA | _mAOD).toModify(selectedPatCandidateSummary.candidates, + func = lambda list: list.remove(cms.InputTag("selectedPatLowPtElectrons")) ) diff --git a/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py b/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py index 92bae6c6e2d46..31b44e0e4b561 100644 --- a/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py @@ -89,6 +89,9 @@ 'keep recoTracks_displacedStandAloneMuons__*', # L1 prefiring weights 'keep *_prefiringweight_*_*', + # patLowPtElectrons + 'keep *_slimmedLowPtElectrons_*_*', + 'keep *_gsfTracksOpenConversions_*_*', ) ) @@ -115,20 +118,6 @@ ) ) -# Only add low pT electrons for bParking era -from Configuration.Eras.Modifier_bParking_cff import bParking -_bParking_extraCommands = ['keep *_slimmedLowPtElectrons_*_*', - 'keep recoGsfElectronCores_lowPtGsfElectronCores_*_*', - 'keep recoSuperClusters_lowPtGsfElectronSuperClusters_*_*', - 'keep recoCaloClusters_lowPtGsfElectronSuperClusters_*_*', - 'keep recoGsfTracks_lowPtGsfEleGsfTracks_*_*', - 'keep floatedmValueMap_lowPtGsfElectronSeedValueMaps_*_*', - 'keep floatedmValueMap_lowPtGsfElectronID_*_*', - 'keep *_lowPtGsfLinks_*_*', - 'keep *_gsfTracksOpenConversions_*_*', - ] -bParking.toModify(MicroEventContent, outputCommands = MicroEventContent.outputCommands + _bParking_extraCommands) - # --- Only for 2018 data & MC _run2_HCAL_2018_extraCommands = ["keep *_packedPFCandidates_hcalDepthEnergyFractions_*"] from Configuration.Eras.Modifier_run2_HCAL_2018_cff import run2_HCAL_2018 diff --git a/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cff.py b/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cff.py new file mode 100644 index 0000000000000..73e67d9d1a08c --- /dev/null +++ b/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cff.py @@ -0,0 +1,11 @@ +import FWCore.ParameterSet.Config as cms + +from PhysicsTools.PatAlgos.slimming.slimmedLowPtElectrons_cfi import slimmedLowPtElectrons +from PhysicsTools.PatAlgos.slimming.lowPtGsfLinks_cfi import lowPtGsfLinks +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronID_cfi import lowPtGsfElectronID + +# Task +slimmedLowPtElectronsTask = cms.Task( + lowPtGsfLinks, + slimmedLowPtElectrons, +) diff --git a/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py b/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py index 1095a37703ca2..684512a4ce0fd 100644 --- a/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py +++ b/PhysicsTools/PatAlgos/python/slimming/slimmedLowPtElectrons_cfi.py @@ -20,7 +20,26 @@ saveNonZSClusterShapes = cms.string("1"), # save additional user floats: (sigmaIetaIeta,sigmaIphiIphi,sigmaIetaIphi,r9,e1x5_over_e5x5)_NoZS reducedBarrelRecHitCollection = cms.InputTag("reducedEcalRecHitsEB"), reducedEndcapRecHitCollection = cms.InputTag("reducedEcalRecHitsEE"), - modifyElectrons = cms.bool(False), - modifierConfig = cms.PSet( modifications = cms.VPSet() ) + modifyElectrons = cms.bool(True), + modifierConfig = cms.PSet( + modifications = cms.VPSet( + cms.PSet( + electron_config = cms.PSet( + ele2packed = cms.InputTag("lowPtGsfLinks:ele2packed"), + electronSrc = cms.InputTag("selectedPatLowPtElectrons"), + ), + modifierName = cms.string('EGExtraInfoModifierFromPackedCandPtrValueMaps'), + photon_config = cms.PSet() + ), + cms.PSet( + electron_config = cms.PSet( + ele2lost = cms.InputTag("lowPtGsfLinks:ele2lost"), + electronSrc = cms.InputTag("selectedPatLowPtElectrons"), + ), + modifierName = cms.string('EGExtraInfoModifierFromPackedCandPtrValueMaps'), + photon_config = cms.PSet() + ), + ) + ) ) diff --git a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py index 4b97dbfae42ff..454ed98e007c7 100644 --- a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py @@ -13,8 +13,7 @@ from PhysicsTools.PatAlgos.slimming.slimmedCaloJets_cfi import * from PhysicsTools.PatAlgos.slimming.slimmedGenJets_cfi import * from PhysicsTools.PatAlgos.slimming.slimmedElectrons_cfi import * -from PhysicsTools.PatAlgos.slimming.slimmedLowPtElectrons_cfi import * -from PhysicsTools.PatAlgos.slimming.lowPtGsfLinks_cfi import * +from PhysicsTools.PatAlgos.slimming.slimmedLowPtElectrons_cff import * from PhysicsTools.PatAlgos.slimming.slimmedTrackExtras_cff import * from PhysicsTools.PatAlgos.slimming.slimmedMuons_cfi import * from PhysicsTools.PatAlgos.slimming.slimmedPhotons_cfi import * @@ -45,8 +44,7 @@ slimmedGenJets, slimmedGenJetsAK8, slimmedElectrons, - slimmedLowPtElectrons, - lowPtGsfLinks, + slimmedLowPtElectronsTask, slimmedMuonTrackExtras, slimmedMuons, slimmedPhotons, @@ -65,6 +63,13 @@ from Configuration.ProcessModifiers.pp_on_AA_cff import pp_on_AA pp_on_AA.toReplaceWith(slimmingTask, slimmingTask.copyAndExclude([slimmedOOTPhotons])) + +from Configuration.Eras.Modifier_run2_miniAOD_94XFall17_cff import run2_miniAOD_94XFall17 +from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy +_mAOD = (run2_miniAOD_94XFall17 | run2_miniAOD_80XLegacy) +(pp_on_AA | _mAOD).toReplaceWith(slimmingTask, + slimmingTask.copyAndExclude([slimmedLowPtElectronsTask])) + from PhysicsTools.PatAlgos.slimming.hiPixelTracks_cfi import hiPixelTracks from RecoHI.HiEvtPlaneAlgos.HiEvtPlane_cfi import hiEvtPlane from RecoHI.HiEvtPlaneAlgos.hiEvtPlaneFlat_cfi import hiEvtPlaneFlat diff --git a/PhysicsTools/PatAlgos/test/IntegrationTest_cfg.py b/PhysicsTools/PatAlgos/test/IntegrationTest_cfg.py index 4808105fcf8fa..fa1bd5897bfd7 100644 --- a/PhysicsTools/PatAlgos/test/IntegrationTest_cfg.py +++ b/PhysicsTools/PatAlgos/test/IntegrationTest_cfg.py @@ -22,6 +22,7 @@ process.selectedPatCandidates ) +process.patLowPtElectrons.addElectronID = False process.patLowPtElectrons.electronSource = "gedGsfElectrons" process.patLowPtElectrons.genParticleMatch = "electronMatch" process.selectedPatLowPtElectrons.cut = "pt>99999." diff --git a/PhysicsTools/PatAlgos/test/patTuple_addTriggerInfo_cfg.py b/PhysicsTools/PatAlgos/test/patTuple_addTriggerInfo_cfg.py index 313f0465bb950..9e7c2fe6ffa7e 100644 --- a/PhysicsTools/PatAlgos/test/patTuple_addTriggerInfo_cfg.py +++ b/PhysicsTools/PatAlgos/test/patTuple_addTriggerInfo_cfg.py @@ -18,6 +18,7 @@ process.selectedPatCandidates ) +process.patLowPtElectrons.addElectronID = False process.patLowPtElectrons.electronSource = "gedGsfElectrons" process.patLowPtElectrons.genParticleMatch = "electronMatch" process.selectedPatLowPtElectrons.cut = "pt>99999." diff --git a/RecoEgamma/Configuration/python/RecoEgamma_EventContent_cff.py b/RecoEgamma/Configuration/python/RecoEgamma_EventContent_cff.py index 42358054685bf..f14d59c228b0f 100644 --- a/RecoEgamma/Configuration/python/RecoEgamma_EventContent_cff.py +++ b/RecoEgamma/Configuration/python/RecoEgamma_EventContent_cff.py @@ -46,7 +46,8 @@ 'keep *_lowPtGsfToTrackLinks_*_*', 'keep recoSuperClusters_lowPtGsfElectronSuperClusters_*_*', 'keep floatedmValueMap_lowPtGsfElectronSeedValueMaps_*_*', - 'keep floatedmValueMap_lowPtGsfElectronID_*_*') + 'keep floatedmValueMap_rekeyLowPtGsfElectronSeedValueMaps_*_*', + 'keep floatedmValueMap_lowPtGsfElectronID_*_*') ) # mods for HGCAL _phase2_hgcal_RecoEgamma_tokeep = [ 'keep *_ecalDrivenGsfElectronCores_*_*', diff --git a/RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h b/RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h index f2d64aff698ba..89e81b2b60d17 100644 --- a/RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h +++ b/RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h @@ -23,11 +23,29 @@ namespace lowptgsfeleseed { namespace lowptgsfeleid { - std::vector features(edm::Ptr const& ele, float rho, float unbiased); - - std::vector features(edm::Ref > const& ele, float rho, float unbiased); - - std::vector features(edm::Ref > const& ele, float rho, float unbiased); + // feature list for new model (2019Sept15) + std::vector features_V1(reco::GsfElectron const& ele, float rho, float unbiased, float field_z); + + // feature list for original models (2019Aug07 and earlier) + std::vector features_V0(reco::GsfElectron const& ele, float rho, float unbiased); + + // Find most energetic clusters + void findEnergeticClusters(reco::SuperCluster const&, int&, float&, float&, int&, int&); + + // Track-cluster matching for most energetic clusters + void trackClusterMatching(reco::SuperCluster const&, + reco::GsfTrack const&, + bool const&, + GlobalPoint const&, + float&, + float&, + float&, + float&, + float&, + float&, + float&, + float&, + float&); } // namespace lowptgsfeleid diff --git a/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml b/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml index dd6603c1f6a81..12d72e6b58257 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml +++ b/RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml @@ -4,6 +4,7 @@ + diff --git a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc index 8eba1d052f8d0..6b2d8cf39a3c2 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc +++ b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc @@ -11,11 +11,14 @@ #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include "CommonTools/MVAUtils/interface/GBRForestTools.h" #include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/PatCandidates/interface/Electron.h" #include "DataFormats/EgammaReco/interface/SuperCluster.h" #include "DataFormats/GsfTrackReco/interface/GsfTrack.h" #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/Math/interface/LorentzVector.h" #include "FWCore/Framework/interface/Event.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" #include "RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h" #include #include @@ -31,9 +34,10 @@ class LowPtGsfElectronIDProducer final : public edm::global::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions&); private: - double eval(const std::string& name, const edm::Ptr&, double rho, float unbiased) const; + double eval( + const std::string& name, const edm::Ptr&, double rho, float unbiased, float field_z) const; - const edm::EDGetTokenT > gsfElectrons_; + const edm::EDGetTokenT > electrons_; const edm::EDGetTokenT rho_; const edm::EDGetTokenT > unbiased_; const std::vector names_; @@ -42,19 +46,21 @@ class LowPtGsfElectronIDProducer final : public edm::global::EDProducer<> { const double maxPtThreshold_; std::vector > models_; const std::vector thresholds_; + const std::string version_; }; //////////////////////////////////////////////////////////////////////////////// // LowPtGsfElectronIDProducer::LowPtGsfElectronIDProducer(const edm::ParameterSet& conf) - : gsfElectrons_(consumes >(conf.getParameter("electrons"))), + : electrons_(consumes >(conf.getParameter("electrons"))), rho_(consumes(conf.getParameter("rho"))), unbiased_(consumes >(conf.getParameter("unbiased"))), names_(conf.getParameter >("ModelNames")), passThrough_(conf.getParameter("PassThrough")), minPtThreshold_(conf.getParameter("MinPtThreshold")), maxPtThreshold_(conf.getParameter("MaxPtThreshold")), - thresholds_(conf.getParameter >("ModelThresholds")) { + thresholds_(conf.getParameter >("ModelThresholds")), + version_(conf.getParameter("Version")) { for (auto& weights : conf.getParameter >("ModelWeights")) { models_.push_back(createGBRForest(edm::FileInPath(weights))); } @@ -66,6 +72,9 @@ LowPtGsfElectronIDProducer::LowPtGsfElectronIDProducer(const edm::ParameterSet& throw cms::Exception("Incorrect configuration") << "'ModelWeights' size (" << models_.size() << ") != 'ModelThresholds' size (" << thresholds_.size() << ").\n"; } + if (version_ != "V0" && version_ != "V1") { + throw cms::Exception("Incorrect configuration") << "Unknown Version: " << version_ << "\n"; + } for (const auto& name : names_) { produces >(name); } @@ -73,44 +82,41 @@ LowPtGsfElectronIDProducer::LowPtGsfElectronIDProducer(const edm::ParameterSet& //////////////////////////////////////////////////////////////////////////////// // -void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const edm::EventSetup&) const { +void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& setup) const { + // Get z-component of B field + edm::ESHandle field; + setup.get().get(field); + math::XYZVector zfield(field->inTesla(GlobalPoint(0, 0, 0))); + // Pileup edm::Handle rho; event.getByToken(rho_, rho); if (!rho.isValid()) { std::ostringstream os; os << "Problem accessing rho collection for low-pT electrons" << std::endl; - edm::LogError("InvalidHandle") << os.str(); throw cms::Exception("InvalidHandle", os.str()); } // Retrieve GsfElectrons from Event - edm::Handle > gsfElectrons; - event.getByToken(gsfElectrons_, gsfElectrons); - if (!gsfElectrons.isValid()) { + edm::Handle > electrons; + event.getByToken(electrons_, electrons); + if (!electrons.isValid()) { std::ostringstream os; - os << "Problem accessing low-pT gsfElectrons collection" << std::endl; - edm::LogError("InvalidHandle") << os.str(); + os << "Problem accessing low-pT electrons collection" << std::endl; throw cms::Exception("InvalidHandle", os.str()); } // ElectronSeed unbiased BDT edm::Handle > unbiasedH; event.getByToken(unbiased_, unbiasedH); - if (!unbiasedH.isValid()) { - std::ostringstream os; - os << "Problem accessing low-pT 'unbiased' ElectronSeed collection" << std::endl; - edm::LogError("InvalidHandle") << os.str(); - throw cms::Exception("InvalidHandle", os.str()); - } // Iterate through Electrons, evaluate BDT, and store result std::vector > output; for (unsigned int iname = 0; iname < names_.size(); ++iname) { - output.emplace_back(gsfElectrons->size(), -999.); + output.emplace_back(electrons->size(), -999.); } - for (unsigned int iele = 0; iele < gsfElectrons->size(); iele++) { - edm::Ptr ele(gsfElectrons, iele); + for (unsigned int iele = 0; iele < electrons->size(); iele++) { + edm::Ptr ele(electrons, iele); if (ele->core().isNull()) { continue; @@ -123,7 +129,7 @@ void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const //if ( !passThrough_ && ( ele->pt() < minPtThreshold_ ) ) { continue; } for (unsigned int iname = 0; iname < names_.size(); ++iname) { - output[iname][iele] = eval(names_[iname], ele, *rho, unbiased); + output[iname][iele] = eval(names_[iname], ele, *rho, unbiased, zfield.z()); } } @@ -131,7 +137,7 @@ 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(gsfElectrons, output[iname].begin(), output[iname].end()); + filler.insert(electrons, output[iname].begin(), output[iname].end()); filler.fill(); event.put(std::move(ptr), names_[iname]); } @@ -139,14 +145,17 @@ void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const ////////////////////////////////////////////////////////////////////////////////////////// // -double LowPtGsfElectronIDProducer::eval(const std::string& name, - const edm::Ptr& ele, - double rho, - float unbiased) const { +double LowPtGsfElectronIDProducer::eval( + const std::string& name, const edm::Ptr& ele, double rho, float unbiased, float field_z) const { auto iter = std::find(names_.begin(), names_.end(), name); if (iter != names_.end()) { int index = std::distance(names_.begin(), iter); - std::vector inputs = lowptgsfeleid::features(ele, rho, unbiased); + std::vector inputs; + if (version_ == "V0") { + inputs = lowptgsfeleid::features_V0(*ele, rho, unbiased); + } else if (version_ == "V1") { + inputs = lowptgsfeleid::features_V1(*ele, rho, unbiased, field_z); + } return models_.at(index)->GetResponse(inputs.data()); } else { throw cms::Exception("Unknown model name") << "'Name given: '" << name << "'. Check against configuration file.\n"; @@ -160,15 +169,15 @@ void LowPtGsfElectronIDProducer::fillDescriptions(edm::ConfigurationDescriptions edm::ParameterSetDescription desc; desc.add("electrons", edm::InputTag("lowPtGsfElectrons")); desc.add("unbiased", edm::InputTag("lowPtGsfElectronSeedValueMaps:unbiased")); - desc.add("rho", edm::InputTag("fixedGridRhoFastjetAllTmp")); + desc.add("rho", edm::InputTag("fixedGridRhoFastjetAll")); desc.add >("ModelNames", {""}); desc.add >( - "ModelWeights", - {"RecoEgamma/ElectronIdentification/data/LowPtElectrons/RunII_Autumn18_LowPtElectrons_mva_id.xml.gz"}); - desc.add >("ModelThresholds", {-10.}); + "ModelWeights", {"RecoEgamma/ElectronIdentification/data/LowPtElectrons/LowPtElectrons_ID_2020Sept15.root"}); + desc.add >("ModelThresholds", {-99.}); desc.add("PassThrough", false); desc.add("MinPtThreshold", 0.5); desc.add("MaxPtThreshold", 15.); + desc.add("Version", "V1"); descriptions.add("lowPtGsfElectronID", desc); } diff --git a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedValueMapsProducer.cc b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedValueMapsProducer.cc index 94340b6592a48..876b71a6a143b 100644 --- a/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedValueMapsProducer.cc +++ b/RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronSeedValueMapsProducer.cc @@ -1,4 +1,5 @@ #include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" #include "DataFormats/EgammaReco/interface/ElectronSeed.h" #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h" #include "DataFormats/GsfTrackReco/interface/GsfTrack.h" @@ -14,6 +15,7 @@ #include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/stream/EDProducer.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/Utilities/interface/transform.h" #include #include @@ -27,66 +29,112 @@ class LowPtGsfElectronSeedValueMapsProducer : public edm::stream::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions&); private: - const edm::EDGetTokenT gsfTracks_; - const edm::EDGetTokenT > preIdsValueMap_; - const std::vector names_; + edm::EDGetTokenT gsfTracks_; + edm::EDGetTokenT > preIdsValueMap_; + std::vector names_; + const bool rekey_; + edm::EDGetTokenT gsfElectrons_; + std::vector > > floatValueMaps_; }; //////////////////////////////////////////////////////////////////////////////// // LowPtGsfElectronSeedValueMapsProducer::LowPtGsfElectronSeedValueMapsProducer(const edm::ParameterSet& conf) - : gsfTracks_(consumes(conf.getParameter("gsfTracks"))), - preIdsValueMap_(consumes >(conf.getParameter("preIdsValueMap"))), - names_(conf.getParameter >("ModelNames")) { - for (const auto& name : names_) { - produces >(name); + : gsfTracks_(), + preIdsValueMap_(), + names_(), + rekey_(conf.getParameter("rekey")), + gsfElectrons_(), + floatValueMaps_() { + if (rekey_) { + gsfElectrons_ = consumes(conf.getParameter("gsfElectrons")); + std::vector tags = conf.getParameter >("floatValueMaps"); + for (const auto& tag : tags) { + floatValueMaps_ = edm::vector_transform( + tags, [this](edm::InputTag const& tag) { return consumes >(tag); }); + names_.push_back(tag.instance()); + produces >(tag.instance()); + } + } else { + gsfTracks_ = consumes(conf.getParameter("gsfTracks")); + preIdsValueMap_ = consumes >(conf.getParameter("preIdsValueMap")); + names_ = conf.getParameter >("ModelNames"); + for (const auto& name : names_) { + produces >(name); + } } } //////////////////////////////////////////////////////////////////////////////// // void LowPtGsfElectronSeedValueMapsProducer::produce(edm::Event& event, const edm::EventSetup&) { - // Retrieve GsfTracks from Event - edm::Handle gsfTracks; - event.getByToken(gsfTracks_, gsfTracks); - if (!gsfTracks.isValid()) { - edm::LogError("Problem with gsfTracks handle"); - } + if (rekey_ == false) { + // TRANSFORM VALUEMAP OF PREID OBJECTS KEYED BY KF TRACK ... + // .. INTO VALUEMAP OF FLOATS (BDT SCORE) KEYED BY GSF TRACK ... - // Retrieve PreIds from Event - edm::Handle > preIdsValueMap; - event.getByToken(preIdsValueMap_, preIdsValueMap); - if (!preIdsValueMap.isValid()) { - edm::LogError("Problem with preIdsValueMap handle"); - } + // Retrieve GsfTracks from Event + auto gsfTracks = event.getHandle(gsfTracks_); - // Iterate through GsfTracks, extract BDT output, and store result in ValueMap for each model - std::vector > output; - for (unsigned int iname = 0; iname < names_.size(); ++iname) { - output.push_back(std::vector(gsfTracks->size(), -999.)); - } - for (unsigned int igsf = 0; igsf < gsfTracks->size(); igsf++) { - reco::GsfTrackRef gsf(gsfTracks, igsf); - if (gsf.isNonnull() && gsf->extra().isNonnull() && gsf->extra()->seedRef().isNonnull()) { - reco::ElectronSeedRef seed = gsf->extra()->seedRef().castTo(); - if (seed.isNonnull() && seed->ctfTrack().isNonnull()) { - const reco::PreIdRef preid = (*preIdsValueMap)[seed->ctfTrack()]; - if (preid.isNonnull()) { - for (unsigned int iname = 0; iname < names_.size(); ++iname) { - output[iname][igsf] = preid->mva(iname); + // Retrieve PreIds from Event + auto preIdsValueMap = event.getHandle(preIdsValueMap_); + + // Iterate through GsfTracks, extract BDT output, and store result in ValueMap for each model + std::vector > output; + for (unsigned int iname = 0; iname < names_.size(); ++iname) { + output.push_back(std::vector(gsfTracks->size(), -999.)); + } + auto const& gsfTracksV = *gsfTracks; + for (unsigned int igsf = 0; igsf < gsfTracksV.size(); igsf++) { + const reco::GsfTrack& gsf = gsfTracksV[igsf]; + if (gsf.extra().isNonnull() && gsf.extra()->seedRef().isNonnull()) { + reco::ElectronSeedRef seed = gsf.extra()->seedRef().castTo(); + if (seed.isNonnull() && seed->ctfTrack().isNonnull()) { + const reco::PreIdRef preid = (*preIdsValueMap)[seed->ctfTrack()]; + if (preid.isNonnull()) { + for (unsigned int iname = 0; iname < names_.size(); ++iname) { + output[iname][igsf] = preid->mva(iname); + } } } } } - } - // Create and put ValueMap in Event - for (unsigned int iname = 0; iname < names_.size(); ++iname) { - auto ptr = std::make_unique >(edm::ValueMap()); - edm::ValueMap::Filler filler(*ptr); - filler.insert(gsfTracks, output[iname].begin(), output[iname].end()); - filler.fill(); - event.put(std::move(ptr), names_[iname]); + // Create and put ValueMap in Event + for (unsigned int iname = 0; iname < names_.size(); ++iname) { + auto ptr = std::make_unique >(edm::ValueMap()); + edm::ValueMap::Filler filler(*ptr); + filler.insert(gsfTracks, output[iname].begin(), output[iname].end()); + filler.fill(); + event.put(std::move(ptr), names_[iname]); + } + + } else { + // TRANSFORM VALUEMAP OF FLOATS (BDT SCORE) KEYED BY GSF TRACK ... + // .. INTO VALUEMAP OF FLOATS (BDT SCORE) KEYED BY GSF ELECTRON ... + + // Retrieve GsfElectrons from Event + auto gsfElectrons = event.getHandle(gsfElectrons_); + + // Retrieve float ValueMaps from Event + for (unsigned int idx = 0; idx < names_.size(); ++idx) { + // Extract ValueMap from Event + auto const& floatValueMap = event.get(floatValueMaps_[idx]); + + // Store BDT scores in vector + std::vector output(gsfElectrons->size(), -99.); + auto const& gsfElectronsV = *gsfElectrons; + for (unsigned int iele = 0; iele < gsfElectronsV.size(); iele++) { + const reco::GsfElectron& ele = gsfElectronsV[iele]; + reco::GsfTrackRef gsf = ele.gsfTrack(); + output[iele] = floatValueMap[gsf]; + } + // Create and put ValueMap in Event + auto ptr = std::make_unique >(edm::ValueMap()); + edm::ValueMap::Filler filler(*ptr); + filler.insert(gsfElectrons, output.begin(), output.end()); + filler.fill(); + event.put(std::move(ptr), names_[idx]); + } } } @@ -97,6 +145,9 @@ void LowPtGsfElectronSeedValueMapsProducer::fillDescriptions(edm::ConfigurationD desc.add("gsfTracks", edm::InputTag("lowPtGsfEleGsfTracks")); desc.add("preIdsValueMap", edm::InputTag("lowPtGsfElectronSeeds")); desc.add >("ModelNames", {"unbiased", "ptbiased"}); + desc.add("rekey", false); + desc.add("gsfElectrons", edm::InputTag()); + desc.add >("floatValueMaps", std::vector()); descriptions.add("lowPtGsfElectronSeedValueMaps", desc); } diff --git a/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectronSeedValueMaps_cff.py b/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectronSeedValueMaps_cff.py new file mode 100644 index 0000000000000..4ee7a8a749796 --- /dev/null +++ b/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectronSeedValueMaps_cff.py @@ -0,0 +1,11 @@ +import FWCore.ParameterSet.Config as cms +# Low pT Electron value maps +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronSeedValueMaps_cfi import lowPtGsfElectronSeedValueMaps + +# Low pT Electron value maps, rekeyed by reco::GsfElectron +rekeyLowPtGsfElectronSeedValueMaps = lowPtGsfElectronSeedValueMaps.clone( + rekey=True, + gsfElectrons="lowPtGsfElectrons", + floatValueMaps=["lowPtGsfElectronSeedValueMaps:unbiased", + "lowPtGsfElectronSeedValueMaps:ptbiased"], +) diff --git a/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectronSequence_cff.py b/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectronSequence_cff.py index f8c021ce25f8d..270f26685a3a3 100644 --- a/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectronSequence_cff.py +++ b/RecoEgamma/EgammaElectronProducers/python/lowPtGsfElectronSequence_cff.py @@ -69,12 +69,12 @@ from RecoEgamma.EgammaElectronProducers.lowPtGsfElectrons_cfi import * # Low pT Electron value maps -from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronSeedValueMaps_cfi import lowPtGsfElectronSeedValueMaps +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronSeedValueMaps_cff import lowPtGsfElectronSeedValueMaps +from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronSeedValueMaps_cff import rekeyLowPtGsfElectronSeedValueMaps # Low pT Electron ID from RecoEgamma.EgammaElectronProducers.lowPtGsfElectronID_cfi import lowPtGsfElectronID - # Full sequence lowPtGsfElectronTask = cms.Task(lowPtGsfElePfTracks, lowPtGsfElectronSeeds, @@ -86,6 +86,7 @@ lowPtGsfElectronCores, lowPtGsfElectrons, lowPtGsfElectronSeedValueMaps, + rekeyLowPtGsfElectronSeedValueMaps, lowPtGsfElectronID ) lowPtGsfElectronSequence = cms.Sequence(lowPtGsfElectronTask) diff --git a/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc b/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc index a346367b7900f..b09f7037b16cf 100644 --- a/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc +++ b/RecoEgamma/EgammaElectronProducers/src/LowPtGsfElectronFeatures.cc @@ -6,6 +6,8 @@ #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "TVector3.h" +#include +#include namespace lowptgsfeleseed { @@ -129,7 +131,7 @@ namespace lowptgsfeleseed { namespace lowptgsfeleid { - std::vector features(edm::Ptr const& ele, float rho, float unbiased) { + std::vector features_V1(reco::GsfElectron const& ele, float rho, float unbiased, float field_z) { float eid_rho = -999.; float eid_sc_eta = -999.; float eid_shape_full5x5_r9 = -999.; @@ -165,8 +167,8 @@ namespace lowptgsfeleid { float sc_clus2_E_ov_p = -999.; // KF tracks - if (ele->core().isNonnull()) { - const auto& trk = ele->closestCtfTrackRef(); // reco::TrackRef + if (ele.core().isNonnull()) { + reco::TrackRef trk = ele.closestCtfTrackRef(); if (trk.isNonnull()) { eid_trk_p = (float)trk->p(); eid_trk_nhits = (float)trk->found(); @@ -174,14 +176,14 @@ namespace lowptgsfeleid { TVector3 trkTV3(0, 0, 0); trkTV3.SetPtEtaPhi(trk->pt(), trk->eta(), trk->phi()); TVector3 eleTV3(0, 0, 0); - eleTV3.SetPtEtaPhi(ele->pt(), ele->eta(), ele->phi()); + eleTV3.SetPtEtaPhi(ele.pt(), ele.eta(), ele.phi()); trk_dr = eleTV3.DeltaR(trkTV3); } } // GSF tracks - if (ele->core().isNonnull()) { - const auto& gsf = ele->core()->gsfTrack(); // reco::GsfTrackRef + if (ele.core().isNonnull()) { + reco::GsfTrackRef gsf = ele.core()->gsfTrack(); if (gsf.isNonnull()) { gsf_mode_p = gsf->pMode(); eid_gsf_nhits = (float)gsf->found(); @@ -189,59 +191,282 @@ namespace lowptgsfeleid { TVector3 gsfTV3(0, 0, 0); gsfTV3.SetPtEtaPhi(gsf->ptMode(), gsf->etaMode(), gsf->phiMode()); TVector3 eleTV3(0, 0, 0); - eleTV3.SetPtEtaPhi(ele->pt(), ele->eta(), ele->phi()); + eleTV3.SetPtEtaPhi(ele.pt(), ele.eta(), ele.phi()); gsf_dr = eleTV3.DeltaR(gsfTV3); } } // Super clusters - if (ele->core().isNonnull()) { - const auto& sc = ele->core()->superCluster(); // reco::SuperClusterRef + if (ele.core().isNonnull()) { + reco::SuperClusterRef sc = ele.core()->superCluster(); if (sc.isNonnull()) { eid_sc_E = sc->energy(); eid_sc_eta = sc->eta(); eid_sc_etaWidth = sc->etaWidth(); eid_sc_phiWidth = sc->phiWidth(); - sc_Nclus = (float)sc->clustersSize(); + sc_Nclus = sc->clustersSize(); } } // Track-cluster matching - if (ele.isNonnull()) { - eid_match_seed_dEta = ele->deltaEtaSeedClusterTrackAtCalo(); - eid_match_eclu_EoverP = (1. / ele->ecalEnergy()) - (1. / ele->p()); - eid_match_SC_EoverP = ele->eSuperClusterOverP(); - eid_match_SC_dEta = ele->deltaEtaSuperClusterTrackAtVtx(); - eid_match_SC_dPhi = ele->deltaPhiSuperClusterTrackAtVtx(); - } + eid_match_seed_dEta = ele.deltaEtaSeedClusterTrackAtCalo(); + eid_match_eclu_EoverP = (1. / ele.ecalEnergy()) - (1. / ele.p()); + eid_match_SC_EoverP = ele.eSuperClusterOverP(); + eid_match_SC_dEta = ele.deltaEtaSuperClusterTrackAtVtx(); + eid_match_SC_dPhi = ele.deltaPhiSuperClusterTrackAtVtx(); // Shower shape vars - if (ele.isNonnull()) { - eid_shape_full5x5_HoverE = ele->full5x5_hcalOverEcal(); - eid_shape_full5x5_r9 = ele->full5x5_r9(); - } + eid_shape_full5x5_HoverE = ele.full5x5_hcalOverEcal(); + eid_shape_full5x5_r9 = ele.full5x5_r9(); // Misc eid_rho = rho; - if (ele.isNonnull()) { - eid_brem_frac = ele->fbrem(); - core_shFracHits = (float)ele->shFracInnerHits(); + eid_brem_frac = ele.fbrem(); + core_shFracHits = ele.shFracInnerHits(); + + // Unbiased BDT from ElectronSeed + gsf_bdtout1 = unbiased; + + // Clusters + if (ele.core().isNonnull()) { + reco::GsfTrackRef gsf = ele.core()->gsfTrack(); + if (gsf.isNonnull()) { + reco::SuperClusterRef sc = ele.core()->superCluster(); + if (sc.isNonnull()) { + // Propagate electron track to ECAL surface + double mass2 = 0.000511 * 0.000511; + float p2 = pow(gsf->p(), 2); + float energy = sqrt(mass2 + p2); + XYZTLorentzVector mom = XYZTLorentzVector(gsf->px(), gsf->py(), gsf->pz(), energy); + XYZTLorentzVector pos = XYZTLorentzVector(gsf->vx(), gsf->vy(), gsf->vz(), 0.); + BaseParticlePropagator propagator(RawParticle(mom, pos, gsf->charge()), 0, 0, field_z); + propagator.propagateToEcalEntrance(true); // true only first half loop , false more than one loop + bool reach_ECAL = propagator.getSuccess(); // 0 does not reach ECAL, 1 yes barrel, 2 yes endcaps + // ECAL entry point for track + GlobalPoint ecal_pos(propagator.particle().x(), propagator.particle().y(), propagator.particle().z()); + + // Track-cluster matching for most energetic clusters + sc_clus1_nxtal = -999; + sc_clus1_dphi = -999.; + sc_clus2_dphi = -999.; + sc_clus1_deta = -999.; + sc_clus2_deta = -999.; + sc_clus1_E = -999.; + sc_clus2_E = -999.; + sc_clus1_E_ov_p = -999.; + sc_clus2_E_ov_p = -999.; + trackClusterMatching(*sc, + *gsf, + reach_ECAL, + ecal_pos, + sc_clus1_nxtal, + sc_clus1_dphi, + sc_clus2_dphi, + sc_clus1_deta, + sc_clus2_deta, + sc_clus1_E, + sc_clus2_E, + sc_clus1_E_ov_p, + sc_clus2_E_ov_p); + sc_clus1_nxtal = (int)sc_clus1_nxtal; + + } // sc.isNonnull() + } // gsf.isNonnull() + } // clusters + + // Out-of-range + eid_rho = std::clamp(eid_rho, 0.f, 100.f); + eid_sc_eta = std::clamp(eid_sc_eta, -5.f, 5.f); + eid_shape_full5x5_r9 = std::clamp(eid_shape_full5x5_r9, 0.f, 2.f); + eid_sc_etaWidth = std::clamp(eid_sc_etaWidth, 0.f, 3.14f); + eid_sc_phiWidth = std::clamp(eid_sc_phiWidth, 0.f, 3.14f); + eid_shape_full5x5_HoverE = std::clamp(eid_shape_full5x5_HoverE, 0.f, 50.f); + eid_trk_nhits = std::clamp(eid_trk_nhits, -1.f, 50.f); + eid_trk_chi2red = std::clamp(eid_trk_chi2red, -1.f, 50.f); + eid_gsf_chi2red = std::clamp(eid_gsf_chi2red, -1.f, 100.f); + if (eid_brem_frac < 0.) + eid_brem_frac = -1.; // + if (eid_brem_frac > 1.) + eid_brem_frac = 1.; // + eid_gsf_nhits = std::clamp(eid_gsf_nhits, -1.f, 50.f); + eid_match_SC_EoverP = std::clamp(eid_match_SC_EoverP, 0.f, 100.f); + eid_match_eclu_EoverP = std::clamp(eid_match_eclu_EoverP, -1.f, 1.f); + eid_match_SC_dEta = std::clamp(eid_match_SC_dEta, -10.f, 10.f); + eid_match_SC_dPhi = std::clamp(eid_match_SC_dPhi, -3.14f, 3.14f); + eid_match_seed_dEta = std::clamp(eid_match_seed_dEta, -10.f, 10.f); + eid_sc_E = std::clamp(eid_sc_E, 0.f, 1000.f); + eid_trk_p = std::clamp(eid_trk_p, -1.f, 1000.f); + gsf_mode_p = std::clamp(gsf_mode_p, 0.f, 1000.f); + core_shFracHits = std::clamp(core_shFracHits, 0.f, 1.f); + gsf_bdtout1 = std::clamp(gsf_bdtout1, -20.f, 20.f); + if (gsf_dr < 0.) + gsf_dr = 5.; // + if (gsf_dr > 5.) + gsf_dr = 5.; // + if (trk_dr < 0.) + trk_dr = 5.; // + if (trk_dr > 5.) + trk_dr = 5.; // + sc_Nclus = std::clamp(sc_Nclus, 0.f, 20.f); + sc_clus1_nxtal = std::clamp(sc_clus1_nxtal, 0.f, 100.f); + sc_clus1_dphi = std::clamp(sc_clus1_dphi, -3.14f, 3.14f); + sc_clus2_dphi = std::clamp(sc_clus2_dphi, -3.14f, 3.14f); + sc_clus1_deta = std::clamp(sc_clus1_deta, -5.f, 5.f); + sc_clus2_deta = std::clamp(sc_clus2_deta, -5.f, 5.f); + sc_clus1_E = std::clamp(sc_clus1_E, 0.f, 1000.f); + sc_clus2_E = std::clamp(sc_clus2_E, 0.f, 1000.f); + if (sc_clus1_E_ov_p < 0.) + sc_clus1_E_ov_p = -1.; // + if (sc_clus2_E_ov_p < 0.) + sc_clus2_E_ov_p = -1.; // + + // Set contents of vector + std::vector output = {eid_rho, + eid_sc_eta, + eid_shape_full5x5_r9, + eid_sc_etaWidth, + eid_sc_phiWidth, + eid_shape_full5x5_HoverE, + eid_trk_nhits, + eid_trk_chi2red, + eid_gsf_chi2red, + eid_brem_frac, + eid_gsf_nhits, + eid_match_SC_EoverP, + eid_match_eclu_EoverP, + eid_match_SC_dEta, + eid_match_SC_dPhi, + eid_match_seed_dEta, + eid_sc_E, + eid_trk_p, + gsf_mode_p, + core_shFracHits, + gsf_bdtout1, + gsf_dr, + trk_dr, + sc_Nclus, + sc_clus1_nxtal, + sc_clus1_dphi, + sc_clus2_dphi, + sc_clus1_deta, + sc_clus2_deta, + sc_clus1_E, + sc_clus2_E, + sc_clus1_E_ov_p, + sc_clus2_E_ov_p}; + return output; + } + + //////////////////////////////////////////////////////////////////////////////// + // feature list for original models (2019Aug07 and earlier) + std::vector features_V0(reco::GsfElectron const& ele, float rho, float unbiased) { + float eid_rho = -999.; + float eid_sc_eta = -999.; + float eid_shape_full5x5_r9 = -999.; + float eid_sc_etaWidth = -999.; + float eid_sc_phiWidth = -999.; + float eid_shape_full5x5_HoverE = -999.; + float eid_trk_nhits = -999.; + float eid_trk_chi2red = -999.; + float eid_gsf_chi2red = -999.; + float eid_brem_frac = -999.; + float eid_gsf_nhits = -999.; + float eid_match_SC_EoverP = -999.; + float eid_match_eclu_EoverP = -999.; + float eid_match_SC_dEta = -999.; + float eid_match_SC_dPhi = -999.; + float eid_match_seed_dEta = -999.; + float eid_sc_E = -999.; + float eid_trk_p = -999.; + float gsf_mode_p = -999.; + float core_shFracHits = -999.; + float gsf_bdtout1 = -999.; + float gsf_dr = -999.; + float trk_dr = -999.; + float sc_Nclus = -999.; + float sc_clus1_nxtal = -999.; + float sc_clus1_dphi = -999.; + float sc_clus2_dphi = -999.; + float sc_clus1_deta = -999.; + float sc_clus2_deta = -999.; + float sc_clus1_E = -999.; + float sc_clus2_E = -999.; + float sc_clus1_E_ov_p = -999.; + float sc_clus2_E_ov_p = -999.; + + // KF tracks + if (ele.core().isNonnull()) { + const auto& trk = ele.closestCtfTrackRef(); // reco::TrackRef + if (trk.isNonnull()) { + eid_trk_p = (float)trk->p(); + eid_trk_nhits = (float)trk->found(); + eid_trk_chi2red = (float)trk->normalizedChi2(); + TVector3 trkTV3(0, 0, 0); + trkTV3.SetPtEtaPhi(trk->pt(), trk->eta(), trk->phi()); + TVector3 eleTV3(0, 0, 0); + eleTV3.SetPtEtaPhi(ele.pt(), ele.eta(), ele.phi()); + trk_dr = eleTV3.DeltaR(trkTV3); + } + } + + // GSF tracks + if (ele.core().isNonnull()) { + const auto& gsf = ele.core()->gsfTrack(); // reco::GsfTrackRef + if (gsf.isNonnull()) { + gsf_mode_p = gsf->pMode(); + eid_gsf_nhits = (float)gsf->found(); + eid_gsf_chi2red = gsf->normalizedChi2(); + TVector3 gsfTV3(0, 0, 0); + gsfTV3.SetPtEtaPhi(gsf->ptMode(), gsf->etaMode(), gsf->phiMode()); + TVector3 eleTV3(0, 0, 0); + eleTV3.SetPtEtaPhi(ele.pt(), ele.eta(), ele.phi()); + gsf_dr = eleTV3.DeltaR(gsfTV3); + } + } + + // Super clusters + if (ele.core().isNonnull()) { + const auto& sc = ele.core()->superCluster(); // reco::SuperClusterRef + if (sc.isNonnull()) { + eid_sc_E = sc->energy(); + eid_sc_eta = sc->eta(); + eid_sc_etaWidth = sc->etaWidth(); + eid_sc_phiWidth = sc->phiWidth(); + sc_Nclus = (float)sc->clustersSize(); + } } + // Track-cluster matching + eid_match_seed_dEta = ele.deltaEtaSeedClusterTrackAtCalo(); + eid_match_eclu_EoverP = (1. / ele.ecalEnergy()) - (1. / ele.p()); + eid_match_SC_EoverP = ele.eSuperClusterOverP(); + eid_match_SC_dEta = ele.deltaEtaSuperClusterTrackAtVtx(); + eid_match_SC_dPhi = ele.deltaPhiSuperClusterTrackAtVtx(); + + // Shower shape vars + eid_shape_full5x5_HoverE = ele.full5x5_hcalOverEcal(); + eid_shape_full5x5_r9 = ele.full5x5_r9(); + + // Misc + eid_rho = rho; + + eid_brem_frac = ele.fbrem(); + core_shFracHits = (float)ele.shFracInnerHits(); + // Unbiased BDT from ElectronSeed gsf_bdtout1 = unbiased; // Clusters - if (ele->core().isNonnull()) { - const auto& gsf = ele->core()->gsfTrack(); // reco::GsfTrackRef + if (ele.core().isNonnull()) { + const auto& gsf = ele.core()->gsfTrack(); // reco::GsfTrackRef if (gsf.isNonnull()) { - const auto& sc = ele->core()->superCluster(); // reco::SuperClusterRef + const auto& sc = ele.core()->superCluster(); // reco::SuperClusterRef if (sc.isNonnull()) { // Propagate electron track to ECAL surface - double mass_ = 0.000511 * 0.000511; + double mass2 = 0.000511 * 0.000511; float p2 = pow(gsf->p(), 2); - float energy = sqrt(mass_ + p2); + float energy = sqrt(mass2 + p2); math::XYZTLorentzVector mom = math::XYZTLorentzVector(gsf->px(), gsf->py(), gsf->pz(), energy); math::XYZTLorentzVector pos = math::XYZTLorentzVector(gsf->vx(), gsf->vy(), gsf->vz(), 0.); float field_z = 3.8; @@ -253,41 +478,7 @@ namespace lowptgsfeleid { GlobalPoint ecal_pos( mypart.particle().vertex().x(), mypart.particle().vertex().y(), mypart.particle().vertex().z()); - // Iterate through ECAL clusters and sort in energy - int clusNum = 0; - float maxEne1 = -1; - float maxEne2 = -1; - int i1 = -1; - int i2 = -1; - try { - if (sc->clustersSize() > 0 && sc->clustersBegin() != sc->clustersEnd()) { - for (const auto& cluster : sc->clusters()) { - if (cluster->energy() > maxEne1) { - maxEne1 = cluster->energy(); - i1 = clusNum; - } - clusNum++; - } - if (sc->clustersSize() > 1) { - clusNum = 0; - for (const auto& cluster : sc->clusters()) { - if (clusNum != i1) { - if (cluster->energy() > maxEne2) { - maxEne2 = cluster->energy(); - i2 = clusNum; - } - } - clusNum++; - } - } - } // loop over clusters - } catch (...) { - edm::LogError("SuperClusters") << "Problem accessing SC constituent clusters:" - << " clusNum=" << clusNum << " clustersSize=" << sc->clustersSize() - << " energy=" << sc->energy() << std::endl; - } - - // Initializations + // Track-cluster matching for most energetic clusters sc_clus1_nxtal = -999.; sc_clus1_dphi = -999.; sc_clus2_dphi = -999.; @@ -297,178 +488,74 @@ namespace lowptgsfeleid { sc_clus2_E = -999.; sc_clus1_E_ov_p = -999.; sc_clus2_E_ov_p = -999.; + trackClusterMatching(*sc, + *gsf, + reach_ECAL, + ecal_pos, + sc_clus1_nxtal, + sc_clus1_dphi, + sc_clus2_dphi, + sc_clus1_deta, + sc_clus2_deta, + sc_clus1_E, + sc_clus2_E, + sc_clus1_E_ov_p, + sc_clus2_E_ov_p); - // track-clusters match - clusNum = 0; - try { - if (sc->clustersSize() > 0 && sc->clustersBegin() != sc->clustersEnd()) { - for (const auto& cluster : sc->clusters()) { - float deta = std::fabs(ecal_pos.eta() - cluster->eta()); - float dphi = std::fabs(ecal_pos.phi() - cluster->phi()); - if (dphi > M_PI) - dphi -= 2 * M_PI; - if (ecal_pos.phi() - cluster->phi() < 0) - dphi = -dphi; - if (ecal_pos.eta() - cluster->eta() < 0) - deta = -deta; - - if (clusNum == i1) { - sc_clus1_E = cluster->energy(); - if (gsf->pMode() > 0) - sc_clus1_E_ov_p = cluster->energy() / gsf->pMode(); - sc_clus1_nxtal = (float)cluster->size(); - if (reach_ECAL > 0) { - sc_clus1_deta = deta; - sc_clus1_dphi = dphi; - } - } else if (clusNum == i2) { - sc_clus2_E = cluster->energy(); - if (gsf->pMode() > 0) - sc_clus2_E_ov_p = cluster->energy() / gsf->pMode(); - if (reach_ECAL > 0) { - sc_clus2_deta = deta; - sc_clus2_dphi = dphi; - } - } - clusNum++; - } - } - } catch (...) { - edm::LogError("SuperClusters") << "Problem with track-cluster matching" << std::endl; - } - } - } - } // clusters + } // sc.isNonnull() + } // gsf.isNonnull() + } // clusters // Out-of-range - if (eid_rho < 0) - eid_rho = 0; - if (eid_rho > 100) - eid_rho = 100; - if (eid_sc_eta < -5) - eid_sc_eta = -5; - if (eid_sc_eta > 5) - eid_sc_eta = 5; - if (eid_shape_full5x5_r9 < 0) - eid_shape_full5x5_r9 = 0; - if (eid_shape_full5x5_r9 > 2) - eid_shape_full5x5_r9 = 2; - if (eid_sc_etaWidth < 0) - eid_sc_etaWidth = 0; - if (eid_sc_etaWidth > 3.14) - eid_sc_etaWidth = 3.14; - if (eid_sc_phiWidth < 0) - eid_sc_phiWidth = 0; - if (eid_sc_phiWidth > 3.14) - eid_sc_phiWidth = 3.14; - if (eid_shape_full5x5_HoverE < 0) - eid_shape_full5x5_HoverE = 0; - if (eid_shape_full5x5_HoverE > 50) - eid_shape_full5x5_HoverE = 50; - if (eid_trk_nhits < -1) - eid_trk_nhits = -1; - if (eid_trk_nhits > 50) - eid_trk_nhits = 50; - if (eid_trk_chi2red < -1) - eid_trk_chi2red = -1; - if (eid_trk_chi2red > 50) - eid_trk_chi2red = 50; - if (eid_gsf_chi2red < -1) - eid_gsf_chi2red = -1; - if (eid_gsf_chi2red > 100) - eid_gsf_chi2red = 100; - if (eid_brem_frac < 0) - eid_brem_frac = -1; - if (eid_brem_frac > 1) - eid_brem_frac = 1; - if (eid_gsf_nhits < -1) - eid_gsf_nhits = -1; - if (eid_gsf_nhits > 50) - eid_gsf_nhits = 50; - if (eid_match_SC_EoverP < 0) - eid_match_SC_EoverP = 0; - if (eid_match_SC_EoverP > 100) - eid_match_SC_EoverP = 100; - if (eid_match_eclu_EoverP < -1.) - eid_match_eclu_EoverP = -1.; - if (eid_match_eclu_EoverP > 1.) - eid_match_eclu_EoverP = 1.; - if (eid_match_SC_dEta < -10) - eid_match_SC_dEta = -10; - if (eid_match_SC_dEta > 10) - eid_match_SC_dEta = 10; - if (eid_match_SC_dPhi < -3.14) - eid_match_SC_dPhi = -3.14; - if (eid_match_SC_dPhi > 3.14) - eid_match_SC_dPhi = 3.14; - if (eid_match_seed_dEta < -10) - eid_match_seed_dEta = -10; - if (eid_match_seed_dEta > 10) - eid_match_seed_dEta = 10; - if (eid_sc_E < 0) - eid_sc_E = 0; - if (eid_sc_E > 1000) - eid_sc_E = 1000; - if (eid_trk_p < -1) - eid_trk_p = -1; - if (eid_trk_p > 1000) - eid_trk_p = 1000; - if (gsf_mode_p < 0) - gsf_mode_p = 0; - if (gsf_mode_p > 1000) - gsf_mode_p = 1000; - if (core_shFracHits < 0) - core_shFracHits = 0; - if (core_shFracHits > 1) - core_shFracHits = 1; - if (gsf_bdtout1 < -20) - gsf_bdtout1 = -20; - if (gsf_bdtout1 > 20) - gsf_bdtout1 = 20; - if (gsf_dr < 0) - gsf_dr = 5; - if (gsf_dr > 5) - gsf_dr = 5; - if (trk_dr < 0) - trk_dr = 5; - if (trk_dr > 5) - trk_dr = 5; - if (sc_Nclus < 0) - sc_Nclus = 0; - if (sc_Nclus > 20) - sc_Nclus = 20; - if (sc_clus1_nxtal < 0) - sc_clus1_nxtal = 0; - if (sc_clus1_nxtal > 100) - sc_clus1_nxtal = 100; + eid_sc_eta = std::clamp(eid_sc_eta, -5.f, 5.f); + eid_shape_full5x5_r9 = std::clamp(eid_shape_full5x5_r9, 0.f, 2.f); + eid_sc_etaWidth = std::clamp(eid_sc_etaWidth, 0.f, 3.14f); + eid_sc_phiWidth = std::clamp(eid_sc_phiWidth, 0.f, 3.14f); + eid_shape_full5x5_HoverE = std::clamp(eid_shape_full5x5_HoverE, 0.f, 50.f); + eid_trk_nhits = std::clamp(eid_trk_nhits, -1.f, 50.f); + eid_trk_chi2red = std::clamp(eid_trk_chi2red, -1.f, 50.f); + eid_gsf_chi2red = std::clamp(eid_gsf_chi2red, -1.f, 100.f); + if (eid_brem_frac < 0.) + eid_brem_frac = -1.; // + if (eid_brem_frac > 1.) + eid_brem_frac = 1.; // + eid_gsf_nhits = std::clamp(eid_gsf_nhits, -1.f, 50.f); + eid_match_SC_EoverP = std::clamp(eid_match_SC_EoverP, 0.f, 100.f); + eid_match_eclu_EoverP = std::clamp(eid_match_eclu_EoverP, -1.f, 1.f); + eid_match_SC_dEta = std::clamp(eid_match_SC_dEta, -10.f, 10.f); + eid_match_SC_dPhi = std::clamp(eid_match_SC_dPhi, -3.14f, 3.14f); + eid_match_seed_dEta = std::clamp(eid_match_seed_dEta, -10.f, 10.f); + eid_sc_E = std::clamp(eid_sc_E, 0.f, 1000.f); + eid_trk_p = std::clamp(eid_trk_p, -1.f, 1000.f); + gsf_mode_p = std::clamp(gsf_mode_p, 0.f, 1000.f); + core_shFracHits = std::clamp(core_shFracHits, 0.f, 1.f); + gsf_bdtout1 = std::clamp(gsf_bdtout1, -20.f, 20.f); + if (gsf_dr < 0.) + gsf_dr = 5.; // + if (gsf_dr > 5.) + gsf_dr = 5.; // + if (trk_dr < 0.) + trk_dr = 5.; // + if (trk_dr > 5.) + trk_dr = 5.; // + sc_Nclus = std::clamp(sc_Nclus, 0.f, 20.f); + sc_clus1_nxtal = std::clamp(sc_clus1_nxtal, 0.f, 100.f); if (sc_clus1_dphi < -3.14) - sc_clus1_dphi = -5; + sc_clus1_dphi = -5.; // if (sc_clus1_dphi > 3.14) - sc_clus1_dphi = 5; + sc_clus1_dphi = 5.; // if (sc_clus2_dphi < -3.14) - sc_clus2_dphi = -5; + sc_clus2_dphi = -5.; // if (sc_clus2_dphi > 3.14) - sc_clus2_dphi = 5; - if (sc_clus1_deta < -5) - sc_clus1_deta = -5; - if (sc_clus1_deta > 5) - sc_clus1_deta = 5; - if (sc_clus2_deta < -5) - sc_clus2_deta = -5; - if (sc_clus2_deta > 5) - sc_clus2_deta = 5; - if (sc_clus1_E < 0) - sc_clus1_E = 0; - if (sc_clus1_E > 1000) - sc_clus1_E = 1000; - if (sc_clus2_E < 0) - sc_clus2_E = 0; - if (sc_clus2_E > 1000) - sc_clus2_E = 1000; - if (sc_clus1_E_ov_p < 0) - sc_clus1_E_ov_p = -1; - if (sc_clus2_E_ov_p < 0) - sc_clus2_E_ov_p = -1; + sc_clus2_dphi = 5.; // + sc_clus1_deta = std::clamp(sc_clus1_deta, -5.f, 5.f); + sc_clus2_deta = std::clamp(sc_clus2_deta, -5.f, 5.f); + sc_clus1_E = std::clamp(sc_clus1_E, 0.f, 1000.f); + sc_clus2_E = std::clamp(sc_clus2_E, 0.f, 1000.f); + if (sc_clus1_E_ov_p < 0.) + sc_clus1_E_ov_p = -1.; // + if (sc_clus2_E_ov_p < 0.) + sc_clus2_E_ov_p = -1.; // // Set contents of vector std::vector output = {eid_rho, @@ -508,15 +595,82 @@ namespace lowptgsfeleid { } //////////////////////////////////////////////////////////////////////////////// - // - std::vector features(edm::Ref > const& ele, float rho, float unbiased) { - return features(edm::refToPtr(ele), rho, unbiased); + // Find most energetic clusters + void findEnergeticClusters( + reco::SuperCluster const& sc, int& clusNum, float& maxEne1, float& maxEne2, int& i1, int& i2) { + if (sc.clustersSize() > 0 && sc.clustersBegin() != sc.clustersEnd()) { + for (auto const& cluster : sc.clusters()) { + if (cluster->energy() > maxEne1) { + maxEne1 = cluster->energy(); + i1 = clusNum; + } + clusNum++; + } + if (sc.clustersSize() > 1) { + clusNum = 0; + for (auto const& cluster : sc.clusters()) { + if (clusNum != i1) { + if (cluster->energy() > maxEne2) { + maxEne2 = cluster->energy(); + i2 = clusNum; + } + } + clusNum++; + } + } + } // loop over clusters } //////////////////////////////////////////////////////////////////////////////// - // - std::vector features(edm::Ref > const& ele, float rho, float unbiased) { - return features(edm::refToPtr(ele), rho, unbiased); + // Track-cluster matching for most energetic clusters + void trackClusterMatching(reco::SuperCluster const& sc, + reco::GsfTrack const& gsf, + bool const& reach_ECAL, + GlobalPoint const& ecal_pos, + float& sc_clus1_nxtal, + float& sc_clus1_dphi, + float& sc_clus2_dphi, + float& sc_clus1_deta, + float& sc_clus2_deta, + float& sc_clus1_E, + float& sc_clus2_E, + float& sc_clus1_E_ov_p, + float& sc_clus2_E_ov_p) { + // Iterate through ECAL clusters and sort in energy + int clusNum = 0; + float maxEne1 = -1; + float maxEne2 = -1; + int i1 = -1; + int i2 = -1; + findEnergeticClusters(sc, clusNum, maxEne1, maxEne2, i1, i2); + + // track-clusters match + clusNum = 0; + if (sc.clustersSize() > 0 && sc.clustersBegin() != sc.clustersEnd()) { + for (auto const& cluster : sc.clusters()) { + float deta = ecal_pos.eta() - cluster->eta(); + float dphi = reco::deltaPhi(ecal_pos.phi(), cluster->phi()); + if (clusNum == i1) { + sc_clus1_E = cluster->energy(); + if (gsf.pMode() > 0) + sc_clus1_E_ov_p = cluster->energy() / gsf.pMode(); + sc_clus1_nxtal = (float)cluster->size(); + if (reach_ECAL > 0) { + sc_clus1_deta = deta; + sc_clus1_dphi = dphi; + } + } else if (clusNum == i2) { + sc_clus2_E = cluster->energy(); + if (gsf.pMode() > 0) + sc_clus2_E_ov_p = cluster->energy() / gsf.pMode(); + if (reach_ECAL > 0) { + sc_clus2_deta = deta; + sc_clus2_dphi = dphi; + } + } + clusNum++; + } + } } } // namespace lowptgsfeleid diff --git a/RecoEgamma/EgammaTools/plugins/EGExtraInfoModifierFromValueMaps.cc b/RecoEgamma/EgammaTools/plugins/EGExtraInfoModifierFromValueMaps.cc index 3ae841c81750e..c0ea0bce9953b 100644 --- a/RecoEgamma/EgammaTools/plugins/EGExtraInfoModifierFromValueMaps.cc +++ b/RecoEgamma/EgammaTools/plugins/EGExtraInfoModifierFromValueMaps.cc @@ -1,6 +1,7 @@ #include "RecoEgamma/EgammaTools/interface/EGExtraInfoModifierFromValueMaps.h" #include "DataFormats/EgammaCandidates/interface/HIPhotonIsolation.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" using EGExtraInfoModifierFromFloatValueMaps = EGExtraInfoModifierFromValueMaps; DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, @@ -45,3 +46,8 @@ using EGExtraInfoModifierFromEGIDValueMaps = EGExtraInfoModifierFromValueMaps >; +DEFINE_EDM_PLUGIN(ModifyObjectValueFactory, + EGExtraInfoModifierFromPackedCandPtrValueMaps, + "EGExtraInfoModifierFromPackedCandPtrValueMaps");