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 produce(edm::Event &, const edm::EventSetup &) override;
0041   std::array<double, 4> getHitP4(const DetId &detId, const double hitE, const CaloGeometry &caloGeometry) const;
0042   bool passedHcalNoiseCut(const HBHERecHit &hit, const HcalPFCuts *) const;
0043   bool passedEcalNoiseCut(const EcalRecHit &hit, const EcalPFRecHitThresholds &thresholds) const;
0044
0045   fastjet::GridMedianBackgroundEstimator bge_;
0046   const edm::EDGetTokenT<HBHERecHitCollection> hbheRecHitsTag_;
0047   const edm::EDGetTokenT<EcalRecHitCollection> ebRecHitsTag_;
0048   const edm::EDGetTokenT<EcalRecHitCollection> eeRecHitsTag_;
0049
0050   const EgammaHcalIsolation::arrayHB eThresHB_;
0051   const EgammaHcalIsolation::arrayHE eThresHE_;
0052
0053   //Muon HLT currently use ECAL-only rho for ECAL isolation,
0054   //and HCAL-only rho for HCAL isolation. So, this skipping option is needed.
0055   const bool skipHCAL_;
0056   const bool skipECAL_;
0057
0058   const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> ecalPFRecHitThresholdsToken_;
0059   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0060
0061   // following are needed to grab HCal thresholds from GT
0062   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0063   const bool cutsFromDB_;
0064   HcalPFCuts const *paramPF_ = nullptr;
0065 };
0066
0067 FixedGridRhoProducerFastjetFromRecHit::FixedGridRhoProducerFastjetFromRecHit(const edm::ParameterSet &iConfig)
0068     : bge_(iConfig.getParameter<double>("maxRapidity"), iConfig.getParameter<double>("gridSpacing")),
0069       hbheRecHitsTag_(consumes(iConfig.getParameter<edm::InputTag>("hbheRecHitsTag"))),
0070       ebRecHitsTag_(consumes(iConfig.getParameter<edm::InputTag>("ebRecHitsTag"))),
0071       eeRecHitsTag_(consumes(iConfig.getParameter<edm::InputTag>("eeRecHitsTag"))),
0072       eThresHB_(iConfig.getParameter<EgammaHcalIsolation::arrayHB>("eThresHB")),
0073       eThresHE_(iConfig.getParameter<EgammaHcalIsolation::arrayHE>("eThresHE")),
0074       skipHCAL_(iConfig.getParameter<bool>("skipHCAL")),
0075       skipECAL_(iConfig.getParameter<bool>("skipECAL")),
0076       ecalPFRecHitThresholdsToken_{esConsumes()},
0077       caloGeometryToken_{esConsumes()},
0078       cutsFromDB_(iConfig.getParameter<bool>("usePFThresholdsFromDB")) {
0079   if (skipHCAL_ && skipECAL_) {
0080     throw cms::Exception("FixedGridRhoProducerFastjetFromRecHit")
0081         << "skipHCAL and skipECAL both can't be True. Please make at least one of them False.";
0082   }
0083   if (cutsFromDB_) {
0084     hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"));
0085   }
0086   produces<double>();
0087 }
0088
0089 void FixedGridRhoProducerFastjetFromRecHit::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0090   edm::ParameterSetDescription desc;
0091   //We use raw recHits and not the PF recHits, because in Phase1 muon and egamma HLT,
0092   //PF recHit collection is regional, while the full raw recHit collection is available.
0093   desc.add<edm::InputTag>("hbheRecHitsTag", edm::InputTag("hltHbhereco"));
0094   desc.add<edm::InputTag>("ebRecHitsTag", edm::InputTag("hltEcalRecHit", "EcalRecHitsEB"));
0095   desc.add<edm::InputTag>("eeRecHitsTag", edm::InputTag("hltEcalRecHit", "EcalRecHitsEE"));
0096   desc.add<bool>("skipHCAL", false);
0097   desc.add<bool>("skipECAL", false);
0098   //eThresHB/HE are from RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py
0099   desc.add<std::vector<double> >("eThresHB", {0.1, 0.2, 0.3, 0.3});
0100   desc.add<std::vector<double> >("eThresHE", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0101   desc.add<double>("maxRapidity", 2.5);
0102   desc.add<double>("gridSpacing", 0.55);
0103   desc.add<bool>("usePFThresholdsFromDB", true);
0104   descriptions.addWithDefaultLabel(desc);
0105 }
0106
0107 FixedGridRhoProducerFastjetFromRecHit::~FixedGridRhoProducerFastjetFromRecHit() = default;
0108
0109 void FixedGridRhoProducerFastjetFromRecHit::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0110   if (cutsFromDB_) {
0111     paramPF_ = &iSetup.getData(hcalCutsToken_);
0112   }
0113
0114   std::vector<fastjet::PseudoJet> inputs;
0115   auto const &thresholds = iSetup.getData(ecalPFRecHitThresholdsToken_);
0116   auto const &caloGeometry = iSetup.getData(caloGeometryToken_);
0117
0118   if (!skipHCAL_) {
0119     auto const &hbheRecHits = iEvent.get(hbheRecHitsTag_);
0120     inputs.reserve(inputs.size() + hbheRecHits.size());
0121     for (const auto &hit : hbheRecHits) {
0122       if (passedHcalNoiseCut(hit, paramPF_)) {
0123         const auto &hitp4 = getHitP4(hit.id(), hit.energy(), caloGeometry);
0124         inputs.emplace_back(fastjet::PseudoJet(hitp4[0], hitp4[1], hitp4[2], hitp4[3]));
0125       }
0126     }
0127   }
0128
0129   if (!skipECAL_) {
0130     auto const &ebRecHits = iEvent.get(ebRecHitsTag_);
0131     inputs.reserve(inputs.size() + ebRecHits.size());
0132     for (const auto &hit : ebRecHits) {
0133       if (passedEcalNoiseCut(hit, thresholds)) {
0134         const auto &hitp4 = getHitP4(hit.id(), hit.energy(), caloGeometry);
0135         inputs.emplace_back(fastjet::PseudoJet(hitp4[0], hitp4[1], hitp4[2], hitp4[3]));
0136       }
0137     }
0138
0139     auto const &eeRecHits = iEvent.get(eeRecHitsTag_);
0140     inputs.reserve(inputs.size() + eeRecHits.size());
0141     for (const auto &hit : eeRecHits) {
0142       if (passedEcalNoiseCut(hit, thresholds)) {
0143         const auto &hitp4 = getHitP4(hit.id(), hit.energy(), caloGeometry);
0144         inputs.emplace_back(fastjet::PseudoJet(hitp4[0], hitp4[1], hitp4[2], hitp4[3]));
0145       }
0146     }
0147   }
0148
0149   bge_.set_particles(inputs);
0150   iEvent.put(std::make_unique<double>(bge_.rho()));
0151 }
0152
0153 std::array<double, 4> FixedGridRhoProducerFastjetFromRecHit::getHitP4(const DetId &detId,
0154                                                                       const double hitE,
0155                                                                       const CaloGeometry &caloGeometry) const {
0156   const CaloSubdetectorGeometry *subDetGeom = caloGeometry.getSubdetectorGeometry(detId);
0157   const auto &gpPos = subDetGeom->getGeometry(detId)->repPos();
0158   const double thispt = hitE / cosh(gpPos.eta());
0159   const double thispx = thispt * cos(gpPos.phi());
0160   const double thispy = thispt * sin(gpPos.phi());
0161   const double thispz = thispt * sinh(gpPos.eta());
0162   return std::array<double, 4>{{thispx, thispy, thispz, hitE}};
0163 }
0164
0165 //HCAL noise cleaning cuts.
0166 bool FixedGridRhoProducerFastjetFromRecHit::passedHcalNoiseCut(const HBHERecHit &hit,
0167                                                                const HcalPFCuts *hcalcuts) const {
0168   if (hcalcuts != nullptr) {  // using hcal cuts from DB
0169     const HcalPFCut *item = hcalcuts->getValues(hit.id().rawId());
0170     return (hit.energy() > item->noiseThreshold());
0171   } else {  // using hcal cuts from config file
0172     const auto thisDetId = hit.id();
0173     const auto thisDepth = thisDetId.depth();
0174     return (thisDetId.subdet() == HcalBarrel && hit.energy() > eThresHB_[thisDepth - 1]) ||
0175            (thisDetId.subdet() == HcalEndcap && hit.energy() > eThresHE_[thisDepth - 1]);
0176   }
0177 }
0178
0179 //ECAL noise cleaning cuts using per-crystal PF-recHit thresholds.
0180 bool FixedGridRhoProducerFastjetFromRecHit::passedEcalNoiseCut(const EcalRecHit &hit,
0181                                                                const EcalPFRecHitThresholds &thresholds) const {
0182   return (hit.energy() > thresholds[hit.detid()]);
0183 }
0184
0185 DEFINE_FWK_MODULE(FixedGridRhoProducerFastjetFromRecHit);