JetProducers/plugins/FixedGridRhoProducerFastjetFromRecHit.cc

0001 /*
0002
0003 Author: Swagata Mukherjee
0004
0005 Date: November 2021
0006
0007 This producer computes rho from ECAL and HCAL recHits.
0008 The current plan is to use it in egamma and muon HLT, who currently use
0009 the other producer called FixedGridRhoProducerFastjet.
0010 At HLT, FixedGridRhoProducerFastjet takes calotowers as input and computes rho.
0011 But calotowers are expected to be phased out sometime in Run3.
0012 So this recHit-based rho producer, FixedGridRhoProducerFastjetFromRecHit, can be used as an alternative.
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   //Muon HLT currently use ECAL-only rho for ECAL isolation,
0055   //and HCAL-only rho for HCAL isolation. So, this skipping option is needed.
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   // following are needed to grab HCal thresholds from GT
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   //We use raw recHits and not the PF recHits, because in Phase1 muon and egamma HLT,
0093   //PF recHit collection is regional, while the full raw recHit collection is available.
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   //eThresHB/HE are from RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py
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 //HCAL noise cleaning cuts.
0169 bool FixedGridRhoProducerFastjetFromRecHit::passedHcalNoiseCut(const HBHERecHit &hit,
0170                                                                const HcalPFCuts *hcalcuts) const {
0171   if (hcalcuts != nullptr) {  // using hcal cuts from DB
0172     const HcalPFCut *item = hcalcuts->getValues(hit.id().rawId());
0173     return (hit.energy() > item->noiseThreshold());
0174   } else {  // using hcal cuts from config file
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 //ECAL noise cleaning cuts using per-crystal PF-recHit thresholds.
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);