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