Producer and filter for jet timing trigger

HLTrigger/JetMET/plugins/HLTJetTimingFilter.cc

@@ -0,0 +1,8 @@
+#include "HLTJetTimingFilter.h"
+
+typedef HLTJetTimingFilter<reco::CaloJet> HLTCaloJetTimingFilter;
+typedef HLTJetTimingFilter<reco::PFJet> HLTPFJetTimingFilter;
+
+// declare classes as framework plugins
+DEFINE_FWK_MODULE(HLTCaloJetTimingFilter);
+DEFINE_FWK_MODULE(HLTPFJetTimingFilter);

HLTrigger/JetMET/plugins/HLTJetTimingFilter.h

@@ -0,0 +1,123 @@
+/** \class HLTJetTimingFilter
+ *
+ *  \brief  This makes selections on the timing and associated ecal cells 
+ *  produced by HLTJetTimingProducer
+ *  \author Matthew Citron
+ *
+ */
+#ifndef HLTrigger_JetMET_plugins_HLTJetTimingFilter_h
+#define HLTrigger_JetMET_plugins_HLTJetTimingFilter_h
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "HLTrigger/HLTcore/interface/HLTFilter.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
+
+namespace edm {
+  class ConfigurationDescriptions;
+}
+
+//
+// class declaration
+//
+template <typename T>
+class HLTJetTimingFilter : public HLTFilter {
+public:
+  explicit HLTJetTimingFilter(const edm::ParameterSet& iConfig);
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+  bool hltFilter(edm::Event&,
+                 const edm::EventSetup&,
+                 trigger::TriggerFilterObjectWithRefs& filterproduct) const override;
+
+private:
+  // Input collections
+  const edm::InputTag jetInput_;
+  const edm::EDGetTokenT<std::vector<T>> jetInputToken_;
+  const edm::EDGetTokenT<edm::ValueMap<float>> jetTimesInputToken_;
+  const edm::EDGetTokenT<edm::ValueMap<unsigned int>> jetCellsForTimingInputToken_;
+  const edm::EDGetTokenT<edm::ValueMap<float>> jetEcalEtForTimingInputToken_;
+
+  // Thresholds for selection
+  const unsigned int minJets_;
+  const double jetTimeThresh_;
+  const double jetEcalEtForTimingThresh_;
+  const unsigned int jetCellsForTimingThresh_;
+  const double minPt_;
+};
+
+//Constructor
+template <typename T>
+HLTJetTimingFilter<T>::HLTJetTimingFilter(const edm::ParameterSet& iConfig)
+    : HLTFilter(iConfig),
+      jetInput_{iConfig.getParameter<edm::InputTag>("jets")},
+      jetInputToken_{consumes<std::vector<T>>(jetInput_)},
+      jetTimesInputToken_{consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("jetTimes"))},
+      jetCellsForTimingInputToken_{
+          consumes<edm::ValueMap<unsigned int>>(iConfig.getParameter<edm::InputTag>("jetCellsForTiming"))},
+      jetEcalEtForTimingInputToken_{
+          consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("jetEcalEtForTiming"))},
+      minJets_{iConfig.getParameter<unsigned int>("minJets")},
+      jetTimeThresh_{iConfig.getParameter<double>("jetTimeThresh")},
+      jetEcalEtForTimingThresh_{iConfig.getParameter<double>("jetEcalEtForTimingThresh")},
+      jetCellsForTimingThresh_{iConfig.getParameter<unsigned int>("jetCellsForTimingThresh")},
+      minPt_{iConfig.getParameter<double>("minJetPt")} {}
+
+//Filter
+template <typename T>
+bool HLTJetTimingFilter<T>::hltFilter(edm::Event& iEvent,
+                                      const edm::EventSetup& iSetup,
+                                      trigger::TriggerFilterObjectWithRefs& filterproduct) const {
+  // typedef vector<T> TCollection;
+  // typedef edm::Ref<TCollection> TRef;
+  if (saveTags())
+    filterproduct.addCollectionTag(jetInput_);
+
+  auto const jets = iEvent.getHandle(jetInputToken_);
+  auto const& jetTimes = iEvent.get(jetTimesInputToken_);
+  auto const& jetCellsForTiming = iEvent.get(jetCellsForTimingInputToken_);
+  auto const& jetEcalEtForTiming = iEvent.get(jetEcalEtForTimingInputToken_);
+
+  uint njets = 0;
+  for (auto iterJet = jets->begin(); iterJet != jets->end(); ++iterJet) {
+    edm::Ref<vector<T>> const caloJetRef(jets, std::distance(jets->begin(), iterJet));
+    if (iterJet->pt() > minPt_ and jetTimes[caloJetRef] > jetTimeThresh_ and
+        jetEcalEtForTiming[caloJetRef] > jetEcalEtForTimingThresh_ and
+        jetCellsForTiming[caloJetRef] > jetCellsForTimingThresh_) {
+      // add caloJetRef to the event
+      filterproduct.addObject(trigger::TriggerJet, caloJetRef);
+      ++njets;
+    }
+  }
+
+  return njets >= minJets_;
+}
+
+// Fill descriptions
+template <typename T>
+void HLTJetTimingFilter<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  makeHLTFilterDescription(desc);
+  desc.add<edm::InputTag>("jets", edm::InputTag("hltDisplacedHLTCaloJetCollectionProducerMidPt"));
+  desc.add<edm::InputTag>("jetTimes", edm::InputTag("hltDisplacedHLTCaloJetCollectionProducerMidPtTiming"));
+  desc.add<edm::InputTag>("jetCellsForTiming",
+                          edm::InputTag("hltDisplacedHLTCaloJetCollectionProducerMidPtTiming", "jetCellsForTiming"));
+  desc.add<edm::InputTag>("jetEcalEtForTiming",
+                          edm::InputTag("hltDisplacedHLTCaloJetCollectionProducerMidPtTiming", "jetEcalEtForTiming"));
+  desc.add<unsigned int>("minJets", 1);
+  desc.add<double>("jetTimeThresh", 1.);
+  desc.add<unsigned int>("jetCellsForTimingThresh", 5);
+  desc.add<double>("jetEcalEtForTimingThresh", 10.);
+  desc.add<double>("minJetPt", 40.);
+  descriptions.addWithDefaultLabel(desc);
+}
+
+#endif  // HLTrigger_JetMET_plugins_HLTJetTimingFilter_h

HLTrigger/JetMET/plugins/HLTJetTimingProducer.cc

@@ -0,0 +1,8 @@
+#include "HLTJetTimingProducer.h"
+
+typedef HLTJetTimingProducer<reco::CaloJet> HLTCaloJetTimingProducer;
+typedef HLTJetTimingProducer<reco::PFJet> HLTPFJetTimingProducer;
+
+// declare classes as framework plugins
+DEFINE_FWK_MODULE(HLTCaloJetTimingProducer);
+DEFINE_FWK_MODULE(HLTPFJetTimingProducer);

HLTrigger/JetMET/plugins/HLTJetTimingProducer.h

@@ -0,0 +1,186 @@
+/** \class HLTJetTimingProducer
+ *
+ *  \brief  This produces timing and associated ecal cell information for calo jets 
+ *  \author Matthew Citron
+ *
+ */
+#ifndef HLTrigger_JetMET_plugins_HLTJetTimingProducer_h
+#define HLTrigger_JetMET_plugins_HLTJetTimingProducer_h
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "DataFormats/Common/interface/ValueMap.h"
+
+#include "DataFormats/JetReco/interface/CaloJetCollection.h"
+#include "DataFormats/JetReco/interface/PFJetCollection.h"
+
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
+#include "DataFormats/Math/interface/deltaR.h"
+
+//
+// class declaration
+//
+template <typename T>
+class HLTJetTimingProducer : public edm::stream::EDProducer<> {
+public:
+  explicit HLTJetTimingProducer(const edm::ParameterSet&);
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+  void jetTimeFromEcalCells(const T&,
+                            const edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>&,
+                            const edm::ESHandle<CaloGeometry>&,
+                            float&,
+                            float&,
+                            uint&);
+
+  const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
+  // Input collections
+  const edm::EDGetTokenT<std::vector<T>> jetInputToken_;
+  const edm::EDGetTokenT<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEBToken_;
+  const edm::EDGetTokenT<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEEToken_;
+
+  // Include barrel, endcap jets or both
+  const bool barrelJets_;
+  const bool endcapJets_;
+
+  // Configurables for timing definition
+  const double ecalCellEnergyThresh_;
+  const double ecalCellTimeThresh_;
+  const double ecalCellTimeErrorThresh_;
+  const double matchingRadius2_;
+};
+
+//Constructor
+template <typename T>
+HLTJetTimingProducer<T>::HLTJetTimingProducer(const edm::ParameterSet& iConfig)
+    : caloGeometryToken_(esConsumes()),
+      jetInputToken_{consumes<std::vector<T>>(iConfig.getParameter<edm::InputTag>("jets"))},
+      ecalRecHitsEBToken_{consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>>(
+          iConfig.getParameter<edm::InputTag>("ebRecHitsColl"))},
+      ecalRecHitsEEToken_{consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>>(
+          iConfig.getParameter<edm::InputTag>("eeRecHitsColl"))},
+      barrelJets_{iConfig.getParameter<bool>("barrelJets")},
+      endcapJets_{iConfig.getParameter<bool>("endcapJets")},
+      ecalCellEnergyThresh_{iConfig.getParameter<double>("ecalCellEnergyThresh")},
+      ecalCellTimeThresh_{iConfig.getParameter<double>("ecalCellTimeThresh")},
+      ecalCellTimeErrorThresh_{iConfig.getParameter<double>("ecalCellTimeErrorThresh")},
+      matchingRadius2_{std::pow(iConfig.getParameter<double>("matchingRadius"), 2)} {
+  produces<edm::ValueMap<float>>("");
+  produces<edm::ValueMap<unsigned int>>("jetCellsForTiming");
+  produces<edm::ValueMap<float>>("jetEcalEtForTiming");
+}
+
+//calculate jet time
+template <typename T>
+void HLTJetTimingProducer<T>::jetTimeFromEcalCells(
+    const T& jet,
+    const edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>& ecalRecHits,
+    const edm::ESHandle<CaloGeometry>& caloGeometry,
+    float& weightedTimeCell,
+    float& totalEmEnergyCell,
+    uint& nCells) {
+  for (auto const& ecalRH : ecalRecHits) {
+    if (ecalRH.checkFlag(EcalRecHit::kSaturated) || ecalRH.checkFlag(EcalRecHit::kLeadingEdgeRecovered) ||
+        ecalRH.checkFlag(EcalRecHit::kPoorReco) || ecalRH.checkFlag(EcalRecHit::kWeird) ||
+        ecalRH.checkFlag(EcalRecHit::kDiWeird))
+      continue;
+    if (ecalRH.energy() < ecalCellEnergyThresh_)
+      continue;
+    if (ecalRH.timeError() <= 0. || ecalRH.timeError() > ecalCellTimeErrorThresh_)
+      continue;
+    if (fabs(ecalRH.time()) > ecalCellTimeThresh_)
+      continue;
+    auto const pos = caloGeometry->getPosition(ecalRH.detid());
+    if (reco::deltaR2(jet, pos) > matchingRadius2_)
+      continue;
+    weightedTimeCell += ecalRH.time() * ecalRH.energy() * sin(pos.theta());
+    totalEmEnergyCell += ecalRH.energy() * sin(pos.theta());
+    nCells++;
+  }
+  if (totalEmEnergyCell > 0) {
+    weightedTimeCell /= totalEmEnergyCell;
+  }
+}
+
+//Producer
+template <typename T>
+void HLTJetTimingProducer<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  edm::ESHandle<CaloGeometry> caloGeometry = iSetup.getHandle(caloGeometryToken_);
+  auto const jets = iEvent.getHandle(jetInputToken_);
+  auto const& ecalRecHitsEB = iEvent.get(ecalRecHitsEBToken_);
+  auto const& ecalRecHitsEE = iEvent.get(ecalRecHitsEEToken_);
+
+  std::vector<float> jetTimings;
+  std::vector<unsigned int> jetCellsForTiming;
+  std::vector<float> jetEcalEtForTiming;
+
+  jetTimings.reserve(jets->size());
+  jetEcalEtForTiming.reserve(jets->size());
+  jetCellsForTiming.reserve(jets->size());
+
+  for (auto const& jet : *jets) {
+    float weightedTimeCell = 0;
+    float totalEmEnergyCell = 0;
+    unsigned int nCells = 0;
+    if (barrelJets_)
+      jetTimeFromEcalCells(jet, ecalRecHitsEB, caloGeometry, weightedTimeCell, totalEmEnergyCell, nCells);
+    if (endcapJets_) {
+      weightedTimeCell *= totalEmEnergyCell;
+      jetTimeFromEcalCells(jet, ecalRecHitsEE, caloGeometry, weightedTimeCell, totalEmEnergyCell, nCells);
+    }
+
+    // If there is at least one ecal cell passing selection, calculate timing
+    jetTimings.emplace_back(totalEmEnergyCell > 0 ? weightedTimeCell : -50);
+    jetEcalEtForTiming.emplace_back(totalEmEnergyCell);
+    jetCellsForTiming.emplace_back(nCells);
+  }
+
+  std::unique_ptr<edm::ValueMap<float>> jetTimings_out(new edm::ValueMap<float>());
+  edm::ValueMap<float>::Filler jetTimings_filler(*jetTimings_out);
+  jetTimings_filler.insert(jets, jetTimings.begin(), jetTimings.end());
+  jetTimings_filler.fill();
+  iEvent.put(std::move(jetTimings_out), "");
+
+  std::unique_ptr<edm::ValueMap<float>> jetEcalEtForTiming_out(new edm::ValueMap<float>());
+  edm::ValueMap<float>::Filler jetEcalEtForTiming_filler(*jetEcalEtForTiming_out);
+  jetEcalEtForTiming_filler.insert(jets, jetEcalEtForTiming.begin(), jetEcalEtForTiming.end());
+  jetEcalEtForTiming_filler.fill();
+  iEvent.put(std::move(jetEcalEtForTiming_out), "jetEcalEtForTiming");
+
+  std::unique_ptr<edm::ValueMap<unsigned int>> jetCellsForTiming_out(new edm::ValueMap<unsigned int>());
+  edm::ValueMap<unsigned int>::Filler jetCellsForTiming_filler(*jetCellsForTiming_out);
+  jetCellsForTiming_filler.insert(jets, jetCellsForTiming.begin(), jetCellsForTiming.end());
+  jetCellsForTiming_filler.fill();
+  iEvent.put(std::move(jetCellsForTiming_out), "jetCellsForTiming");
+}
+
+// Fill descriptions
+template <typename T>
+void HLTJetTimingProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("jets", edm::InputTag(""));
+  desc.add<bool>("barrelJets", false);
+  desc.add<bool>("endcapJets", false);
+  desc.add<double>("ecalCellEnergyThresh", 0.5);
+  desc.add<double>("ecalCellTimeThresh", 12.5);
+  desc.add<double>("ecalCellTimeErrorThresh", 100.);
+  desc.add<double>("matchingRadius", 0.4);
+  desc.add<edm::InputTag>("ebRecHitsColl", edm::InputTag("hltEcalRecHit", "EcalRecHitsEB"));
+  desc.add<edm::InputTag>("eeRecHitsColl", edm::InputTag("hltEcalRecHit", "EcalRecHitsEE"));
+  descriptions.addWithDefaultLabel(desc);
+}
+
+#endif  // HLTrigger_JetMET_plugins_HLTJetTimingProducer_h

HLTrigger/JetMET/test/hltCaloJetTimingFilter_cfg.py

@@ -0,0 +1,53 @@
+import FWCore.ParameterSet.Config as cms
+
+process = cms.Process("Demo")
+
+process.load("FWCore.MessageService.MessageLogger_cfi")
+process.load("Geometry.CMSCommonData.cmsIdealGeometryXML_cfi");
+process.load("Geometry.CaloEventSetup.CaloGeometry_cfi");
+process.load("Geometry.CaloEventSetup.CaloTopology_cfi");
+process.load("TrackingTools/TransientTrack/TransientTrackBuilder_cfi")
+process.load("Configuration.Geometry.GeometryIdeal_cff")
+process.load("Configuration.StandardSequences.MagneticField_cff")
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+
+process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) )
+
+process.source = cms.Source("PoolSource",
+    fileNames = cms.untracked.vstring(
+        '/store/mc/Run3Winter21DRMiniAOD/HTo2LongLivedTo4b_MH-250_MFF-60_CTau-1000mm_TuneCP5_14TeV-pythia8/GEN-SIM-RECO/FlatPU30to80FEVT_112X_mcRun3_2021_realistic_v16-v2/130000/058470ca-8aaa-4727-80d8-f42621bafd39.root'
+    )
+)
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag.globaltag = '121X_mcRun3_2021_realistic_v1'
+
+process.hltTimingProducer = cms.EDProducer('HLTCaloJetTimingProducer',
+        jets = cms.InputTag( "ak4CaloJets" ),
+        ebRecHitsColl = cms.InputTag( 'ecalRecHit','EcalRecHitsEB' ),
+        eeRecHitsColl = cms.InputTag( 'ecalRecHit','EcalRecHitsEE' ),
+        barrelJets = cms.bool(True),
+        endcapJets = cms.bool(False),
+        ecalCellEnergyThresh =cms.double(0.5),
+        ecalCellTimeThresh = cms.double(12.5),
+        ecalCellTimeErrorThresh = cms.double(100.),
+        matchingRadius = cms.double(0.4),
+)
+
+process.hltTimingFilter = cms.EDFilter('HLTCaloJetTimingFilter',
+    saveTags = cms.bool( True ),
+    minJets = cms.uint32(1),
+    minJetPt = cms.double(40.0),
+    jetTimeThresh = cms.double(1.),
+    jetCellsForTimingThresh = cms.uint32(5),
+    jetEcalEtForTimingThresh = cms.double(10.),
+    jets = cms.InputTag( "ak4CaloJets" ),
+    jetTimes = cms.InputTag( "hltTimingProducer" ),
+    jetEcalEtForTiming = cms.InputTag( "hltTimingProducer" ,"jetEcalEtForTiming"),
+    jetCellsForTiming = cms.InputTag( "hltTimingProducer" ,"jetCellsForTiming"),
+)
+process.output = cms.OutputModule( "PoolOutputModule",
+    fileName = cms.untracked.string( "timingOutput.root" ),
+)
+
+process.p = cms.Path(process.hltTimingProducer+process.hltTimingFilter)
+process.Output = cms.EndPath(process.output)