From d84f4dea2fb5a70746b2ec39f3140c97d6c8c325 Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Thu, 10 Oct 2024 19:23:16 +0200 Subject: [PATCH 01/10] add unbiased SCs pt > 5 GeV to miniAOD --- .../PatAlgos/python/slimming/slimming_cff.py | 1 + .../python/reducedEgamma_cfi.py | 11 ++++++++ .../src/ReducedEGProducer.cc | 27 ++++++++++++++++++- 3 files changed, 38 insertions(+), 1 deletion(-) diff --git a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py index 17aeba0c5b5fe..92499ab0b21d1 100644 --- a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py @@ -66,6 +66,7 @@ slimmedLambdaVertices, slimmedMETs, metFilterPathsTask, + superClusterMerger, reducedEgamma, slimmedHcalRecHits, bunchSpacingProducer, diff --git a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py index 6001c246d1e38..458f0f6868a7b 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py @@ -4,6 +4,9 @@ from RecoJets.Configuration.CaloTowersES_cfi import * reducedEgamma = cms.EDProducer("ReducedEGProducer", + keepPfSuperclusterPtMin = cms.double(5.), + keepPfSuperclusterAbsetaMax = cms.double(2.5), + relinkSupercluster = cms.bool(True), keepPhotons = cms.string("hadTowOverEm()<0.15 && pt>10 && (pt>14 || chargedHadronIso()<10)"), #keep in output slimRelinkPhotons = cms.string("hadTowOverEm()<0.15 && pt>10 && (pt>14 || chargedHadronIso()<10)"), #keep only slimmed SuperCluster plus seed cluster relinkPhotons = cms.string("(r9()>0.8 || chargedHadronIso()<20 || chargedHadronIso()<0.3*pt())"), #keep all associated clusters/rechits/conversions @@ -13,6 +16,7 @@ keepGsfElectrons = cms.string(""), #keep in output slimRelinkGsfElectrons = cms.string(""), #keep only slimmed SuperCluster plus seed cluster relinkGsfElectrons = cms.string("pt>5"), #keep all associated clusters/rechits/conversions + pflowSuperclusters = cms.InputTag("superClusterMerger"), photons = cms.InputTag("gedPhotons"), ootPhotons = cms.InputTag("ootPhotons"), gsfElectrons = cms.InputTag("gedGsfElectrons"), @@ -50,6 +54,13 @@ hcalHitSel = interestingEgammaIsoHCALSel ) +superClusterMerger = cms.EDProducer("EgammaSuperClusterMerger", + src = cms.VInputTag( + cms.InputTag("particleFlowSuperClusterECAL:particleFlowSuperClusterECALBarrel"), + cms.InputTag("particleFlowSuperClusterECAL:particleFlowSuperClusterECALEndcapWithPreshower"), + ) +) + from Configuration.Eras.Modifier_phase2_common_cff import phase2_common phase2_common.toModify(reducedEgamma, preshowerEcalHits = "", diff --git a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc index 56a87fd6abbbd..61f8e38125c9b 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc @@ -153,6 +153,7 @@ class ReducedEGProducer : public edm::stream::EDProducer<> { } //tokens for input collections + const edm::EDGetTokenT superclusterT_; const edm::EDGetTokenT photonT_; edm::EDGetTokenT ootPhotonT_; const edm::EDGetTokenT gsfElectronT_; @@ -178,6 +179,10 @@ class ReducedEGProducer : public edm::stream::EDProducer<> { const edm::EDGetTokenT recoHIPhotonIsolationMapInputToken_; edm::EDPutTokenT recoHIPhotonIsolationMapOutputName_; + const double scPtMin_; + const double scAbsetaMax_; + const bool relinkSupercluster_; + const bool applyPhotonCalibOnData_; const bool applyPhotonCalibOnMC_; const bool applyGsfElectronCalibOnData_; @@ -263,7 +268,8 @@ namespace { } // namespace ReducedEGProducer::ReducedEGProducer(const edm::ParameterSet& config) - : photonT_(consumes(config.getParameter("photons"))), + : superclusterT_(consumes(config.getParameter("pflowSuperclusters"))), + photonT_(consumes(config.getParameter("photons"))), gsfElectronT_(consumes(config.getParameter("gsfElectrons"))), conversionT_(consumes(config.getParameter("conversions"))), singleConversionT_(consumes(config.getParameter("singleConversions"))), @@ -280,6 +286,9 @@ ReducedEGProducer::ReducedEGProducer(const edm::ParameterSet& config) !config.getParameter("hiPhotonIsolationMapInput").label().empty() ? consumes(config.getParameter("hiPhotonIsolationMapInput")) : edm::EDGetTokenT{}}, + scPtMin_(config.getParameter("keepPfSuperclusterPtMin")), + scAbsetaMax_(config.getParameter("keepPfSuperclusterAbsetaMax")), + relinkSupercluster_(config.getParameter("relinkSupercluster")), //calibration flags applyPhotonCalibOnData_(config.getParameter("applyPhotonCalibOnData")), applyPhotonCalibOnMC_(config.getParameter("applyPhotonCalibOnMC")), @@ -380,6 +389,7 @@ void ReducedEGProducer::beginRun(edm::Run const& run, const edm::EventSetup& iSe void ReducedEGProducer::produce(edm::Event& event, const edm::EventSetup& eventSetup) { //get input collections + auto scHandle = event.getHandle(superclusterT_); auto photonHandle = event.getHandle(photonT_); auto ootPhotonHandle = @@ -682,6 +692,21 @@ void ReducedEGProducer::produce(edm::Event& event, const edm::EventSetup& eventS linkHcalHits(*gsfElectron.superCluster(), *hbheHitHandle, hcalRechitMap); } + //loop over input SuperClusters + index = -1; + for (const auto& superCluster : *scHandle) { + index++; + + if ( superCluster.energy()/std::cosh(superCluster.eta()) < scPtMin_ ) + continue; + + if ( std::abs(superCluster.eta()) > scAbsetaMax_ ) + continue; + + reco::SuperClusterRef superClusterRef(scHandle, index); + linkSuperCluster(superClusterRef, superClusterMap, superClusters, relinkSupercluster_, superClusterFullRelinkMap); + } + //loop over output SuperClusters and fill maps index = 0; for (auto& superCluster : superClusters) { From 9c6091bd8495d9b9f682004fa0f9566d16409586 Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Wed, 16 Oct 2024 01:06:35 +0200 Subject: [PATCH 02/10] differentiate supercluster Pt cut and linking Pt cut --- .../python/reducedEgamma_cfi.py | 2 +- .../EgammaPhotonProducers/src/ReducedEGProducer.cc | 12 ++++++++---- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py index 458f0f6868a7b..335965e57a9be 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py @@ -6,7 +6,7 @@ reducedEgamma = cms.EDProducer("ReducedEGProducer", keepPfSuperclusterPtMin = cms.double(5.), keepPfSuperclusterAbsetaMax = cms.double(2.5), - relinkSupercluster = cms.bool(True), + relinkSuperclusterPtMin = cms.double(10.), keepPhotons = cms.string("hadTowOverEm()<0.15 && pt>10 && (pt>14 || chargedHadronIso()<10)"), #keep in output slimRelinkPhotons = cms.string("hadTowOverEm()<0.15 && pt>10 && (pt>14 || chargedHadronIso()<10)"), #keep only slimmed SuperCluster plus seed cluster relinkPhotons = cms.string("(r9()>0.8 || chargedHadronIso()<20 || chargedHadronIso()<0.3*pt())"), #keep all associated clusters/rechits/conversions diff --git a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc index 61f8e38125c9b..1fb9e3d5f29a3 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc @@ -181,7 +181,7 @@ class ReducedEGProducer : public edm::stream::EDProducer<> { const double scPtMin_; const double scAbsetaMax_; - const bool relinkSupercluster_; + const double relinkSuperclusterPtMin_; const bool applyPhotonCalibOnData_; const bool applyPhotonCalibOnMC_; @@ -288,7 +288,7 @@ ReducedEGProducer::ReducedEGProducer(const edm::ParameterSet& config) : edm::EDGetTokenT{}}, scPtMin_(config.getParameter("keepPfSuperclusterPtMin")), scAbsetaMax_(config.getParameter("keepPfSuperclusterAbsetaMax")), - relinkSupercluster_(config.getParameter("relinkSupercluster")), + relinkSuperclusterPtMin_(config.getParameter("relinkSuperclusterPtMin")), //calibration flags applyPhotonCalibOnData_(config.getParameter("applyPhotonCalibOnData")), applyPhotonCalibOnMC_(config.getParameter("applyPhotonCalibOnMC")), @@ -697,14 +697,18 @@ void ReducedEGProducer::produce(edm::Event& event, const edm::EventSetup& eventS for (const auto& superCluster : *scHandle) { index++; - if ( superCluster.energy()/std::cosh(superCluster.eta()) < scPtMin_ ) + const double superclusPt = superCluster.energy()/std::cosh(superCluster.eta()); + + if ( superclusPt < scPtMin_ ) continue; if ( std::abs(superCluster.eta()) > scAbsetaMax_ ) continue; + bool relinkSupercluster = superclusPt > relinkSuperclusterPtMin_; + reco::SuperClusterRef superClusterRef(scHandle, index); - linkSuperCluster(superClusterRef, superClusterMap, superClusters, relinkSupercluster_, superClusterFullRelinkMap); + linkSuperCluster(superClusterRef, superClusterMap, superClusters, relinkSupercluster, superClusterFullRelinkMap); } //loop over output SuperClusters and fill maps From af375a2d22c5becd4a05130d8751f8b15bedce62 Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Wed, 16 Oct 2024 12:14:00 +0200 Subject: [PATCH 03/10] propagate SC collection to the EGM flavor nano --- .../plugins/SimpleFlatTableProducerPlugins.cc | 4 ++++ .../NanoAOD/python/egamma_custom_cff.py | 24 +++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc b/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc index 53fc6f16ee47e..c3c3b1e769575 100644 --- a/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc +++ b/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc @@ -6,6 +6,9 @@ typedef SimpleFlatTableProducer SimpleCandidateFlatTableProduce #include "DataFormats/TrackReco/interface/Track.h" typedef SimpleFlatTableProducer SimpleTrackFlatTableProducer; +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +typedef SimpleFlatTableProducer SimpleSuperclusterFlatTableProducer; + #include "DataFormats/JetReco/interface/PFJet.h" typedef SimpleFlatTableProducer SimplePFJetFlatTableProducer; @@ -48,6 +51,7 @@ typedef EventSingletonSimpleFlatTableProducer SimpleBeamspotFlat #include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(SimpleCandidateFlatTableProducer); DEFINE_FWK_MODULE(SimpleTrackFlatTableProducer); +DEFINE_FWK_MODULE(SimpleSuperclusterFlatTableProducer); DEFINE_FWK_MODULE(SimplePFJetFlatTableProducer); DEFINE_FWK_MODULE(SimpleVertexFlatTableProducer); DEFINE_FWK_MODULE(SimpleSecondaryVertexFlatTableProducer); diff --git a/PhysicsTools/NanoAOD/python/egamma_custom_cff.py b/PhysicsTools/NanoAOD/python/egamma_custom_cff.py index 176115a1b1cbf..af39bae36bf20 100644 --- a/PhysicsTools/NanoAOD/python/egamma_custom_cff.py +++ b/PhysicsTools/NanoAOD/python/egamma_custom_cff.py @@ -81,6 +81,24 @@ ) ) +superclusterTable = cms.EDProducer("SimpleSuperclusterFlatTableProducer", + src = cms.InputTag("reducedEgamma","reducedSuperClusters"), + name = cms.string("Supercluster"), + doc = cms.string("Supercluster collection"), + variables = cms.PSet( + energy = Var("energy()",float,doc="supercluster energy",precision=10), + eta = Var("eta()",float,doc="supercluster eta",precision=10), + phi = Var("phi()",float,doc="supercluster phi",precision=10), + rawEnergy = Var("rawEnergy()",float,doc="sum of basic clusters energy",precision=10), + preshowerEnergy = Var("preshowerEnergy()",float,doc="sum of energy in preshower",precision=10), + etaWidth = Var("etaWidth()",float,doc="supercluster eta width",precision=10), + phiWidth = Var("etaWidth()",float,doc="supercluster phi width",precision=10), + seedClusEnergy = Var("seed().energy()",float,doc="seed cluster energy",precision=10), + seedClusterEta = Var("seed().eta()",float,doc="seed cluster eta",precision=10), + seedClusterPhi = Var("seed().phi()",float,doc="seed cluster phi",precision=10), + ) +) + def addExtraEGammaVarsCustomize(process): #photon process.finalPhotons.cut = cms.string("pt > 1 ") @@ -97,5 +115,11 @@ def addExtraEGammaVarsCustomize(process): if hasattr(process,'nanoDQM'): process.nanoDQM.vplots.Electron.plots = _Electron_extra_plots + + #superCluster + process.superclusterTable = superclusterTable + + process.superclusterTask = cms.Task(process.superclusterTable) + process.nanoTableTaskCommon.add(process.superclusterTask) return process From d0d570d56279b4a72feed80420f11ee4c4ac1e36 Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Thu, 17 Oct 2024 14:28:26 +0200 Subject: [PATCH 04/10] minimize the file size increase of miniAOD --- RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py index 335965e57a9be..6be63221991a8 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py @@ -4,9 +4,9 @@ from RecoJets.Configuration.CaloTowersES_cfi import * reducedEgamma = cms.EDProducer("ReducedEGProducer", - keepPfSuperclusterPtMin = cms.double(5.), + keepPfSuperclusterPtMin = cms.double(10.), keepPfSuperclusterAbsetaMax = cms.double(2.5), - relinkSuperclusterPtMin = cms.double(10.), + relinkSuperclusterPtMin = cms.double(99999.), # no SC linking keepPhotons = cms.string("hadTowOverEm()<0.15 && pt>10 && (pt>14 || chargedHadronIso()<10)"), #keep in output slimRelinkPhotons = cms.string("hadTowOverEm()<0.15 && pt>10 && (pt>14 || chargedHadronIso()<10)"), #keep only slimmed SuperCluster plus seed cluster relinkPhotons = cms.string("(r9()>0.8 || chargedHadronIso()<20 || chargedHadronIso()<0.3*pt())"), #keep all associated clusters/rechits/conversions From 3f7ea8287ccd934f17c97dd86dd86b8a24dc1d59 Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Thu, 17 Oct 2024 14:58:44 +0200 Subject: [PATCH 05/10] replace the ratio of corr/uncorr variable to the uncorrected pt (to avoid issues like #46046) --- PhysicsTools/NanoAOD/python/electrons_cff.py | 2 +- PhysicsTools/NanoAOD/python/photons_cff.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/PhysicsTools/NanoAOD/python/electrons_cff.py b/PhysicsTools/NanoAOD/python/electrons_cff.py index 0ae00ce1e79d4..7bc2fb06467c3 100644 --- a/PhysicsTools/NanoAOD/python/electrons_cff.py +++ b/PhysicsTools/NanoAOD/python/electrons_cff.py @@ -433,7 +433,7 @@ def _get_bitmapVIDForEle_docstring(modules,WorkingPoints): electronTable.variables, pt = Var("pt*userFloat('ecalTrkEnergyPostCorrNew')/userFloat('ecalTrkEnergyPreCorrNew')", float, precision=-1, doc="p_{T}"), energyErr = Var("userFloat('ecalTrkEnergyErrPostCorrNew')", float, precision=6, doc="energy error of the cluster-track combination"), - eCorr = Var("userFloat('ecalTrkEnergyPostCorrNew')/userFloat('ecalTrkEnergyPreCorrNew')", float, doc="ratio of the calibrated energy/miniaod energy"), + ptPreCorr = Var("pt", float, doc="pt of the electron before energy corrections"), scEtOverPt = Var("(superCluster().energy()/(pt*userFloat('ecalTrkEnergyPostCorrNew')/userFloat('ecalTrkEnergyPreCorrNew')*cosh(superCluster().eta())))-1",float,doc="(supercluster transverse energy)/pt-1",precision=8), dEscaleUp=Var("userFloat('ecalTrkEnergyPostCorrNew')-userFloat('energyScaleUpNew')", float, doc="ecal energy scale shifted 1 sigma up(adding gain/stat/syst in quadrature)", precision=8), dEscaleDown=Var("userFloat('ecalTrkEnergyPostCorrNew')-userFloat('energyScaleDownNew')", float, doc="ecal energy scale shifted 1 sigma down (adding gain/stat/syst in quadrature)", precision=8), diff --git a/PhysicsTools/NanoAOD/python/photons_cff.py b/PhysicsTools/NanoAOD/python/photons_cff.py index be2c620dba4ed..bd9ce7ebd0c2d 100644 --- a/PhysicsTools/NanoAOD/python/photons_cff.py +++ b/PhysicsTools/NanoAOD/python/photons_cff.py @@ -283,7 +283,7 @@ def make_bitmapVID_docstring(id_modules_working_points_pset): photonTable.variables, pt = Var("pt*userFloat('ecalEnergyPostCorrNew')/userFloat('ecalEnergyPreCorrNew')", float, precision=-1, doc="p_{T}"), energyErr = Var("userFloat('ecalEnergyErrPostCorrNew')",float,doc="energy error of the cluster from regression",precision=6), - eCorr = Var("userFloat('ecalEnergyPostCorrNew')/userFloat('ecalEnergyPreCorrNew')",float,doc="ratio of the calibrated energy/miniaod energy"), + ptPreCorr = Var("pt",float,doc="pt of the photon before energy corrections"), cutBased = Var( "userInt('cutBasedID_Fall17V2_loose')+userInt('cutBasedID_Fall17V2_medium')+userInt('cutBasedID_Fall17V2_tight')", "uint8", From d4762a9632816174d337167d1ed82a3a676f998d Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Thu, 17 Oct 2024 15:13:38 +0200 Subject: [PATCH 06/10] apply code-checks --- RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc index 1fb9e3d5f29a3..7afae9aab9dbf 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc @@ -697,12 +697,12 @@ void ReducedEGProducer::produce(edm::Event& event, const edm::EventSetup& eventS for (const auto& superCluster : *scHandle) { index++; - const double superclusPt = superCluster.energy()/std::cosh(superCluster.eta()); + const double superclusPt = superCluster.energy() / std::cosh(superCluster.eta()); - if ( superclusPt < scPtMin_ ) + if (superclusPt < scPtMin_) continue; - if ( std::abs(superCluster.eta()) > scAbsetaMax_ ) + if (std::abs(superCluster.eta()) > scAbsetaMax_) continue; bool relinkSupercluster = superclusPt > relinkSuperclusterPtMin_; From 3e9b4b40c60be1b5db34a435bda9b9685d4305cc Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Sat, 19 Oct 2024 17:54:19 +0200 Subject: [PATCH 07/10] reflect changes in #3f7ea82 to nanoDQM plots --- PhysicsTools/NanoAOD/python/nanoDQM_cff.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PhysicsTools/NanoAOD/python/nanoDQM_cff.py b/PhysicsTools/NanoAOD/python/nanoDQM_cff.py index 83133f9c88104..6772cd3e72bc6 100644 --- a/PhysicsTools/NanoAOD/python/nanoDQM_cff.py +++ b/PhysicsTools/NanoAOD/python/nanoDQM_cff.py @@ -123,7 +123,7 @@ Plot1D('dEscaleDown', 'dEscaleDown', 100, -0.01, 0.01, '#Delta E scaleDown'), Plot1D('dEsigmaUp', 'dEsigmaUp', 100, -0.1, 0.1, '#Delta E sigmaUp'), Plot1D('dEsigmaDown', 'dEsigmaDown', 100, -0.1, 0.1, '#Delta E sigmaDown'), - Plot1D('eCorr', 'eCorr', 20, 0.8, 1.2, 'ratio of the calibrated energy/miniaod energy'), + Plot1D('ptPreCorr', 'ptPreCorr', 100, 0., 500., 'Pt before scale & smearing energy corrections'), ]) run2_egamma.toModify( nanoDQM.vplots.Electron, @@ -146,7 +146,7 @@ def _match(name): Plot1D('dEscaleDown', 'dEscaleDown', 100, -0.01, 0.01, '#Delta E scaleDown'), Plot1D('dEsigmaUp', 'dEsigmaUp', 100, -0.1, 0.1, '#Delta E sigmaUp'), Plot1D('dEsigmaDown', 'dEsigmaDown', 100, -0.1, 0.1, '#Delta E sigmaDown'), - Plot1D('eCorr', 'eCorr', 20, 0.8, 1.2, 'ratio of the calibrated energy/miniaod energy'), + Plot1D('ptPreCorr', 'ptPreCorr', 100, 0., 500., 'Pt before scale & smearing energy corrections'), ]) run2_egamma.toModify( nanoDQM.vplots.Photon, From 7275316cc5a35bdd5f472cd644d7e2e5ea20fc8f Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Sun, 20 Oct 2024 01:55:59 +0200 Subject: [PATCH 08/10] incorporate SC trk isolation from #42007 --- .../NanoAOD/python/egamma_custom_cff.py | 8 +- .../interface/EleTkIsolFromCands.h | 12 +- .../interface/SuperclusTkIsolFromCands.h | 27 ++++ .../plugins/SuperclusValueMapProducer.cc | 128 ++++++++++++++++++ .../python/superclusValueMapProducer_cfi.py | 33 +++++ .../src/SuperclusTkIsolFromCands.cc | 32 +++++ 6 files changed, 233 insertions(+), 7 deletions(-) create mode 100644 RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h create mode 100644 RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc create mode 100644 RecoEgamma/EgammaIsolationAlgos/python/superclusValueMapProducer_cfi.py create mode 100644 RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc diff --git a/PhysicsTools/NanoAOD/python/egamma_custom_cff.py b/PhysicsTools/NanoAOD/python/egamma_custom_cff.py index af39bae36bf20..a47e6484bfb1e 100644 --- a/PhysicsTools/NanoAOD/python/egamma_custom_cff.py +++ b/PhysicsTools/NanoAOD/python/egamma_custom_cff.py @@ -8,6 +8,7 @@ from PhysicsTools.NanoAOD.nanoDQM_cfi import nanoDQM from PhysicsTools.NanoAOD.nanoDQM_cff import _Photon_extra_plots, _Electron_extra_plots from PhysicsTools.NanoAOD.triggerObjects_cff import triggerObjectTable, mksel +from RecoEgamma.EgammaIsolationAlgos.superclusValueMapProducer_cfi import superclusValueMaps customElectronFilterBits = cms.PSet( doc = cms.string("PixelMatched e/gamma"), # this may also select photons! @@ -96,6 +97,9 @@ seedClusEnergy = Var("seed().energy()",float,doc="seed cluster energy",precision=10), seedClusterEta = Var("seed().eta()",float,doc="seed cluster eta",precision=10), seedClusterPhi = Var("seed().phi()",float,doc="seed cluster phi",precision=10), + ), + externalVariables = cms.PSet( + trkIso = ExtVar("superclusValueMaps:superclusTkIso",float,doc="supercluster track iso within 0.06 < dR < 0.4 & |dEta| > 0.03",precision=10), ) ) @@ -117,9 +121,11 @@ def addExtraEGammaVarsCustomize(process): process.nanoDQM.vplots.Electron.plots = _Electron_extra_plots #superCluster + process.superclusValueMaps = superclusValueMaps process.superclusterTable = superclusterTable - process.superclusterTask = cms.Task(process.superclusterTable) + process.superclusterTask = cms.Task(process.superclusValueMaps) + process.superclusterTask.add(process.superclusterTable) process.nanoTableTaskCommon.add(process.superclusterTask) return process diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h b/RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h index c8695909a214d..db0c96bb7fa24 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h @@ -92,6 +92,12 @@ class EleTkIsolFromCands { Output operator()(const reco::TrackBase& electronTrack); +protected: + using TrackTable = edm::soa::Table; + TrackTable const& getPreselectedTracks(bool isBarrel); + + Configuration const& cfg_; + private: // For each electron, we want to try out which tracks are in a cone around // it. However, this will get expensive if there are many electrons and @@ -103,8 +109,6 @@ class EleTkIsolFromCands { // the electron. Note that this has to be done twice, because the required // preselection is different for barrel and endcap electrons. - using TrackTable = edm::soa::Table; - static bool passPIDVeto(const int pdgId, const EleTkIsolFromCands::PIDVeto pidVeto); static TrackTable preselectTracks(reco::TrackCollection const& tracks, TrkCuts const& cuts); @@ -118,10 +122,6 @@ class EleTkIsolFromCands { static bool passQual(const reco::TrackBase& trk, const std::vector& quals); static bool passAlgo(const reco::TrackBase& trk, const std::vector& algosToRej); - TrackTable const& getPreselectedTracks(bool isBarrel); - - Configuration const& cfg_; - // All of these member variables are related to the caching of preselected tracks reco::TrackCollection const* tracks_ = nullptr; pat::PackedCandidateCollection const* cands_ = nullptr; diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h b/RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h new file mode 100644 index 0000000000000..b05a49ec40d11 --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h @@ -0,0 +1,27 @@ +#ifndef RecoEgamma_EgammaIsolationAlgos_SuperclusTkIsolFromCands_h +#define RecoEgamma_EgammaIsolationAlgos_SuperclusTkIsolFromCands_h + +#include "CommonTools/Utils/interface/KinematicColumns.h" +#include "DataFormats/TrackReco/interface/TrackBase.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "DataFormats/Math/interface/Point3D.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "FWCore/SOA/interface/Table.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h" + +class SuperclusTkIsolFromCands : public EleTkIsolFromCands { +public: + explicit SuperclusTkIsolFromCands(Configuration const& cfg, reco::TrackCollection const& tracks) + : EleTkIsolFromCands(cfg,tracks) {} + explicit SuperclusTkIsolFromCands(Configuration const& cfg, + pat::PackedCandidateCollection const& cands, + PIDVeto pidVeto = PIDVeto::NONE) + : EleTkIsolFromCands(cfg,cands,pidVeto) {} + + Output operator() (const reco::SuperCluster& sc, const math::XYZPoint& vtx); +}; + +#endif diff --git a/RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc b/RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc new file mode 100644 index 0000000000000..94e36c2a9544d --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc @@ -0,0 +1,128 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/Common/interface/ValueMap.h" +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Math/interface/Point3D.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" + +#include "RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h" + +// copy-pasted from ElectronHEEPIDValueMapProducer + +class SuperclusValueMapProducer : public edm::stream::EDProducer<> { +public: + explicit SuperclusValueMapProducer(const edm::ParameterSet&); + ~SuperclusValueMapProducer()=default; + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + std::vector> setTokens(const std::vector& tags); + + template + static void writeValueMap(edm::Event& iEvent, + const edm::Handle>& handle, + const std::vector& values, + const std::string& label); + + std::vector candVetos_; + + const std::vector candTags_; + const std::vector> candTokens_; + const edm::EDGetTokenT> scToken_; + const edm::EDGetTokenT> pvToken_; + const edm::EDGetTokenT bsToken_; + + const SuperclusTkIsolFromCands::Configuration trkIsoCalcCfg_; + + const std::string superclusTkIsoLabel_ = "superclusTkIso"; +}; + +SuperclusValueMapProducer::SuperclusValueMapProducer(const edm::ParameterSet& iConfig) +: candTags_(iConfig.getParameter>("cands")), + candTokens_(setTokens(candTags_)), + scToken_(consumes>(iConfig.getParameter("srcSc"))), + pvToken_(consumes>(iConfig.getParameter("srcPv"))), + bsToken_(consumes(iConfig.getParameter("srcBs"))), + trkIsoCalcCfg_(iConfig.getParameter("trkIsoConfig")) +{ + auto fillVetos = [] (const auto& in, auto& out) { + std::transform(in.begin(), in.end(), std::back_inserter(out), SuperclusTkIsolFromCands::pidVetoFromStr); + }; + + fillVetos(iConfig.getParameter>("candVetos"), candVetos_); + + if (candVetos_.size() != candTags_.size()) + throw cms::Exception("ConfigError") << "Error candVetos should be the same size as cands" << std::endl; + + produces>(superclusTkIsoLabel_); +} + +void SuperclusValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + edm::Handle> scHandle; + iEvent.getByToken(scToken_,scHandle); + + edm::Handle> pvHandle; + iEvent.getByToken(pvToken_,pvHandle); + + edm::Handle bsHandle; + iEvent.getByToken(bsToken_,bsHandle); + + math::XYZPoint pos; + + if (pvHandle.isValid() && !pvHandle->empty()) + pos = pvHandle->front().position(); // first try PV + else + pos = (*bsHandle).position(); // fall back to BS + + std::vector> candHandles(candTokens_.size()); + std::vector> tkIsoCalc; + + for (unsigned idx=0; idx(trkIsoCalcCfg_,*(candHandles.at(idx)),candVetos_.at(idx))); + } + + std::vector vecTkIso; + vecTkIso.reserve(scHandle->size()); + + for (const auto& sc : *scHandle) { + float tkIso = 0.; + + for (auto& calc : tkIsoCalc) + tkIso += (*calc)(sc,pos).ptSum; + + vecTkIso.push_back(tkIso); + } + + writeValueMap(iEvent,scHandle,vecTkIso,superclusTkIsoLabel_); +} + +std::vector> SuperclusValueMapProducer::setTokens(const std::vector& tags) { + std::vector> out; + + for (const auto& tag : tags) + out.push_back(consumes(tag)); + + return out; +} + +template +void SuperclusValueMapProducer::writeValueMap(edm::Event& iEvent, + const edm::Handle>& handle, + const std::vector& values, + const std::string& label) { + std::unique_ptr> valMap(new edm::ValueMap()); + typename edm::ValueMap::Filler filler(*valMap); + filler.insert(handle, values.begin(), values.end()); + filler.fill(); + iEvent.put(std::move(valMap), label); +} + +DEFINE_FWK_MODULE(SuperclusValueMapProducer); diff --git a/RecoEgamma/EgammaIsolationAlgos/python/superclusValueMapProducer_cfi.py b/RecoEgamma/EgammaIsolationAlgos/python/superclusValueMapProducer_cfi.py new file mode 100644 index 0000000000000..5d4fe1a51ff8f --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/python/superclusValueMapProducer_cfi.py @@ -0,0 +1,33 @@ +import FWCore.ParameterSet.Config as cms + +from RecoEgamma.EgammaIsolationAlgos.electronTrackIsolations_cfi import trkIsol04CfgV2 + +# supposed to be calculated from AOD (see https://github.com/cms-sw/cmssw/pull/42007) +# but since we didn't create the new PAT version of SC +# trk iso values are calculated on-the-fly from miniAOD +# parameters are inherited from HEEP trk iso + +# instead, we use fairly loose inner cone & strip veto from EgammaHLTTrackIsolation +# to avoid strong bias due to the tight trk isolation +# e.g. see hltEgammaHollowTrackIsoL1Seeded filter + +scTrkIso04 = cms.PSet( + barrelCuts = trkIsol04CfgV2.barrelCuts.clone( + minDR = 0.06, + minDEta = 0.03 + ), + endcapCuts = trkIsol04CfgV2.endcapCuts.clone( + minDR = 0.06, + minDEta = 0.03 + ) +) + +superclusValueMaps = cms.EDProducer("SuperclusValueMapProducer", + srcBs = cms.InputTag("offlineBeamSpot"), + srcPv = cms.InputTag("offlineSlimmedPrimaryVertices"), + srcSc = cms.InputTag("reducedEgamma:reducedSuperClusters"), + cands = cms.VInputTag("packedPFCandidates", + "lostTracks"), # do not count electron tracks in the trk iso + candVetos = cms.vstring("ELES","NONE"), + trkIsoConfig = scTrkIso04, +) diff --git a/RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc b/RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc new file mode 100644 index 0000000000000..391509e79e2bf --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc @@ -0,0 +1,32 @@ +#include "RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/Math/interface/deltaR.h" + +SuperclusTkIsolFromCands::Output SuperclusTkIsolFromCands::operator() (const reco::SuperCluster& sc, const math::XYZPoint& vtx) { + using namespace edm::soa::col; + + float ptSum = 0.; + int nrTrks = 0; + + const float scEta = sc.eta(); + const float scPhi = sc.phi(); + const float vtxVz = vtx.z(); + + const bool isBarrelSC = std::abs(scEta) < 1.5; + + auto const& preselectedTracks = getPreselectedTracks(isBarrelSC); + auto const& cuts = isBarrelSC ? cfg_.barrelCuts : cfg_.endcapCuts; + + for (auto const& trk : preselectedTracks) { + const float dR2 = reco::deltaR2(scEta, scPhi, trk.get(), trk.get()); + const float dEta = trk.get() - scEta; + const float dZ = vtxVz - trk.get(); + + if (dR2 >= cuts.minDR2 && dR2 <= cuts.maxDR2 && std::abs(dEta) >= cuts.minDEta && std::abs(dZ) < cuts.maxDZ) { + ptSum += trk.get(); + nrTrks++; + } + } + + return {nrTrks, ptSum}; +} From 6340fac927b62e956587d473f6839261de33f431 Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Sun, 20 Oct 2024 02:34:30 +0200 Subject: [PATCH 09/10] apply code-checks --- .../interface/SuperclusTkIsolFromCands.h | 10 ++--- .../plugins/SuperclusValueMapProducer.cc | 43 ++++++++++--------- .../src/SuperclusTkIsolFromCands.cc | 3 +- 3 files changed, 29 insertions(+), 27 deletions(-) diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h b/RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h index b05a49ec40d11..5b4e47bcc4e33 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h @@ -15,13 +15,13 @@ class SuperclusTkIsolFromCands : public EleTkIsolFromCands { public: explicit SuperclusTkIsolFromCands(Configuration const& cfg, reco::TrackCollection const& tracks) - : EleTkIsolFromCands(cfg,tracks) {} + : EleTkIsolFromCands(cfg, tracks) {} explicit SuperclusTkIsolFromCands(Configuration const& cfg, - pat::PackedCandidateCollection const& cands, - PIDVeto pidVeto = PIDVeto::NONE) - : EleTkIsolFromCands(cfg,cands,pidVeto) {} + pat::PackedCandidateCollection const& cands, + PIDVeto pidVeto = PIDVeto::NONE) + : EleTkIsolFromCands(cfg, cands, pidVeto) {} - Output operator() (const reco::SuperCluster& sc, const math::XYZPoint& vtx); + Output operator()(const reco::SuperCluster& sc, const math::XYZPoint& vtx); }; #endif diff --git a/RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc b/RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc index 94e36c2a9544d..28ee8eac7ebe3 100644 --- a/RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc +++ b/RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc @@ -18,10 +18,10 @@ class SuperclusValueMapProducer : public edm::stream::EDProducer<> { public: explicit SuperclusValueMapProducer(const edm::ParameterSet&); - ~SuperclusValueMapProducer()=default; + ~SuperclusValueMapProducer() = default; private: - void produce(edm::Event&, const edm::EventSetup&) override; + void produce(edm::Event&, const edm::EventSetup&) override; std::vector> setTokens(const std::vector& tags); @@ -45,14 +45,13 @@ class SuperclusValueMapProducer : public edm::stream::EDProducer<> { }; SuperclusValueMapProducer::SuperclusValueMapProducer(const edm::ParameterSet& iConfig) -: candTags_(iConfig.getParameter>("cands")), - candTokens_(setTokens(candTags_)), - scToken_(consumes>(iConfig.getParameter("srcSc"))), - pvToken_(consumes>(iConfig.getParameter("srcPv"))), - bsToken_(consumes(iConfig.getParameter("srcBs"))), - trkIsoCalcCfg_(iConfig.getParameter("trkIsoConfig")) -{ - auto fillVetos = [] (const auto& in, auto& out) { + : candTags_(iConfig.getParameter>("cands")), + candTokens_(setTokens(candTags_)), + scToken_(consumes>(iConfig.getParameter("srcSc"))), + pvToken_(consumes>(iConfig.getParameter("srcPv"))), + bsToken_(consumes(iConfig.getParameter("srcBs"))), + trkIsoCalcCfg_(iConfig.getParameter("trkIsoConfig")) { + auto fillVetos = [](const auto& in, auto& out) { std::transform(in.begin(), in.end(), std::back_inserter(out), SuperclusTkIsolFromCands::pidVetoFromStr); }; @@ -66,27 +65,28 @@ SuperclusValueMapProducer::SuperclusValueMapProducer(const edm::ParameterSet& iC void SuperclusValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { edm::Handle> scHandle; - iEvent.getByToken(scToken_,scHandle); + iEvent.getByToken(scToken_, scHandle); edm::Handle> pvHandle; - iEvent.getByToken(pvToken_,pvHandle); + iEvent.getByToken(pvToken_, pvHandle); edm::Handle bsHandle; - iEvent.getByToken(bsToken_,bsHandle); + iEvent.getByToken(bsToken_, bsHandle); math::XYZPoint pos; if (pvHandle.isValid() && !pvHandle->empty()) - pos = pvHandle->front().position(); // first try PV + pos = pvHandle->front().position(); // first try PV else - pos = (*bsHandle).position(); // fall back to BS + pos = (*bsHandle).position(); // fall back to BS std::vector> candHandles(candTokens_.size()); std::vector> tkIsoCalc; - for (unsigned idx=0; idx(trkIsoCalcCfg_,*(candHandles.at(idx)),candVetos_.at(idx))); + for (unsigned idx = 0; idx < candTokens_.size(); idx++) { + iEvent.getByToken(candTokens_.at(idx), candHandles.at(idx)); + tkIsoCalc.push_back( + std::make_unique(trkIsoCalcCfg_, *(candHandles.at(idx)), candVetos_.at(idx))); } std::vector vecTkIso; @@ -96,15 +96,16 @@ void SuperclusValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetu float tkIso = 0.; for (auto& calc : tkIsoCalc) - tkIso += (*calc)(sc,pos).ptSum; + tkIso += (*calc)(sc, pos).ptSum; vecTkIso.push_back(tkIso); } - writeValueMap(iEvent,scHandle,vecTkIso,superclusTkIsoLabel_); + writeValueMap(iEvent, scHandle, vecTkIso, superclusTkIsoLabel_); } -std::vector> SuperclusValueMapProducer::setTokens(const std::vector& tags) { +std::vector> SuperclusValueMapProducer::setTokens( + const std::vector& tags) { std::vector> out; for (const auto& tag : tags) diff --git a/RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc b/RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc index 391509e79e2bf..2168fb848111a 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc @@ -2,7 +2,8 @@ #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/Math/interface/deltaR.h" -SuperclusTkIsolFromCands::Output SuperclusTkIsolFromCands::operator() (const reco::SuperCluster& sc, const math::XYZPoint& vtx) { +SuperclusTkIsolFromCands::Output SuperclusTkIsolFromCands::operator()(const reco::SuperCluster& sc, + const math::XYZPoint& vtx) { using namespace edm::soa::col; float ptSum = 0.; From b2c2152cba2e4c352461df155c8721decfaa1c67 Mon Sep 17 00:00:00 2001 From: Sanghyun Ko Date: Mon, 21 Oct 2024 15:47:06 +0200 Subject: [PATCH 10/10] lower the SC pt threshold to accommodate a possible future development --- RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py index 6be63221991a8..8ab1a376b47bf 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py @@ -4,7 +4,7 @@ from RecoJets.Configuration.CaloTowersES_cfi import * reducedEgamma = cms.EDProducer("ReducedEGProducer", - keepPfSuperclusterPtMin = cms.double(10.), + keepPfSuperclusterPtMin = cms.double(5.), keepPfSuperclusterAbsetaMax = cms.double(2.5), relinkSuperclusterPtMin = cms.double(99999.), # no SC linking keepPhotons = cms.string("hadTowOverEm()<0.15 && pt>10 && (pt>14 || chargedHadronIso()<10)"), #keep in output