diff --git a/DataFormats/EgammaReco/interface/SlimmedSuperCluster.h b/DataFormats/EgammaReco/interface/SlimmedSuperCluster.h new file mode 100644 index 0000000000000..516c83beac215 --- /dev/null +++ b/DataFormats/EgammaReco/interface/SlimmedSuperCluster.h @@ -0,0 +1,51 @@ +#ifndef DataFormats_PatCandidates_SuperCluster_h +#define DataFormats_PatCandidates_SuperCluster_h + +/*** + +This is a slimmed data format to be able to store all EGamma superclusters +in MiniAOD as well as the necessary ID variables for EG workflows + +The goal is to catch the superclusters that dont make it to electrons/photons (of which there are many) to enable efficiency measurements of this stage. + +So initially it was thought to make this a LeafCandidate but a +supercluster is really a point + energy not a momentum so while more awkward not to just treat it as a p4, its more accurate to treat it as a position + energy + +author: Sam Harper, Swagata Mukherjee + +***/ + +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/Math/interface/Point3D.h" +#include + +namespace reco { + class SuperCluster; +} + +namespace reco { + + class SlimmedSuperCluster { + public: + SlimmedSuperCluster() : rawEnergy_(0.), preshowerEnergy_(0.), trkIso_(0.) {} + SlimmedSuperCluster(const reco::SuperCluster&, float trkIso = 0.); + float rawEnergy() const { return rawEnergy_; } + float preshowerEnergy() const { return preshowerEnergy_; } + DetId seedId() const { return clusterSeedIds_.empty() != 0 ? DetId(0) : clusterSeedIds_.front(); } + const std::vector& clusterSeedIds() const { return clusterSeedIds_; } + math::XYZPoint position() const; + float trkIso() const { return trkIso_; } + void setTrkIso(float val) { trkIso_ = val; } + + private: + float correctedEnergy_; + float rawEnergy_; + float preshowerEnergy_; + float rho_; + float eta_; + float phi_; + std::vector clusterSeedIds_; //the first one is always the seed + float trkIso_; + }; +} // namespace reco +#endif diff --git a/DataFormats/EgammaReco/interface/SlimmedSuperClusterFwd.h b/DataFormats/EgammaReco/interface/SlimmedSuperClusterFwd.h new file mode 100644 index 0000000000000..5c2469111a85d --- /dev/null +++ b/DataFormats/EgammaReco/interface/SlimmedSuperClusterFwd.h @@ -0,0 +1,25 @@ +#ifndef EgammaReco_SlimmedSuperClusterFwd_h +#define EgammaReco_SlimmedSuperClusterFwd_h +#include +#include "DataFormats/Common/interface/Ref.h" +#include "DataFormats/Common/interface/RefVector.h" +#include "DataFormats/Common/interface/RefProd.h" + +namespace reco { + class SlimmedSuperCluster; + + /// collection of SuperCluser objectr + typedef std::vector SlimmedSuperClusterCollection; + + /// reference to an object in a collection of SlimmedSuperCluster objects + typedef edm::Ref SlimmedSuperClusterRef; + + /// reference to a collection of SlimmedSuperCluster objects + typedef edm::RefProd SlimmedSuperClusterRefProd; + + /// vector of references to objects in the same colletion of SlimmedSuperCluster objects + typedef edm::RefVector SlimmedSuperClusterRefVector; + +} // namespace reco + +#endif diff --git a/DataFormats/EgammaReco/src/SlimmedSuperCluster.cc b/DataFormats/EgammaReco/src/SlimmedSuperCluster.cc new file mode 100644 index 0000000000000..9ba32779f766a --- /dev/null +++ b/DataFormats/EgammaReco/src/SlimmedSuperCluster.cc @@ -0,0 +1,27 @@ +#include "DataFormats/EgammaReco/interface/SlimmedSuperCluster.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "Math/GenVector/PositionVector3D.h" + +reco::SlimmedSuperCluster::SlimmedSuperCluster(const reco::SuperCluster& sc, float trkIso) + : correctedEnergy_(sc.correctedEnergy()), + rawEnergy_(sc.rawEnergy()), + preshowerEnergy_(sc.preshowerEnergy()), + rho_(sc.position().rho()), + eta_(sc.eta()), + phi_(sc.phi()), + trkIso_(trkIso) { + clusterSeedIds_.push_back(sc.seed()->seed()); + for (const auto& clus : sc.clusters()) { + if (clus != sc.seed()) { + clusterSeedIds_.push_back(clus->seed()); + } + } +} + +math::XYZPoint reco::SlimmedSuperCluster::position() const { + ROOT::Math::PositionVector3D> pos; + pos.SetRho(rho_); + pos.SetEta(eta_); + pos.SetPhi(phi_); + return math::XYZPoint(pos.x(), pos.y(), pos.z()); +} diff --git a/DataFormats/EgammaReco/src/classes.h b/DataFormats/EgammaReco/src/classes.h index 943ef0bd6a5be..5e135180b7dca 100644 --- a/DataFormats/EgammaReco/src/classes.h +++ b/DataFormats/EgammaReco/src/classes.h @@ -7,6 +7,7 @@ #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h" #include "DataFormats/EgammaReco/interface/PreshowerClusterShape.h" #include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/EgammaReco/interface/SlimmedSuperCluster.h" #include "Rtypes.h" #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" #include "DataFormats/EgammaReco/interface/ClusterShape.h" diff --git a/DataFormats/EgammaReco/src/classes_def.xml b/DataFormats/EgammaReco/src/classes_def.xml index 41b86c005dfb1..80ef32d2082a4 100644 --- a/DataFormats/EgammaReco/src/classes_def.xml +++ b/DataFormats/EgammaReco/src/classes_def.xml @@ -159,5 +159,19 @@ + + + + + + + + + + + + + + diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h b/RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h index baa545125d73f..a882a5a74073d 100644 --- a/RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h +++ b/RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h @@ -38,9 +38,18 @@ #include "DataFormats/Math/interface/Point3D.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" class EgammaHLTTrackIsolation { public: + EgammaHLTTrackIsolation(const edm::ParameterSet& iConfig) + : ptMin(iConfig.getParameter("ptMin")), + conesize(iConfig.getParameter("coneSize")), + zspan(iConfig.getParameter("zspan")), + rspan(iConfig.getParameter("rspan")), + vetoConesize(iConfig.getParameter("vetoConeSize")), + stripBarrel(iConfig.getParameter("stripBarrel")), + stripEndcap(iConfig.getParameter("stripEndcap")) {} EgammaHLTTrackIsolation(double egTrkIso_PtMin, double egTrkIso_ConeSize, double egTrkIso_ZSpan, @@ -54,45 +63,46 @@ class EgammaHLTTrackIsolation { rspan(egTrkIso_RSpan), vetoConesize(egTrkIso_VetoConeSize), stripBarrel(egTrkIso_stripBarrel), - stripEndcap(egTrkIso_stripEndcap) { - /* - std::cout << "EgammaHLTTrackIsolation instance:" - << " ptMin=" << ptMin << " " - << " conesize="<< conesize << " " - << " zspan=" << zspan << " " - << " rspan=" << rspan << " " - << " vetoConesize="<< vetoConesize - << std::endl; - */ - } + stripEndcap(egTrkIso_stripEndcap) {} /// Get number of tracks and Pt sum of tracks inside an isolation cone for electrons - std::pair electronIsolation(const reco::Track* const tr, const reco::TrackCollection* isoTracks); + std::pair electronIsolation(const reco::Track* const tr, const reco::TrackCollection* isoTracks) const; std::pair electronIsolation(const reco::Track* const tr, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks); + const reco::TrackCollection* isoTracks) const; std::pair electronIsolation(const reco::Track* const tr, const reco::TrackCollection* isoTracks, - GlobalPoint vertex); + GlobalPoint vertex) const; /// Get number of tracks and Pt sum of tracks inside an isolation cone for photons /// set useVertex=true to use PhotonCandidate vertex from EgammaPhotonVtxFinder /// set useVertex=false to consider all tracks for isolation std::pair photonIsolation(const reco::RecoCandidate* const recocand, const reco::TrackCollection* isoTracks, - bool useVertex); + bool useVertex) const; std::pair photonIsolation(const reco::RecoCandidate* const recocand, const reco::TrackCollection* isoTracks, - GlobalPoint vertex); + GlobalPoint vertex) const; + std::pair photonIsolation(const reco::RecoCandidate::Point& pos, + const reco::TrackCollection* isoTracks, + GlobalPoint vertex = GlobalPoint(0, 0, 0)) const; + std::pair photonIsolation(const reco::RecoCandidate* const recocand, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks); + const reco::TrackCollection* isoTracks) const; + + std::pair photonIsolation(float phoEta, + float phiPhi, + const reco::ElectronCollection* allEle, + const reco::TrackCollection* isoTracks) const; /// Get number of tracks inside an isolation cone for electrons - int electronTrackCount(const reco::Track* const tr, const reco::TrackCollection* isoTracks) { + int electronTrackCount(const reco::Track* const tr, const reco::TrackCollection* isoTracks) const { return electronIsolation(tr, isoTracks).first; } - int electronTrackCount(const reco::Track* const tr, const reco::TrackCollection* isoTracks, GlobalPoint vertex) { + int electronTrackCount(const reco::Track* const tr, + const reco::TrackCollection* isoTracks, + GlobalPoint vertex) const { return electronIsolation(tr, isoTracks, vertex).first; } @@ -101,60 +111,62 @@ class EgammaHLTTrackIsolation { /// set useVertex=false to consider all tracks for isolation int photonTrackCount(const reco::RecoCandidate* const recocand, const reco::TrackCollection* isoTracks, - bool useVertex) { + bool useVertex) const { return photonIsolation(recocand, isoTracks, useVertex).first; } int photonTrackCount(const reco::RecoCandidate* const recocand, const reco::TrackCollection* isoTracks, - GlobalPoint vertex) { + GlobalPoint vertex) const { return photonIsolation(recocand, isoTracks, vertex).first; } int photonTrackCount(const reco::RecoCandidate* const recocand, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks) { + const reco::TrackCollection* isoTracks) const { return photonIsolation(recocand, allEle, isoTracks).first; } /// Get Pt sum of tracks inside an isolation cone for electrons - float electronPtSum(const reco::Track* const tr, const reco::TrackCollection* isoTracks) { + float electronPtSum(const reco::Track* const tr, const reco::TrackCollection* isoTracks) const { return electronIsolation(tr, isoTracks).second; } - float electronPtSum(const reco::Track* const tr, const reco::TrackCollection* isoTracks, GlobalPoint vertex) { + float electronPtSum(const reco::Track* const tr, const reco::TrackCollection* isoTracks, GlobalPoint vertex) const { return electronIsolation(tr, isoTracks, vertex).second; } float electronPtSum(const reco::Track* const tr, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks) { + const reco::TrackCollection* isoTracks) const { return electronIsolation(tr, allEle, isoTracks).second; } /// Get Pt sum of tracks inside an isolation cone for photons /// set useVertex=true to use Photon vertex from EgammaPhotonVtxFinder /// set useVertex=false to consider all tracks for isolation - float photonPtSum(const reco::RecoCandidate* const recocand, const reco::TrackCollection* isoTracks, bool useVertex) { + float photonPtSum(const reco::RecoCandidate* const recocand, + const reco::TrackCollection* isoTracks, + bool useVertex) const { return photonIsolation(recocand, isoTracks, useVertex).second; } float photonPtSum(const reco::RecoCandidate* const recocand, const reco::TrackCollection* isoTracks, - GlobalPoint vertex) { + GlobalPoint vertex) const { return photonIsolation(recocand, isoTracks, vertex).second; } float photonPtSum(const reco::RecoCandidate* const recocand, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks) { + const reco::TrackCollection* isoTracks) const { return photonIsolation(recocand, allEle, isoTracks).second; } /// Get pt cut for itracks. - double getPtMin() { return ptMin; } + double getPtMin() const { return ptMin; } /// Get isolation cone size. - double getConeSize() { return conesize; } + double getConeSize() const { return conesize; } /// Get maximum ivertex z-coordinate spread. - double getZspan() { return zspan; } + double getZspan() const { return zspan; } /// Get maximum transverse distance of ivertex from beam line. - double getRspan() { return rspan; } + double getRspan() const { return rspan; } /// Get veto cone size - double getvetoConesize() { return vetoConesize; } + double getvetoConesize() const { return vetoConesize; } private: // Call track reconstruction @@ -162,11 +174,18 @@ class EgammaHLTTrackIsolation { GlobalPoint vtx, const reco::TrackCollection* isoTracks, bool isElectron, - bool useVertex = true); + bool useVertex = true) const; std::pair findIsoTracksWithoutEle(GlobalVector mom, GlobalPoint vtx, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks); + const reco::TrackCollection* isoTracks) const { + return findIsoTracksWithoutEle(mom.eta(), mom.phi(), vtx, allEle, isoTracks); + } + std::pair findIsoTracksWithoutEle(float centreEta, + float centrePhi, + GlobalPoint vtx, + const reco::ElectronCollection* allEle, + const reco::TrackCollection* isoTracks) const; // Parameters of isolation cone geometry. double ptMin; diff --git a/RecoEgamma/EgammaIsolationAlgos/src/EgammaHLTTrackIsolation.cc b/RecoEgamma/EgammaIsolationAlgos/src/EgammaHLTTrackIsolation.cc index b8fecea64561f..08957e2d7c72c 100644 --- a/RecoEgamma/EgammaIsolationAlgos/src/EgammaHLTTrackIsolation.cc +++ b/RecoEgamma/EgammaIsolationAlgos/src/EgammaHLTTrackIsolation.cc @@ -16,7 +16,7 @@ #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h" std::pair EgammaHLTTrackIsolation::electronIsolation(const reco::Track* const tr, - const reco::TrackCollection* isoTracks) { + const reco::TrackCollection* isoTracks) const { GlobalPoint vtx(0, 0, tr->vertex().z()); const reco::Track::Vector& p = tr->momentum(); GlobalVector mom(p.x(), p.y(), p.z()); @@ -25,7 +25,7 @@ std::pair EgammaHLTTrackIsolation::electronIsolation(const reco::Tra std::pair EgammaHLTTrackIsolation::electronIsolation(const reco::Track* const tr, const reco::TrackCollection* isoTracks, - GlobalPoint zvtx) { + GlobalPoint zvtx) const { // Just to insure consistency with no-vertex-code GlobalPoint vtx(0, 0, zvtx.z()); const reco::Track::Vector& p = tr->momentum(); @@ -35,7 +35,7 @@ std::pair EgammaHLTTrackIsolation::electronIsolation(const reco::Tra std::pair EgammaHLTTrackIsolation::electronIsolation(const reco::Track* const tr, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks) { + const reco::TrackCollection* isoTracks) const { GlobalPoint vtx(0, 0, tr->vertex().z()); const reco::Track::Vector& p = tr->momentum(); GlobalVector mom(p.x(), p.y(), p.z()); @@ -44,7 +44,7 @@ std::pair EgammaHLTTrackIsolation::electronIsolation(const reco::Tra std::pair EgammaHLTTrackIsolation::photonIsolation(const reco::RecoCandidate* const recocandidate, const reco::TrackCollection* isoTracks, - bool useVertex) { + bool useVertex) const { if (useVertex) { GlobalPoint vtx(0, 0, recocandidate->vertex().z()); return photonIsolation(recocandidate, isoTracks, vtx); @@ -57,11 +57,16 @@ std::pair EgammaHLTTrackIsolation::photonIsolation(const reco::RecoC std::pair EgammaHLTTrackIsolation::photonIsolation(const reco::RecoCandidate* const recocandidate, const reco::TrackCollection* isoTracks, - GlobalPoint zvtx) { + GlobalPoint zvtx) const { + return photonIsolation(recocandidate->superCluster()->position(), isoTracks, zvtx); +} + +std::pair EgammaHLTTrackIsolation::photonIsolation(const reco::RecoCandidate::Point& pos, + const reco::TrackCollection* isoTracks, + GlobalPoint zvtx) const { // to insure consistency with no-free-vertex-code GlobalPoint vtx(0, 0, zvtx.z()); - reco::RecoCandidate::Point pos = recocandidate->superCluster()->position(); GlobalVector mom(pos.x() - vtx.x(), pos.y() - vtx.y(), pos.z() - vtx.z()); return findIsoTracks(mom, vtx, isoTracks, false); @@ -69,14 +74,20 @@ std::pair EgammaHLTTrackIsolation::photonIsolation(const reco::RecoC std::pair EgammaHLTTrackIsolation::photonIsolation(const reco::RecoCandidate* const recocandidate, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks) { - reco::RecoCandidate::Point pos = recocandidate->superCluster()->position(); - GlobalVector mom(pos.x(), pos.y(), pos.z()); - return findIsoTracksWithoutEle(mom, GlobalPoint(), allEle, isoTracks); + const reco::TrackCollection* isoTracks) const { + const auto& sc = recocandidate->superCluster(); + return photonIsolation(sc->eta(), sc->phi(), allEle, isoTracks); +} + +std::pair EgammaHLTTrackIsolation::photonIsolation(float phoEta, + float phoPhi, + const reco::ElectronCollection* allEle, + const reco::TrackCollection* isoTracks) const { + return findIsoTracksWithoutEle(phoEta, phoPhi, GlobalPoint(), allEle, isoTracks); } std::pair EgammaHLTTrackIsolation::findIsoTracks( - GlobalVector mom, GlobalPoint vtx, const reco::TrackCollection* isoTracks, bool isElectron, bool useVertex) { + GlobalVector mom, GlobalPoint vtx, const reco::TrackCollection* isoTracks, bool isElectron, bool useVertex) const { // Check that reconstructed tracks fit within cone boundaries, // (Note: tracks will not always stay within boundaries) int ntrack = 0; @@ -113,19 +124,6 @@ std::pair EgammaHLTTrackIsolation::findIsoTracks( float R = sqrt(dphi * dphi + deta * deta); - // Apply boundary cut - // bool selected=false; - - // if (pt > ptMin && R < conesize && - // fabs(dperp) < rspan && fabs(dz) < zspan) selected=true; - - // if (selected) { - // ntrack++; - // if (!isElectron || R > vetoConesize) ptSum+=pt; //to exclude electron track - // } - // float theVetoVar = R; - // if (isElectron) theVetoVar = R; - //hmm how do I figure out if this is barrel or endcap? //abs(mom.eta())<1.5 is the obvious choice but that will be electron not detector eta for electrons //well lets leave it as that for now, its what reco does (well with eta=1.479) @@ -138,15 +136,14 @@ std::pair EgammaHLTTrackIsolation::findIsoTracks( } } - // if (isElectron) ntrack-=1; //to exclude electron track - return (std::pair(ntrack, ptSum)); } -std::pair EgammaHLTTrackIsolation::findIsoTracksWithoutEle(GlobalVector mom, +std::pair EgammaHLTTrackIsolation::findIsoTracksWithoutEle(float centerEta, + float centerPhi, GlobalPoint vtx, const reco::ElectronCollection* allEle, - const reco::TrackCollection* isoTracks) { + const reco::TrackCollection* isoTracks) const { // Check that reconstructed tracks fit within cone boundaries, // (Note: tracks will not always stay within boundaries) int ntrack = 0; @@ -171,8 +168,8 @@ std::pair EgammaHLTTrackIsolation::findIsoTracksWithoutEle(GlobalVec float pt = imom.perp(); float dperp = ivtx.perp() - vtx.perp(); float dz = ivtx.z() - vtx.z(); - float deta = imom.eta() - mom.eta(); - float dphi = imom.phi() - mom.phi(); + float deta = imom.eta() - centerEta; + float dphi = imom.phi() - centerPhi; // Correct dmom_phi's from [-2pi,2pi] to [-pi,pi] if (dphi > M_PI) @@ -187,9 +184,9 @@ std::pair EgammaHLTTrackIsolation::findIsoTracksWithoutEle(GlobalVec bool passedconeveto = true; //hmm how do I figure out if this is barrel or endcap? - //abs(mom.eta())<1.5 is the obvious choice but that will be electron not detector eta for electrons + //abs(centreEta)<1.5 is the obvious choice but that will be electron not detector eta for electrons //well lets leave it as that for now, its what reco does (well with eta=1.479) - double innerStrip = fabs(mom.eta()) < 1.479 ? stripBarrel : stripEndcap; + double innerStrip = fabs(centerEta) < 1.479 ? stripBarrel : stripEndcap; if (pt > ptMin && R < conesize && fabs(dperp) < rspan && fabs(dz) < zspan && fabs(deta) >= innerStrip) selected = true; diff --git a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py index f5b5a03fbfe40..95c23b009f922 100644 --- a/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py +++ b/RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py @@ -46,7 +46,19 @@ gsfElectronCalibEnergyErrSource = cms.InputTag(""), gsfElectronCalibEcalEnergySource = cms.InputTag(""), gsfElectronCalibEcalEnergyErrSource = cms.InputTag(""), - hcalHitSel = interestingEgammaIsoHCALSel + hcalHitSel = interestingEgammaIsoHCALSel, + trksForSCIso = cms.InputTag("generalTracks"), + superClustersToSlim = cms.VInputTag("particleFlowSuperClusterECAL:particleFlowSuperClusterECALBarrel", + "particleFlowSuperClusterECAL:particleFlowSuperClusterECALEndcapWithPreshower"), + scTrkIsol = cms.PSet( + ptMin = cms.double(0.5), + coneSize = cms.double(0.4), + zspan = cms.double(99999.), + rspan = cms.double(99999.), + vetoConeSize = cms.double(0.06), + stripBarrel = cms.double(0.03), + stripEndcap = cms.double(0.03), + ) ) from Configuration.Eras.Modifier_phase2_common_cff import phase2_common diff --git a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc index 56a87fd6abbbd..f1381452992f3 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/ReducedEGProducer.cc @@ -24,6 +24,8 @@ #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h" #include "DataFormats/EgammaReco/interface/ClusterShape.h" #include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/EgammaReco/interface/SlimmedSuperCluster.h" +#include "DataFormats/EgammaReco/interface/SlimmedSuperClusterFwd.h" #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" #include "FWCore/Framework/interface/ESHandle.h" @@ -37,6 +39,7 @@ #include "Geometry/CaloTopology/interface/CaloTopology.h" #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" #include "RecoEgamma/EgammaIsolationAlgos/interface/EGHcalRecHitSelector.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHLTTrackIsolation.h" #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h" #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h" @@ -147,6 +150,8 @@ class ReducedEGProducer : public edm::stream::EDProducer<> { const edm::ValueMap& ecalEnergyMap, const edm::ValueMap& ecalEnergyErrMap); + std::unique_ptr produceSlimmedSuperClusters(const edm::Event& event) const; + template void setToken(edm::EDGetTokenT& token, const edm::ParameterSet& config, const std::string& name) { token = consumes(config.getParameter(name)); @@ -189,6 +194,9 @@ class ReducedEGProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT> gsfElectronCalibEcalEnergyT_; edm::EDGetTokenT> gsfElectronCalibEcalEnergyErrT_; + std::vector> superClustersToSaveTs_; + const edm::EDGetTokenT trksForSCIso_; + edm::ESGetToken caloTopology_; //names for output collections const edm::EDPutTokenT outPhotons_; @@ -217,6 +225,7 @@ class ReducedEGProducer : public edm::stream::EDProducer<> { const std::vector>> outPhotonFloatValueMaps_; std::vector>> outOOTPhotonFloatValueMaps_; const std::vector>> outGsfElectronFloatValueMaps_; + const edm::EDPutTokenT outSlimmedSuperClusters_; const StringCutObjectSelector keepPhotonSel_; const StringCutObjectSelector slimRelinkPhotonSel_; @@ -229,6 +238,7 @@ class ReducedEGProducer : public edm::stream::EDProducer<> { const StringCutObjectSelector relinkGsfElectronSel_; EGHcalRecHitSelector hcalHitSel_; + EgammaHLTTrackIsolation trkIsoCalc_; }; #include "FWCore/Framework/interface/MakerMacros.h" @@ -285,6 +295,7 @@ ReducedEGProducer::ReducedEGProducer(const edm::ParameterSet& config) applyPhotonCalibOnMC_(config.getParameter("applyPhotonCalibOnMC")), applyGsfElectronCalibOnData_(config.getParameter("applyGsfElectronCalibOnData")), applyGsfElectronCalibOnMC_(config.getParameter("applyGsfElectronCalibOnMC")), + trksForSCIso_(consumes(config.getParameter("trksForSCIso"))), //output collections outPhotons_{produces("reducedGedPhotons")}, outPhotonCores_{produces("reducedGedPhotonCores")}, @@ -309,6 +320,7 @@ ReducedEGProducer::ReducedEGProducer(const edm::ParameterSet& config) vproduces>(config.getParameter>("photonFloatValueMapOutput"))}, outGsfElectronFloatValueMaps_{vproduces>( config.getParameter>("gsfElectronFloatValueMapOutput"))}, + outSlimmedSuperClusters_{produces("slimmedSuperClusters")}, keepPhotonSel_(config.getParameter("keepPhotons")), slimRelinkPhotonSel_(config.getParameter("slimRelinkPhotons")), relinkPhotonSel_(config.getParameter("relinkPhotons")), @@ -318,7 +330,8 @@ ReducedEGProducer::ReducedEGProducer(const edm::ParameterSet& config) keepGsfElectronSel_(config.getParameter("keepGsfElectrons")), slimRelinkGsfElectronSel_(config.getParameter("slimRelinkGsfElectrons")), relinkGsfElectronSel_(config.getParameter("relinkGsfElectrons")), - hcalHitSel_(config.getParameter("hcalHitSel"), consumesCollector()) { + hcalHitSel_(config.getParameter("hcalHitSel"), consumesCollector()), + trkIsoCalc_(config.getParameter("scTrkIsol")) { const auto& aTag = config.getParameter("ootPhotons"); caloTopology_ = esConsumes(); if (not aTag.label().empty()) @@ -355,6 +368,10 @@ ReducedEGProducer::ReducedEGProducer(const edm::ParameterSet& config) setToken(gsfElectronCalibEcalEnergyErrT_, config, "gsfElectronCalibEcalEnergyErrSource"); } + for (const auto& tag : config.getParameter>("superClustersToSlim")) { + superClustersToSaveTs_.emplace_back(consumes(tag)); + } + if (!ootPhotonT_.isUninitialized()) { outOOTPhotons_ = produces("reducedOOTPhotons"); outOOTPhotonCores_ = produces("reducedOOTPhotonCores"); @@ -447,6 +464,7 @@ void ReducedEGProducer::produce(edm::Event& event, const edm::EventSetup& eventS EcalRecHitCollection eeRecHits; EcalRecHitCollection esRecHits; HBHERecHitCollection hbheRecHits; + reco::SlimmedSuperClusterCollection slimmedMustacheSuperClusters; edm::ValueMap> photonPfCandMap; edm::ValueMap> gsfElectronPfCandMap; @@ -944,6 +962,8 @@ void ReducedEGProducer::produce(edm::Event& event, const edm::EventSetup& eventS for (auto const& vals : gsfElectronFloatValueMapVals) { emplaceValueMap(outGsfElectronHandle, vals, event, outGsfElectronFloatValueMaps_[index++]); } + + event.put(outSlimmedSuperClusters_, std::move(produceSlimmedSuperClusters(event))); } template @@ -1222,3 +1242,23 @@ void ReducedEGProducer::calibrateElectron(reco::GsfElectron& electron, math::XYZTLorentzVector newP4{oldP4.x() * corr, oldP4.y() * corr, oldP4.z() * corr, newEnergy}; electron.correctMomentum(newP4, electron.trackMomentumError(), newEnergyErr); } + +std::unique_ptr ReducedEGProducer::produceSlimmedSuperClusters( + const edm::Event& event) const { + auto outputSCs = std::make_unique(); + for (const auto& inputSCToken : superClustersToSaveTs_) { + auto inputSCs = event.get(inputSCToken); + for (const auto& inputSC : inputSCs) { + outputSCs->emplace_back(inputSC); + } + } + + auto tracks = event.get(trksForSCIso_); + + for (auto& sc : *outputSCs) { + float trkIso = trkIsoCalc_.photonIsolation(sc.position(), &tracks).second; + sc.setTrkIso(trkIso); + } + + return outputSCs; +}