E/gamma HLT L1 Track Isolation Producer & TTTrack typedef: 11_1_X

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 <vector>
+
+using L1Track = TTTrack<Ref_Phase2TrackerDigi_>;
+using L1TrackCollection = std::vector<L1Track>;
+
+#endif

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 <memory>
+
+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<reco::RecoEcalCandidateCollection> ecalCandsToken_;
+  const edm::EDGetTokenT<reco::ElectronCollection> elesToken_;
+  const edm::EDGetTokenT<L1TrackCollection> l1TrksToken_;
+  EgammaL1TkIsolation isolAlgo_;
+};
+
+EgammaHLTEleL1TrackIsolProducer::EgammaHLTEleL1TrackIsolProducer(const edm::ParameterSet& config)
+    : ecalCandsToken_(consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("ecalCands"))),
+      elesToken_(consumes<reco::ElectronCollection>(config.getParameter<edm::InputTag>("eles"))),
+      l1TrksToken_(consumes<L1TrackCollection>(config.getParameter<edm::InputTag>("l1Tracks"))),
+      isolAlgo_(config.getParameter<edm::ParameterSet>("isolCfg")) {
+  produces<reco::RecoEcalCandidateIsolationMap>();
+}
+
+void EgammaHLTEleL1TrackIsolProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("ecalCands", edm::InputTag("hltEgammaCandidates"));
+  desc.add<edm::InputTag>("eles", edm::InputTag("hltEgammaGsfElectrons"));
+  desc.add<edm::InputTag>("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<reco::RecoEcalCandidateIsolationMap>(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<float>::max();
+
+    recoEcalCandMap->insert(recoEcalCandRef, isol);
+  }
+  iEvent.put(std::move(recoEcalCandMap));
+}
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(EgammaHLTEleL1TrackIsolProducer);

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<int, double> calIsol(const reco::TrackBase& trk, const L1TrackCollection& l1Tks) const;
+
+  std::pair<int, double> 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 <typename... Args>
+  double calIsolPt(Args&&... args) const {
+    return calIsol(std::forward<Args>(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<double> etaBoundaries_;
+  std::vector<TrkCuts> trkCuts_;
+};
+
+#endif

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 <algorithm>
+
+EgammaL1TkIsolation::EgammaL1TkIsolation(const edm::ParameterSet& para)
+    : useAbsEta_(para.getParameter<bool>("useAbsEta")),
+      etaBoundaries_(para.getParameter<std::vector<double>>("etaBoundaries")) {
+  const auto& trkCutParams = para.getParameter<std::vector<edm::ParameterSet>>("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<int, double> EgammaL1TkIsolation::calIsol(const reco::TrackBase& trk, const L1TrackCollection& tracks) const {
+  return calIsol(trk.eta(), trk.phi(), trk.vz(), tracks);
+}
+
+std::pair<int, double> 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<double>("minPt");
+  auto sq = [](double val) { return val * val; };
+  minDR2 = sq(para.getParameter<double>("minDR"));
+  maxDR2 = sq(para.getParameter<double>("maxDR"));
+  minDEta = para.getParameter<double>("minDEta");
+  maxDZ = para.getParameter<double>("maxDZ");
+}
+
+edm::ParameterSetDescription EgammaL1TkIsolation::TrkCuts::makePSetDescription() {
+  edm::ParameterSetDescription desc;
+  desc.add<double>("minPt", 2.0);
+  desc.add<double>("maxDR", 0.3);
+  desc.add<double>("minDR", 0.01);
+  desc.add<double>("minDEta", 0.003);
+  desc.add<double>("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;
+}