File indexing completed on 2023-03-17 11:09:32
0001
0002
0003
0004
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
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
0049 const bool barrelJets_;
0050 const bool endcapJets_;
0051
0052
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
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