diff --git a/DataFormats/L1TrackTrigger/interface/L1Track.h b/DataFormats/L1TrackTrigger/interface/L1Track.h new file mode 100644 index 0000000000000..2cc2d83c79fec --- /dev/null +++ b/DataFormats/L1TrackTrigger/interface/L1Track.h @@ -0,0 +1,11 @@ +#ifndef DataFormats_L1TrackTrigger_L1Track_h +#define DataFormats_L1TrackTrigger_L1Track_h + +#include "DataFormats/L1TrackTrigger/interface/TTTrack.h" +#include "DataFormats/L1TrackTrigger/interface/TTTypes.h" +#include + +using L1Track = TTTrack; +using L1TrackCollection = std::vector; + +#endif diff --git a/RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTEleL1TrackIsolProducer.cc b/RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTEleL1TrackIsolProducer.cc new file mode 100644 index 0000000000000..611ea49c506a5 --- /dev/null +++ b/RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTEleL1TrackIsolProducer.cc @@ -0,0 +1,93 @@ +// Author: Sam Harper (RAL/CERN) +// L1 track isolation producer for the HLT +// A quick primer on how the E/gamma HLT works w.r.t to ID variables +// 1. the supercluster is the primary object +// 2. superclusters get id variables associated to them via association maps keyed +// to the supercluster +// However here we also need to read in electron objects as we need to solve for which +// GsfTrack associated to the supercluster to use for the isolation +// The electron producer solves for this and assigns the electron the best GsfTrack +// which we will use for the vz of the electron +// One thing which Swagata Mukherjee pointed out is that we have to be careful of +// is getting a bad GsfTrack with a bad vertex which will give us a fake vz which then +// leads to a perfectly isolated electron as that random vz is not a vertex + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/EgammaCandidates/interface/Electron.h" +#include "DataFormats/EgammaCandidates/interface/ElectronFwd.h" +#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h" +#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h" +#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/TrackReco/interface/Track.h" + +#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaL1TkIsolation.h" + +#include + +class EgammaHLTEleL1TrackIsolProducer : public edm::global::EDProducer<> { +public: + explicit EgammaHLTEleL1TrackIsolProducer(const edm::ParameterSet&); + ~EgammaHLTEleL1TrackIsolProducer() override = default; + void produce(edm::StreamID sid, edm::Event&, const edm::EventSetup&) const override; + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + const edm::EDGetTokenT ecalCandsToken_; + const edm::EDGetTokenT elesToken_; + const edm::EDGetTokenT l1TrksToken_; + EgammaL1TkIsolation isolAlgo_; +}; + +EgammaHLTEleL1TrackIsolProducer::EgammaHLTEleL1TrackIsolProducer(const edm::ParameterSet& config) + : ecalCandsToken_(consumes(config.getParameter("ecalCands"))), + elesToken_(consumes(config.getParameter("eles"))), + l1TrksToken_(consumes(config.getParameter("l1Tracks"))), + isolAlgo_(config.getParameter("isolCfg")) { + produces(); +} + +void EgammaHLTEleL1TrackIsolProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("ecalCands", edm::InputTag("hltEgammaCandidates")); + desc.add("eles", edm::InputTag("hltEgammaGsfElectrons")); + desc.add("l1Tracks", edm::InputTag("TTTracksFromTrackletEmulation", "Level1TTTracks")); + desc.add("isolCfg", EgammaL1TkIsolation::makePSetDescription()); + descriptions.add("hltEgammaHLTEleL1TrackIsolProducer", desc); +} +void EgammaHLTEleL1TrackIsolProducer::produce(edm::StreamID sid, + edm::Event& iEvent, + const edm::EventSetup& iSetup) const { + auto ecalCands = iEvent.getHandle(ecalCandsToken_); + auto eles = iEvent.getHandle(elesToken_); + auto l1Trks = iEvent.getHandle(l1TrksToken_); + + auto recoEcalCandMap = std::make_unique(ecalCands); + + for (size_t candNr = 0; candNr < ecalCands->size(); candNr++) { + reco::RecoEcalCandidateRef recoEcalCandRef(ecalCands, candNr); + reco::ElectronRef eleRef; + for (size_t eleNr = 0; eleNr < eles->size(); eleNr++) { + if ((*eles)[eleNr].superCluster() == recoEcalCandRef->superCluster()) { + eleRef = reco::ElectronRef(eles, eleNr); + break; + } + } + + float isol = + eleRef.isNonnull() ? isolAlgo_.calIsol(*eleRef->gsfTrack(), *l1Trks).second : std::numeric_limits::max(); + + recoEcalCandMap->insert(recoEcalCandRef, isol); + } + iEvent.put(std::move(recoEcalCandMap)); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(EgammaHLTEleL1TrackIsolProducer); diff --git a/RecoEgamma/EgammaIsolationAlgos/interface/EgammaL1TkIsolation.h b/RecoEgamma/EgammaIsolationAlgos/interface/EgammaL1TkIsolation.h new file mode 100644 index 0000000000000..4ab686e5b9d2f --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/interface/EgammaL1TkIsolation.h @@ -0,0 +1,61 @@ +#ifndef RecoEgamma_EgammaIsolationAlgos_EgammaL1TkIsolation_h +#define RecoEgamma_EgammaIsolationAlgos_EgammaL1TkIsolation_h + +#include "DataFormats/L1TrackTrigger/interface/L1Track.h" +#include "DataFormats/TrackReco/interface/TrackBase.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +//author S. Harper (RAL/CERN) +//based on the work of Swagata Mukherjee and Giulia Sorrentino + +class EgammaL1TkIsolation { +public: + explicit EgammaL1TkIsolation(const edm::ParameterSet& para); + + static void fillPSetDescription(edm::ParameterSetDescription& desc); + static edm::ParameterSetDescription makePSetDescription() { + edm::ParameterSetDescription desc; + fillPSetDescription(desc); + return desc; + } + + std::pair calIsol(const reco::TrackBase& trk, const L1TrackCollection& l1Tks) const; + + std::pair calIsol(const double objEta, + const double objPhi, + const double objZ, + const L1TrackCollection& l1Tks) const; + + //little helper function for the two calIsol functions for it to directly return the pt + template + double calIsolPt(Args&&... args) const { + return calIsol(std::forward(args)...).second; + } + +private: + struct TrkCuts { + float minPt; + float minDR2; + float maxDR2; + float minDEta; + float maxDZ; + explicit TrkCuts(const edm::ParameterSet& para); + static edm::ParameterSetDescription makePSetDescription(); + }; + + size_t etaBinNr(double eta) const; + static bool passTrkSel(const L1Track& trk, + const double trkPt, + const TrkCuts& cuts, + const double objEta, + const double objPhi, + const double objZ); + + bool useAbsEta_; + std::vector etaBoundaries_; + std::vector trkCuts_; +}; + +#endif diff --git a/RecoEgamma/EgammaIsolationAlgos/src/EgammaL1TkIsolation.cc b/RecoEgamma/EgammaIsolationAlgos/src/EgammaL1TkIsolation.cc new file mode 100644 index 0000000000000..1fa0aef29b150 --- /dev/null +++ b/RecoEgamma/EgammaIsolationAlgos/src/EgammaL1TkIsolation.cc @@ -0,0 +1,99 @@ +#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaL1TkIsolation.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/Math/interface/deltaR.h" + +#include + +EgammaL1TkIsolation::EgammaL1TkIsolation(const edm::ParameterSet& para) + : useAbsEta_(para.getParameter("useAbsEta")), + etaBoundaries_(para.getParameter>("etaBoundaries")) { + const auto& trkCutParams = para.getParameter>("trkCuts"); + for (const auto& params : trkCutParams) { + trkCuts_.emplace_back(TrkCuts(params)); + } + if (etaBoundaries_.size() + 1 != trkCuts_.size()) { + throw cms::Exception("ConfigError") << "EgammaL1TkIsolation: etaBoundaries parameters size (" + << etaBoundaries_.size() + << ") should be one less than the size of trkCuts VPSet (" << trkCuts_.size() + << ")"; + } + if (!std::is_sorted(etaBoundaries_.begin(), etaBoundaries_.end())) { + throw cms::Exception("ConfigError") + << "EgammaL1TkIsolation: etaBoundaries parameter's entries should be in increasing value"; + } +} + +void EgammaL1TkIsolation::fillPSetDescription(edm::ParameterSetDescription& desc) { + desc.add("useAbsEta", true); + desc.add("etaBoundaries", std::vector{1.5}); + desc.addVPSet("trkCuts", TrkCuts::makePSetDescription(), {edm::ParameterSet(), edm::ParameterSet()}); +} + +std::pair EgammaL1TkIsolation::calIsol(const reco::TrackBase& trk, const L1TrackCollection& tracks) const { + return calIsol(trk.eta(), trk.phi(), trk.vz(), tracks); +} + +std::pair EgammaL1TkIsolation::calIsol(const double objEta, + const double objPhi, + const double objZ, + const L1TrackCollection& tracks) const { + double ptSum = 0.; + int nrTrks = 0; + + const TrkCuts& cuts = trkCuts_[etaBinNr(objEta)]; + + for (const auto& trk : tracks) { + const float trkPt = trk.momentum().perp(); + if (passTrkSel(trk, trkPt, cuts, objEta, objPhi, objZ)) { + ptSum += trkPt; + nrTrks++; + } + } + return {nrTrks, ptSum}; +} + +EgammaL1TkIsolation::TrkCuts::TrkCuts(const edm::ParameterSet& para) { + minPt = para.getParameter("minPt"); + auto sq = [](double val) { return val * val; }; + minDR2 = sq(para.getParameter("minDR")); + maxDR2 = sq(para.getParameter("maxDR")); + minDEta = para.getParameter("minDEta"); + maxDZ = para.getParameter("maxDZ"); +} + +edm::ParameterSetDescription EgammaL1TkIsolation::TrkCuts::makePSetDescription() { + edm::ParameterSetDescription desc; + desc.add("minPt", 2.0); + desc.add("maxDR", 0.3); + desc.add("minDR", 0.01); + desc.add("minDEta", 0.003); + desc.add("maxDZ", 0.7); + return desc; +} + +//as we have verfied that trkCuts_ size is etaBoundaries_ size +1 +//then this is always a valid binnr for trkCuts_ +size_t EgammaL1TkIsolation::etaBinNr(double eta) const { + if (useAbsEta_) { + eta = std::abs(eta); + } + auto res = std::upper_bound(etaBoundaries_.begin(), etaBoundaries_.end(), eta); + size_t binNr = std::distance(etaBoundaries_.begin(), res); + return binNr; +} + +bool EgammaL1TkIsolation::passTrkSel(const L1Track& trk, + const double trkPt, + const TrkCuts& cuts, + const double objEta, + const double objPhi, + const double objZ) { + if (trkPt > cuts.minPt && std::abs(objZ - trk.z0()) < cuts.maxDZ) { + const float trkEta = trk.eta(); + const float dEta = trkEta - objEta; + const float dR2 = reco::deltaR2(objEta, objPhi, trkEta, trk.phi()); + return dR2 >= cuts.minDR2 && dR2 <= cuts.maxDR2 && std::abs(dEta) >= cuts.minDEta; + } + + return false; +}