File indexing completed on 2023-11-29 01:33:08
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include "FWCore/Framework/interface/stream/EDProducer.h"
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0022 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0023 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0024 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0025 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0026 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0027 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0028 #include "fastjet/tools/GridMedianBackgroundEstimator.hh"
0029 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
0030 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0031 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0032
0033 class FixedGridRhoProducerFastjetFromRecHit : public edm::stream::EDProducer<> {
0034 public:
0035 explicit FixedGridRhoProducerFastjetFromRecHit(const edm::ParameterSet &iConfig);
0036 ~FixedGridRhoProducerFastjetFromRecHit() override;
0037 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0038
0039 private:
0040 void beginRun(edm::Run const &, edm::EventSetup const &) override;
0041 void produce(edm::Event &, const edm::EventSetup &) override;
0042 std::array<double, 4> getHitP4(const DetId &detId, const double hitE, const CaloGeometry &caloGeometry) const;
0043 bool passedHcalNoiseCut(const HBHERecHit &hit, const HcalPFCuts *) const;
0044 bool passedEcalNoiseCut(const EcalRecHit &hit, const EcalPFRecHitThresholds &thresholds) const;
0045
0046 fastjet::GridMedianBackgroundEstimator bge_;
0047 const edm::EDGetTokenT<HBHERecHitCollection> hbheRecHitsTag_;
0048 const edm::EDGetTokenT<EcalRecHitCollection> ebRecHitsTag_;
0049 const edm::EDGetTokenT<EcalRecHitCollection> eeRecHitsTag_;
0050
0051 const EgammaHcalIsolation::arrayHB eThresHB_;
0052 const EgammaHcalIsolation::arrayHE eThresHE_;
0053
0054
0055
0056 const bool skipHCAL_;
0057 const bool skipECAL_;
0058
0059 const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> ecalPFRecHitThresholdsToken_;
0060 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0061
0062
0063 edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0064 const bool cutsFromDB_;
0065 HcalPFCuts const *paramPF_ = nullptr;
0066 };
0067
0068 FixedGridRhoProducerFastjetFromRecHit::FixedGridRhoProducerFastjetFromRecHit(const edm::ParameterSet &iConfig)
0069 : bge_(iConfig.getParameter<double>("maxRapidity"), iConfig.getParameter<double>("gridSpacing")),
0070 hbheRecHitsTag_(consumes(iConfig.getParameter<edm::InputTag>("hbheRecHitsTag"))),
0071 ebRecHitsTag_(consumes(iConfig.getParameter<edm::InputTag>("ebRecHitsTag"))),
0072 eeRecHitsTag_(consumes(iConfig.getParameter<edm::InputTag>("eeRecHitsTag"))),
0073 eThresHB_(iConfig.getParameter<EgammaHcalIsolation::arrayHB>("eThresHB")),
0074 eThresHE_(iConfig.getParameter<EgammaHcalIsolation::arrayHE>("eThresHE")),
0075 skipHCAL_(iConfig.getParameter<bool>("skipHCAL")),
0076 skipECAL_(iConfig.getParameter<bool>("skipECAL")),
0077 ecalPFRecHitThresholdsToken_{esConsumes()},
0078 caloGeometryToken_{esConsumes()},
0079 cutsFromDB_(iConfig.getParameter<bool>("usePFThresholdsFromDB")) {
0080 if (skipHCAL_ && skipECAL_) {
0081 throw cms::Exception("FixedGridRhoProducerFastjetFromRecHit")
0082 << "skipHCAL and skipECAL both can't be True. Please make at least one of them False.";
0083 }
0084 if (cutsFromDB_) {
0085 hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd, edm::Transition::BeginRun>(edm::ESInputTag("", "withTopo"));
0086 }
0087 produces<double>();
0088 }
0089
0090 void FixedGridRhoProducerFastjetFromRecHit::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0091 edm::ParameterSetDescription desc;
0092
0093
0094 desc.add<edm::InputTag>("hbheRecHitsTag", edm::InputTag("hltHbhereco"));
0095 desc.add<edm::InputTag>("ebRecHitsTag", edm::InputTag("hltEcalRecHit", "EcalRecHitsEB"));
0096 desc.add<edm::InputTag>("eeRecHitsTag", edm::InputTag("hltEcalRecHit", "EcalRecHitsEE"));
0097 desc.add<bool>("skipHCAL", false);
0098 desc.add<bool>("skipECAL", false);
0099
0100 desc.add<std::vector<double> >("eThresHB", {0.1, 0.2, 0.3, 0.3});
0101 desc.add<std::vector<double> >("eThresHE", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0102 desc.add<double>("maxRapidity", 2.5);
0103 desc.add<double>("gridSpacing", 0.55);
0104 desc.add<bool>("usePFThresholdsFromDB", true);
0105 descriptions.addWithDefaultLabel(desc);
0106 }
0107
0108 FixedGridRhoProducerFastjetFromRecHit::~FixedGridRhoProducerFastjetFromRecHit() = default;
0109
0110 void FixedGridRhoProducerFastjetFromRecHit::beginRun(edm::Run const &r, edm::EventSetup const &es) {
0111 if (cutsFromDB_) {
0112 paramPF_ = &es.getData(hcalCutsToken_);
0113 }
0114 }
0115
0116 void FixedGridRhoProducerFastjetFromRecHit::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0117 std::vector<fastjet::PseudoJet> inputs;
0118 auto const &thresholds = iSetup.getData(ecalPFRecHitThresholdsToken_);
0119 auto const &caloGeometry = iSetup.getData(caloGeometryToken_);
0120
0121 if (!skipHCAL_) {
0122 auto const &hbheRecHits = iEvent.get(hbheRecHitsTag_);
0123 inputs.reserve(inputs.size() + hbheRecHits.size());
0124 for (const auto &hit : hbheRecHits) {
0125 if (passedHcalNoiseCut(hit, paramPF_)) {
0126 const auto &hitp4 = getHitP4(hit.id(), hit.energy(), caloGeometry);
0127 inputs.emplace_back(fastjet::PseudoJet(hitp4[0], hitp4[1], hitp4[2], hitp4[3]));
0128 }
0129 }
0130 }
0131
0132 if (!skipECAL_) {
0133 auto const &ebRecHits = iEvent.get(ebRecHitsTag_);
0134 inputs.reserve(inputs.size() + ebRecHits.size());
0135 for (const auto &hit : ebRecHits) {
0136 if (passedEcalNoiseCut(hit, thresholds)) {
0137 const auto &hitp4 = getHitP4(hit.id(), hit.energy(), caloGeometry);
0138 inputs.emplace_back(fastjet::PseudoJet(hitp4[0], hitp4[1], hitp4[2], hitp4[3]));
0139 }
0140 }
0141
0142 auto const &eeRecHits = iEvent.get(eeRecHitsTag_);
0143 inputs.reserve(inputs.size() + eeRecHits.size());
0144 for (const auto &hit : eeRecHits) {
0145 if (passedEcalNoiseCut(hit, thresholds)) {
0146 const auto &hitp4 = getHitP4(hit.id(), hit.energy(), caloGeometry);
0147 inputs.emplace_back(fastjet::PseudoJet(hitp4[0], hitp4[1], hitp4[2], hitp4[3]));
0148 }
0149 }
0150 }
0151
0152 bge_.set_particles(inputs);
0153 iEvent.put(std::make_unique<double>(bge_.rho()));
0154 }
0155
0156 std::array<double, 4> FixedGridRhoProducerFastjetFromRecHit::getHitP4(const DetId &detId,
0157 const double hitE,
0158 const CaloGeometry &caloGeometry) const {
0159 const CaloSubdetectorGeometry *subDetGeom = caloGeometry.getSubdetectorGeometry(detId);
0160 const auto &gpPos = subDetGeom->getGeometry(detId)->repPos();
0161 const double thispt = hitE / cosh(gpPos.eta());
0162 const double thispx = thispt * cos(gpPos.phi());
0163 const double thispy = thispt * sin(gpPos.phi());
0164 const double thispz = thispt * sinh(gpPos.eta());
0165 return std::array<double, 4>{{thispx, thispy, thispz, hitE}};
0166 }
0167
0168
0169 bool FixedGridRhoProducerFastjetFromRecHit::passedHcalNoiseCut(const HBHERecHit &hit,
0170 const HcalPFCuts *hcalcuts) const {
0171 if (hcalcuts != nullptr) {
0172 const HcalPFCut *item = hcalcuts->getValues(hit.id().rawId());
0173 return (hit.energy() > item->noiseThreshold());
0174 } else {
0175 const auto thisDetId = hit.id();
0176 const auto thisDepth = thisDetId.depth();
0177 return (thisDetId.subdet() == HcalBarrel && hit.energy() > eThresHB_[thisDepth - 1]) ||
0178 (thisDetId.subdet() == HcalEndcap && hit.energy() > eThresHE_[thisDepth - 1]);
0179 }
0180 }
0181
0182
0183 bool FixedGridRhoProducerFastjetFromRecHit::passedEcalNoiseCut(const EcalRecHit &hit,
0184 const EcalPFRecHitThresholds &thresholds) const {
0185 return (hit.energy() > thresholds[hit.detid()]);
0186 }
0187
0188 DEFINE_FWK_MODULE(FixedGridRhoProducerFastjetFromRecHit);