JetMET/plugins/HLTJetTimingProducer.h

0001 /** \class HLTJetTimingProducer
0002  *
0003  *  \brief  This produces timing and associated ecal cell information for calo jets
0004  *  \author Matthew Citron
0005  *
0006  */
0007 #ifndef HLTrigger_JetMET_plugins_HLTJetTimingProducer_h
0008 #define HLTrigger_JetMET_plugins_HLTJetTimingProducer_h
0009
0010 #include <memory>
0011 #include <cmath>
0012
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/stream/EDProducer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "DataFormats/Common/interface/ValueMap.h"
0020 #include "DataFormats/Common/interface/SortedCollection.h"
0021 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0022 #include "DataFormats/Math/interface/deltaR.h"
0023 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0026
0027 template <typename T>
0028 class HLTJetTimingProducer : public edm::stream::EDProducer<> {
0029 public:
0030   explicit HLTJetTimingProducer(const edm::ParameterSet&);
0031   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0032
0033 private:
0034   void produce(edm::Event&, const edm::EventSetup&) override;
0035   void jetTimeFromEcalCells(const T&,
0036                             const edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>&,
0037                             const CaloGeometry&,
0038                             float&,
0039                             float&,
0040                             uint&);
0041
0042   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0043   // Input collections
0044   const edm::EDGetTokenT<std::vector<T>> jetInputToken_;
0045   const edm::EDGetTokenT<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEBToken_;
0046   const edm::EDGetTokenT<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>> ecalRecHitsEEToken_;
0047
0048   // Include barrel, endcap jets or both
0049   const bool barrelJets_;
0050   const bool endcapJets_;
0051
0052   // Configurables for timing definition
0053   const double ecalCellEnergyThresh_;
0054   const double ecalCellTimeThresh_;
0055   const double ecalCellTimeErrorThresh_;
0056   const double matchingRadius2_;
0057 };
0058
0059 template <typename T>
0060 HLTJetTimingProducer<T>::HLTJetTimingProducer(const edm::ParameterSet& iConfig)
0061     : caloGeometryToken_(esConsumes()),
0062       jetInputToken_{consumes<std::vector<T>>(iConfig.getParameter<edm::InputTag>("jets"))},
0063       ecalRecHitsEBToken_{consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>>(
0064           iConfig.getParameter<edm::InputTag>("ebRecHitsColl"))},
0065       ecalRecHitsEEToken_{consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>>(
0066           iConfig.getParameter<edm::InputTag>("eeRecHitsColl"))},
0067       barrelJets_{iConfig.getParameter<bool>("barrelJets")},
0068       endcapJets_{iConfig.getParameter<bool>("endcapJets")},
0069       ecalCellEnergyThresh_{iConfig.getParameter<double>("ecalCellEnergyThresh")},
0070       ecalCellTimeThresh_{iConfig.getParameter<double>("ecalCellTimeThresh")},
0071       ecalCellTimeErrorThresh_{iConfig.getParameter<double>("ecalCellTimeErrorThresh")},
0072       matchingRadius2_{std::pow(iConfig.getParameter<double>("matchingRadius"), 2)} {
0073   produces<edm::ValueMap<float>>("");
0074   produces<edm::ValueMap<unsigned int>>("jetCellsForTiming");
0075   produces<edm::ValueMap<float>>("jetEcalEtForTiming");
0076 }
0077
0078 template <typename T>
0079 void HLTJetTimingProducer<T>::jetTimeFromEcalCells(
0080     const T& jet,
0081     const edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit>>& ecalRecHits,
0082     const CaloGeometry& caloGeometry,
0083     float& weightedTimeCell,
0084     float& totalEmEnergyCell,
0085     uint& nCells) {
0086   for (auto const& ecalRH : ecalRecHits) {
0087     if (ecalRH.checkFlag(EcalRecHit::kSaturated) || ecalRH.checkFlag(EcalRecHit::kLeadingEdgeRecovered) ||
0088         ecalRH.checkFlag(EcalRecHit::kPoorReco) || ecalRH.checkFlag(EcalRecHit::kWeird) ||
0089         ecalRH.checkFlag(EcalRecHit::kDiWeird))
0090       continue;
0091     if (ecalRH.energy() < ecalCellEnergyThresh_)
0092       continue;
0093     if (ecalRH.timeError() <= 0. || ecalRH.timeError() > ecalCellTimeErrorThresh_)
0094       continue;
0095     if (std::abs(ecalRH.time()) > ecalCellTimeThresh_)
0096       continue;
0097     auto const pos = caloGeometry.getPosition(ecalRH.detid());
0098     if (reco::deltaR2(jet, pos) > matchingRadius2_)
0099       continue;
0100     weightedTimeCell += ecalRH.time() * ecalRH.energy() * sin(pos.theta());
0101     totalEmEnergyCell += ecalRH.energy() * sin(pos.theta());
0102     nCells++;
0103   }
0104   if (totalEmEnergyCell > 0) {
0105     weightedTimeCell /= totalEmEnergyCell;
0106   }
0107 }
0108
0109 template <typename T>
0110 void HLTJetTimingProducer<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0111   auto const& caloGeometry = iSetup.getData(caloGeometryToken_);
0112   auto const jets = iEvent.getHandle(jetInputToken_);
0113   auto const& ecalRecHitsEB = iEvent.get(ecalRecHitsEBToken_);
0114   auto const& ecalRecHitsEE = iEvent.get(ecalRecHitsEEToken_);
0115
0116   std::vector<float> jetTimings;
0117   std::vector<unsigned int> jetCellsForTiming;
0118   std::vector<float> jetEcalEtForTiming;
0119
0120   jetTimings.reserve(jets->size());
0121   jetEcalEtForTiming.reserve(jets->size());
0122   jetCellsForTiming.reserve(jets->size());
0123
0124   for (auto const& jet : *jets) {
0125     float weightedTimeCell = 0;
0126     float totalEmEnergyCell = 0;
0127     unsigned int nCells = 0;
0128     if (barrelJets_)
0129       jetTimeFromEcalCells(jet, ecalRecHitsEB, caloGeometry, weightedTimeCell, totalEmEnergyCell, nCells);
0130     if (endcapJets_) {
0131       weightedTimeCell *= totalEmEnergyCell;
0132       jetTimeFromEcalCells(jet, ecalRecHitsEE, caloGeometry, weightedTimeCell, totalEmEnergyCell, nCells);
0133     }
0134
0135     // If there is at least one ecal cell passing selection, calculate timing
0136     jetTimings.emplace_back(totalEmEnergyCell > 0 ? weightedTimeCell : -50);
0137     jetEcalEtForTiming.emplace_back(totalEmEnergyCell);
0138     jetCellsForTiming.emplace_back(nCells);
0139   }
0140
0141   std::unique_ptr<edm::ValueMap<float>> jetTimings_out(new edm::ValueMap<float>());
0142   edm::ValueMap<float>::Filler jetTimings_filler(*jetTimings_out);
0143   jetTimings_filler.insert(jets, jetTimings.begin(), jetTimings.end());
0144   jetTimings_filler.fill();
0145   iEvent.put(std::move(jetTimings_out), "");
0146
0147   std::unique_ptr<edm::ValueMap<float>> jetEcalEtForTiming_out(new edm::ValueMap<float>());
0148   edm::ValueMap<float>::Filler jetEcalEtForTiming_filler(*jetEcalEtForTiming_out);
0149   jetEcalEtForTiming_filler.insert(jets, jetEcalEtForTiming.begin(), jetEcalEtForTiming.end());
0150   jetEcalEtForTiming_filler.fill();
0151   iEvent.put(std::move(jetEcalEtForTiming_out), "jetEcalEtForTiming");
0152
0153   std::unique_ptr<edm::ValueMap<unsigned int>> jetCellsForTiming_out(new edm::ValueMap<unsigned int>());
0154   edm::ValueMap<unsigned int>::Filler jetCellsForTiming_filler(*jetCellsForTiming_out);
0155   jetCellsForTiming_filler.insert(jets, jetCellsForTiming.begin(), jetCellsForTiming.end());
0156   jetCellsForTiming_filler.fill();
0157   iEvent.put(std::move(jetCellsForTiming_out), "jetCellsForTiming");
0158 }
0159
0160 template <typename T>
0161 void HLTJetTimingProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0162   edm::ParameterSetDescription desc;
0163   desc.add<edm::InputTag>("jets", edm::InputTag(""));
0164   desc.add<bool>("barrelJets", false);
0165   desc.add<bool>("endcapJets", false);
0166   desc.add<double>("ecalCellEnergyThresh", 0.5);
0167   desc.add<double>("ecalCellTimeThresh", 12.5);
0168   desc.add<double>("ecalCellTimeErrorThresh", 100.);
0169   desc.add<double>("matchingRadius", 0.4);
0170   desc.add<edm::InputTag>("ebRecHitsColl", edm::InputTag("hltEcalRecHit", "EcalRecHitsEB"));
0171   desc.add<edm::InputTag>("eeRecHitsColl", edm::InputTag("hltEcalRecHit", "EcalRecHitsEE"));
0172   descriptions.add(defaultModuleLabel<HLTJetTimingProducer<T>>(), desc);
0173 }
0174
0175 #endif  // HLTrigger_JetMET_plugins_HLTJetTimingProducer_h