Producer and filter for jet timing trigger

HLTrigger/JetMET/plugins/HLTCaloJetTimingFilter.cc

@@ -0,0 +1,122 @@
+/** \class HLTCaloJetTimingFilter
+ *
+ *  \brief  This makes selections on the timing and associated ecal cells 
+ *  produced by HLTCaloJetTimingProducer
+ *  \author Matthew Citron
+ *
+ *
+ */
+
+// 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 "FWCore/Utilities/interface/StreamID.h"
+
+#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
+
+namespace edm {
+  class ConfigurationDescriptions;
+}
+
+//
+// class declaration
+//
+class HLTCaloJetTimingFilter : public HLTFilter {
+public:
+  explicit HLTCaloJetTimingFilter(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
+  edm::InputTag jetLabel_;
+  edm::InputTag jetTimeLabel_;
+  edm::InputTag jetCellsForTimingLabel_;
+  edm::InputTag jetEcalEtForTimingLabel_;
+  //Thresholds for selection
+  unsigned int minJets_;
+  double jetTimeThresh_;
+  double jetEcalEtForTimingThresh_;
+  unsigned int jetCellsForTimingThresh_;
+  double minPt_;
+
+  edm::EDGetTokenT<reco::CaloJetCollection> jetInputToken;
+  edm::EDGetTokenT<edm::ValueMap<float>> jetTimesInputToken;
+  edm::EDGetTokenT<edm::ValueMap<unsigned int>> jetCellsForTimingInputToken;
+  edm::EDGetTokenT<edm::ValueMap<float>> jetEcalEtForTimingInputToken;
+};
+
+//Constructor
+HLTCaloJetTimingFilter::HLTCaloJetTimingFilter(const edm::ParameterSet& iConfig) : HLTFilter(iConfig) {
+  jetLabel_ = iConfig.getParameter<edm::InputTag>("jets");
+  jetTimeLabel_ = iConfig.getParameter<edm::InputTag>("jetTimes");
+  jetCellsForTimingLabel_ = iConfig.getParameter<edm::InputTag>("jetCellsForTiming");
+  jetEcalEtForTimingLabel_ = iConfig.getParameter<edm::InputTag>("jetEcalEtForTiming");
+  minJets_ = iConfig.getParameter<unsigned int>("minJets");
+  jetTimeThresh_ = iConfig.getParameter<double>("jetTimeThresh");
+  jetCellsForTimingThresh_ = iConfig.getParameter<unsigned int>("jetCellsForTimingThresh");
+  jetEcalEtForTimingThresh_ = iConfig.getParameter<double>("jetEcalEtForTimingThresh");
+  minPt_ = iConfig.getParameter<double>("minJetPt");
+  jetInputToken = consumes<std::vector<reco::CaloJet>>(jetLabel_);
+  jetTimesInputToken = consumes<edm::ValueMap<float>>(jetTimeLabel_);
+  jetCellsForTimingInputToken = consumes<edm::ValueMap<unsigned int>>(jetCellsForTimingLabel_);
+  jetEcalEtForTimingInputToken = consumes<edm::ValueMap<float>>(jetEcalEtForTimingLabel_);
+  //now do what ever initialization is needed
+}
+
+//Filter
+bool HLTCaloJetTimingFilter::hltFilter(edm::Event& iEvent,
+                                       const edm::EventSetup& iSetup,
+                                       trigger::TriggerFilterObjectWithRefs& filterproduct) const {
+  bool accept = false;
+  int ijet = 0;
+  edm::Handle<reco::CaloJetCollection> jets;
+  iEvent.getByToken(jetInputToken, jets);
+  edm::Handle<edm::ValueMap<float>> jetTimes;
+  iEvent.getByToken(jetTimesInputToken, jetTimes);
+  edm::Handle<edm::ValueMap<unsigned int>> jetCellsForTiming;
+  iEvent.getByToken(jetCellsForTimingInputToken, jetCellsForTiming);
+  edm::Handle<edm::ValueMap<float>> jetEcalEtForTiming;
+  iEvent.getByToken(jetEcalEtForTimingInputToken, jetEcalEtForTiming);
+  unsigned int njets = 0;
+  for (auto const& c : *jets) {
+    reco::CaloJetRef calojetref(jets, ijet);
+    if ((*jetTimes)[calojetref] > jetTimeThresh_ && (*jetEcalEtForTiming)[calojetref] > jetEcalEtForTimingThresh_ &&
+        (*jetCellsForTiming)[calojetref] > jetCellsForTimingThresh_ && c.pt() > minPt_)
+      njets++;
+    ijet++;
+  }
+  accept = njets >= minJets_;
+  return accept;
+}
+
+// Fill descriptions
+void HLTCaloJetTimingFilter::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.add("caloJetTimingFilter", desc);
+}
+
+// declare this class as a framework plugin
+DEFINE_FWK_MODULE(HLTCaloJetTimingFilter);

HLTrigger/JetMET/plugins/HLTCaloJetTimingProducer.cc

@@ -0,0 +1,172 @@
+/** \class HLTCaloJetTimingProducer
+ *
+ *  \brief  This produces timing and associated ecal cell information for calo jets 
+ *  \author Matthew Citron
+ *
+ *
+ */
+
+// 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 "FWCore/Utilities/interface/StreamID.h"
+#include "DataFormats/Common/interface/ValueMap.h"
+
+#include "DataFormats/JetReco/interface/CaloJetCollection.h"
+
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
+#include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
+#include "TLorentzVector.h"
+#include "DataFormats/Math/interface/deltaR.h"
+
+//
+// class declaration
+//
+class HLTCaloJetTimingProducer : public edm::stream::EDProducer<> {
+public:
+  explicit HLTCaloJetTimingProducer(const edm::ParameterSet&);
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+  // ----------member data ---------------------------
+  // Input collections
+  edm::InputTag jetLabel_;
+  edm::InputTag ecalEBLabel_;
+  edm::InputTag ecalEELabel_;
+  // Include endcap jets or only barrel
+  bool barrelOnly_;
+
+  edm::EDGetTokenT<reco::CaloJetCollection> jetInputToken;
+  edm::EDGetTokenT<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEBToken;
+  edm::EDGetTokenT<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEEToken;
+};
+
+//Constructor
+HLTCaloJetTimingProducer::HLTCaloJetTimingProducer(const edm::ParameterSet& iConfig) {
+  produces<edm::ValueMap<float>>("");
+  produces<edm::ValueMap<unsigned int>>("jetCellsForTiming");
+  produces<edm::ValueMap<float>>("jetEcalEtForTiming");
+  jetLabel_ = iConfig.getParameter<edm::InputTag>("jets");
+  ecalEBLabel_ = iConfig.getParameter<edm::InputTag>("ebRecHitsColl");
+  ecalEELabel_ = iConfig.getParameter<edm::InputTag>("eeRecHitsColl");
+  barrelOnly_ = iConfig.getParameter<bool>("barrelOnly");
+  jetInputToken = consumes<std::vector<reco::CaloJet>>(jetLabel_);
+  ecalRecHitsEBToken = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>>(ecalEBLabel_);
+  ecalRecHitsEEToken = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>>(ecalEELabel_);
+}
+
+//Producer
+void HLTCaloJetTimingProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  using namespace edm;
+  std::vector<float> jetTimings;
+  std::vector<unsigned int> jetCellsForTiming;
+  std::vector<float> jetEcalEtForTiming;
+  Handle<reco::CaloJetCollection> jets;
+  iEvent.getByToken(jetInputToken, jets);
+  Handle<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEB;
+  iEvent.getByToken(ecalRecHitsEBToken, ecalRecHitsEB);
+  Handle<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEE;
+  iEvent.getByToken(ecalRecHitsEEToken, ecalRecHitsEE);
+  edm::ESHandle<CaloGeometry> pG;
+  iSetup.get<CaloGeometryRecord>().get(pG);
+  for (auto const& c : *jets) {
+    int iCell = -1;
+    float weightedTimeCell = 0;
+    float totalEmEnergyCell = 0;
+    unsigned int nCells = 0;
+    for (EcalRecHitCollection::const_iterator i = ecalRecHitsEB->begin(); i != ecalRecHitsEB->end(); i++) {
+      iCell++;
+      if ((i->checkFlag(EcalRecHit::kSaturated) || i->checkFlag(EcalRecHit::kLeadingEdgeRecovered) ||
+           i->checkFlag(EcalRecHit::kPoorReco) || i->checkFlag(EcalRecHit::kWeird) ||
+           i->checkFlag(EcalRecHit::kDiWeird)))
+        continue;
+      if (i->energy() < 0.5)
+        continue;
+      if (i->timeError() <= 0. || i->timeError() > 100)
+        continue;
+      if (i->time() < -12.5 || i->time() > 12.5)
+        continue;
+      GlobalPoint p = pG->getPosition(i->detid());
+      if (reco::deltaR(c, p) > 0.4)
+        continue;
+      weightedTimeCell += i->time() * i->energy() * sin(p.theta());
+      totalEmEnergyCell += i->energy() * sin(p.theta());
+      nCells++;
+    }
+    iCell = -1;
+    if (!barrelOnly_) {
+      for (EcalRecHitCollection::const_iterator i = ecalRecHitsEE->begin(); i != ecalRecHitsEE->end(); i++) {
+        iCell++;
+        if ((i->checkFlag(EcalRecHit::kSaturated) || i->checkFlag(EcalRecHit::kLeadingEdgeRecovered) ||
+             i->checkFlag(EcalRecHit::kPoorReco) || i->checkFlag(EcalRecHit::kWeird) ||
+             i->checkFlag(EcalRecHit::kDiWeird)))
+          continue;
+        if (i->energy() < 0.5)
+          continue;
+        if (i->timeError() <= 0. || i->timeError() > 100)
+          continue;
+        if (i->time() < -12.5 || i->time() > 12.5)
+          continue;
+        GlobalPoint p = pG->getPosition(i->detid());
+        if (reco::deltaR(c, p) > 0.4)
+          continue;
+        weightedTimeCell += i->time() * i->energy() * sin(p.theta());
+        totalEmEnergyCell += i->energy() * sin(p.theta());
+        nCells++;
+      }
+    }
+    //If there is at least one ecal cell passing selection
+    // calculate timing
+    if (totalEmEnergyCell > 0) {
+      jetTimings.push_back(weightedTimeCell / totalEmEnergyCell);
+      jetEcalEtForTiming.push_back(totalEmEnergyCell);
+      jetCellsForTiming.push_back(nCells);
+    } else {
+      jetTimings.push_back(-50);
+      jetEcalEtForTiming.push_back(totalEmEnergyCell);
+      jetCellsForTiming.push_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
+void HLTCaloJetTimingProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("jets", edm::InputTag(""));
+  desc.add<bool>("barrelOnly", false);
+  desc.add<edm::InputTag>("ebRecHitsColl", edm::InputTag("hltEcalRecHit", "EcalRecHitsEB"));
+  desc.add<edm::InputTag>("eeRecHitsColl", edm::InputTag("hltEcalRecHit", "EcalRecHitsEE"));
+  descriptions.add("caloJetTimingProducer", desc);
+}
+
+// declare this class as a framework plugin
+DEFINE_FWK_MODULE(HLTCaloJetTimingProducer);

HLTrigger/JetMET/test/hltCaloJetTimingFilter_cfg.py

@@ -0,0 +1,50 @@
+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("FWCore.MessageLogger.MessageLogger_cfi")
+process.load("FWCore.MessageLogger.MessageLogger_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 = '120X_mcRun3_2021_realistic_v2'
+
+process.hltTimingProducer = cms.EDProducer('HLTCaloJetTimingProducer',
+        jets = cms.InputTag( "ak4CaloJets" ),
+        ebRecHitsColl = cms.InputTag( 'ecalRecHit','EcalRecHitsEB' ),
+        eeRecHitsColl = cms.InputTag( 'ecalRecHit','EcalRecHitsEE' ),
+        barrelOnly = cms.bool(True),
+)
+
+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)