EgammaHLTProducers/plugins/HLTRechitInRegionsProducer.cc

0001 // C/C++ headers
0002 #include <vector>
0003 #include <memory>
0004
0005 // Framework
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0011 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0012
0013 // Reconstruction Classes
0014 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0015 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0016 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0018 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0019 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0020 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0021 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0022
0023 // Geometry
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0026 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0027 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0028 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0029 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0030 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0031
0032 // Level 1 Trigger
0033 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
0034 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
0035 #include "CondFormats/L1TObjects/interface/L1CaloGeometry.h"
0036 #include "CondFormats/DataRecord/interface/L1CaloGeometryRecord.h"
0037
0038 #include "FWCore/Framework/interface/Frameworkfwd.h"
0039 #include "FWCore/Framework/interface/stream/EDProducer.h"
0040 #include "FWCore/Framework/interface/Event.h"
0041 #include "FWCore/Framework/interface/EventSetup.h"
0042 #include "FWCore/Framework/interface/MakerMacros.h"
0043
0044 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0045 #include "FWCore/Utilities/interface/EDGetToken.h"
0046
0047 // Reco candidates
0048 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0049 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0050
0051 #include "DataFormats/Math/interface/RectangularEtaPhiRegion.h"
0052 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0053
0054 // Geometry and topology
0055 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0056 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0057 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0058 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0059 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0060 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0061 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0062
0063 // Level 1 Trigger
0064 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
0065 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
0066 #include "CondFormats/L1TObjects/interface/L1CaloGeometry.h"
0067 #include "CondFormats/DataRecord/interface/L1CaloGeometryRecord.h"
0068
0069 template <typename T1>
0070 class HLTRechitInRegionsProducer : public edm::stream::EDProducer<> {
0071   typedef std::vector<T1> T1Collection;
0072   typedef typename T1::const_iterator T1iterator;
0073
0074 public:
0075   HLTRechitInRegionsProducer(const edm::ParameterSet& ps);
0076   ~HLTRechitInRegionsProducer() override = default;
0077
0078   void produce(edm::Event&, edm::EventSetup const&) override;
0079   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0080
0081 private:
0082   void eventSetupConsumes();
0083
0084   void getEtaPhiRegions(std::vector<RectangularEtaPhiRegion>*, T1Collection const&, edm::EventSetup const&, bool const);
0085
0086   const bool useUncalib_;
0087
0088   const bool doIsolated_;
0089
0090   const edm::EDGetTokenT<T1Collection> l1TokenIsolated_;
0091   const edm::EDGetTokenT<T1Collection> l1TokenNonIsolated_;
0092   const double l1LowerThr_;
0093   const double l1UpperThr_;
0094   const double l1LowerThrIgnoreIsolation_;
0095
0096   const double regionEtaMargin_;
0097   const double regionPhiMargin_;
0098
0099   const std::vector<edm::InputTag> hitLabels;
0100   const std::vector<std::string> productLabels;
0101
0102   std::vector<edm::EDGetTokenT<EcalRecHitCollection>> hitTokens;
0103   std::vector<edm::EDGetTokenT<EcalUncalibratedRecHitCollection>> uncalibHitTokens;
0104
0105   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0106   edm::ESGetToken<L1CaloGeometry, L1CaloGeometryRecord> l1CaloGeometryToken_;
0107 };
0108
0109 template <typename T1>
0110 HLTRechitInRegionsProducer<T1>::HLTRechitInRegionsProducer(const edm::ParameterSet& ps)
0111     : useUncalib_(ps.getParameter<bool>("useUncalib")),
0112       doIsolated_(ps.getParameter<bool>("doIsolated")),
0113       l1TokenIsolated_(doIsolated_ ? consumes<T1Collection>(ps.getParameter<edm::InputTag>("l1TagIsolated"))
0114                                    : edm::EDGetTokenT<T1Collection>()),
0115       l1TokenNonIsolated_(consumes<T1Collection>(ps.getParameter<edm::InputTag>("l1TagNonIsolated"))),
0116       l1LowerThr_(ps.getParameter<double>("l1LowerThr")),
0117       l1UpperThr_(ps.getParameter<double>("l1UpperThr")),
0118       l1LowerThrIgnoreIsolation_(ps.getParameter<double>("l1LowerThrIgnoreIsolation")),
0119       regionEtaMargin_(ps.getParameter<double>("regionEtaMargin")),
0120       regionPhiMargin_(ps.getParameter<double>("regionPhiMargin")),
0121       hitLabels(ps.getParameter<std::vector<edm::InputTag>>("ecalhitLabels")),
0122       productLabels(ps.getParameter<std::vector<std::string>>("productLabels")) {
0123   eventSetupConsumes();
0124
0125   if (useUncalib_) {
0126     for (unsigned int i = 0; i < hitLabels.size(); i++) {
0127       uncalibHitTokens.push_back(consumes<EcalUncalibratedRecHitCollection>(hitLabels[i]));
0128       produces<EcalUncalibratedRecHitCollection>(productLabels[i]);
0129     }
0130   } else {
0131     for (unsigned int i = 0; i < hitLabels.size(); i++) {
0132       hitTokens.push_back(consumes<EcalRecHitCollection>(hitLabels[i]));
0133       produces<EcalRecHitCollection>(productLabels[i]);
0134     }
0135   }
0136 }
0137
0138 template <>
0139 void HLTRechitInRegionsProducer<l1extra::L1EmParticle>::eventSetupConsumes() {
0140   caloGeometryToken_ = esConsumes();
0141   l1CaloGeometryToken_ = esConsumes();
0142 }
0143
0144 template <typename T1>
0145 void HLTRechitInRegionsProducer<T1>::eventSetupConsumes() {
0146   caloGeometryToken_ = esConsumes();
0147 }
0148
0149 template <typename T1>
0150 void HLTRechitInRegionsProducer<T1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0151   edm::ParameterSetDescription desc;
0152   std::vector<std::string> productTags;
0153   productTags.push_back("EcalRegionalRecHitsEB");
0154   productTags.push_back("EcalRegionalRecHitsEE");
0155   desc.add<std::vector<std::string>>("productLabels", productTags);
0156   std::vector<edm::InputTag> inputTags;
0157   inputTags.push_back(edm::InputTag("hltEcalRegionalEgammaRecHit:EcalRecHitsEB"));
0158   inputTags.push_back(edm::InputTag("hltEcalRegionalEgammaRecHit:EcalRecHitsEE"));
0159   inputTags.push_back(edm::InputTag("hltESRegionalEgammaRecHit:EcalRecHitsES"));
0160   desc.add<std::vector<edm::InputTag>>("ecalhitLabels", inputTags);
0161   desc.add<edm::InputTag>("l1TagIsolated", edm::InputTag("l1extraParticles", "Isolated"));
0162   desc.add<edm::InputTag>("l1TagNonIsolated", edm::InputTag("l1extraParticles", "NonIsolated"));
0163   desc.add<bool>("useUncalib", true);
0164   desc.add<bool>("doIsolated", true);
0165   desc.add<double>("l1LowerThr", 5.0);
0166   desc.add<double>("l1UpperThr", 999.);
0167   desc.add<double>("l1LowerThrIgnoreIsolation", 0.0);
0168   desc.add<double>("regionEtaMargin", 0.14);
0169   desc.add<double>("regionPhiMargin", 0.4);
0170   descriptions.add(defaultModuleLabel<HLTRechitInRegionsProducer<T1>>(), desc);
0171 }
0172
0173 template <typename T1>
0174 void HLTRechitInRegionsProducer<T1>::produce(edm::Event& evt, edm::EventSetup const& eventSetup) {
0175   // get the collection geometry:
0176   auto const& geometry = eventSetup.getData(caloGeometryToken_);
0177   const CaloSubdetectorGeometry* geometry_p;
0178   std::unique_ptr<const CaloSubdetectorTopology> topology;
0179
0180   //Get the L1 EM Particle Collection
0181   edm::Handle<T1Collection> emIsolColl;
0182   if (doIsolated_) {
0183     evt.getByToken(l1TokenIsolated_, emIsolColl);
0184   }
0185
0186   std::vector<RectangularEtaPhiRegion> regions;
0187   if (doIsolated_)
0188     getEtaPhiRegions(&regions, *emIsolColl, eventSetup, true);
0189
0190   if (!doIsolated_ or (l1LowerThrIgnoreIsolation_ < 64))
0191     getEtaPhiRegions(&regions, evt.get(l1TokenNonIsolated_), eventSetup, false);
0192
0193   if (useUncalib_) {
0194     edm::Handle<EcalUncalibratedRecHitCollection> urhcH[3];
0195     for (unsigned int i = 0; i < hitLabels.size(); i++) {
0196       auto uhits = std::make_unique<EcalUncalibratedRecHitCollection>();
0197
0198       evt.getByToken(uncalibHitTokens[i], urhcH[i]);
0199       if (!(urhcH[i].isValid())) {
0200         edm::LogError("ProductNotFound") << "could not get a handle on the EcalRecHitCollection! ("
0201                                          << hitLabels[i].encode() << ")" << std::endl;
0202         return;
0203       }
0204       const EcalUncalibratedRecHitCollection* uncalibRecHits = urhcH[i].product();
0205
0206       if (!uncalibRecHits->empty()) {
0207         if ((*uncalibRecHits)[0].id().subdetId() == EcalBarrel) {
0208           geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0209           topology = std::make_unique<EcalBarrelTopology>(geometry);
0210         } else if ((*uncalibRecHits)[0].id().subdetId() == EcalEndcap) {
0211           geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0212           topology = std::make_unique<EcalEndcapTopology>(geometry);
0213         } else if ((*uncalibRecHits)[0].id().subdetId() == EcalPreshower) {
0214           geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0215           topology = std::make_unique<EcalPreshowerTopology>();
0216         } else
0217           throw(std::runtime_error("\n\nProducer encountered invalied ecalhitcollection type.\n\n"));
0218
0219         if (!regions.empty()) {
0220           EcalUncalibratedRecHitCollection::const_iterator it;
0221
0222           for (it = uncalibRecHits->begin(); it != uncalibRecHits->end(); it++) {
0223             auto this_cell = geometry_p->getGeometry(it->id());
0224
0225             std::vector<RectangularEtaPhiRegion>::const_iterator region;
0226             for (region = regions.begin(); region != regions.end(); region++) {
0227               if (region->inRegion(this_cell->etaPos(), this_cell->phiPos())) {
0228                 uhits->push_back(*it);
0229                 break;
0230               }
0231             }
0232           }
0233         }
0234       }
0235       evt.put(std::move(uhits), productLabels[i]);
0236     }
0237
0238   } else {
0239     edm::Handle<EcalRecHitCollection> rhcH[3];
0240     for (unsigned int i = 0; i < hitLabels.size(); i++) {
0241       auto hits = std::make_unique<EcalRecHitCollection>();
0242
0243       evt.getByToken(hitTokens[i], rhcH[i]);
0244       if (!(rhcH[i].isValid())) {
0245         edm::LogError("ProductNotFound") << "could not get a handle on the EcalRecHitCollection! ("
0246                                          << hitLabels[i].encode() << ")" << std::endl;
0247         return;
0248       }
0249       const EcalRecHitCollection* recHits = rhcH[i].product();
0250
0251       if (!recHits->empty()) {
0252         if ((*recHits)[0].id().subdetId() == EcalBarrel) {
0253           geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0254           topology = std::make_unique<EcalBarrelTopology>(geometry);
0255         } else if ((*recHits)[0].id().subdetId() == EcalEndcap) {
0256           geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0257           topology = std::make_unique<EcalEndcapTopology>(geometry);
0258         } else if ((*recHits)[0].id().subdetId() == EcalPreshower) {
0259           geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0260           topology = std::make_unique<EcalPreshowerTopology>();
0261         } else
0262           throw(std::runtime_error("\n\nProducer encountered invalied ecalhitcollection type.\n\n"));
0263
0264         if (!regions.empty()) {
0265           EcalRecHitCollection::const_iterator it;
0266           for (it = recHits->begin(); it != recHits->end(); it++) {
0267             auto this_cell = geometry_p->getGeometry(it->id());
0268
0269             std::vector<RectangularEtaPhiRegion>::const_iterator region;
0270             for (region = regions.begin(); region != regions.end(); region++) {
0271               if (region->inRegion(this_cell->etaPos(), this_cell->phiPos())) {
0272                 hits->push_back(*it);
0273                 break;
0274               }
0275             }
0276           }
0277         }
0278       }
0279       evt.put(std::move(hits), productLabels[i]);
0280     }
0281   }
0282 }
0283
0284 template <>
0285 void HLTRechitInRegionsProducer<l1extra::L1EmParticle>::getEtaPhiRegions(
0286     std::vector<RectangularEtaPhiRegion>* theRegions,
0287     T1Collection const& theCandidateCollection,
0288     edm::EventSetup const& eventSetup,
0289     bool const isolatedCase) {
0290   auto const& l1CaloGeom = eventSetup.getData(l1CaloGeometryToken_);
0291
0292   for (unsigned int candItr = 0; candItr < theCandidateCollection.size(); candItr++) {
0293     const l1extra::L1EmParticle& emItr = theCandidateCollection.at(candItr);
0294
0295     if (!isolatedCase) {
0296       if (doIsolated_ and (emItr.et() < l1LowerThrIgnoreIsolation_))
0297         continue;
0298     }
0299
0300     if ((emItr.et() > l1LowerThr_) and (emItr.et() < l1UpperThr_)) {
0301       // Access the GCT hardware object corresponding to the L1Extra EM object.
0302       int etaIndex = emItr.gctEmCand()->etaIndex();
0303       int phiIndex = emItr.gctEmCand()->phiIndex();
0304
0305       // Use the L1CaloGeometry to find the eta, phi bin boundaries.
0306       double etaLow = l1CaloGeom.etaBinLowEdge(etaIndex);
0307       double etaHigh = l1CaloGeom.etaBinHighEdge(etaIndex);
0308       double phiLow = l1CaloGeom.emJetPhiBinLowEdge(phiIndex);
0309       double phiHigh = l1CaloGeom.emJetPhiBinHighEdge(phiIndex);
0310
0311       etaLow -= regionEtaMargin_;
0312       etaHigh += regionEtaMargin_;
0313       phiLow -= regionPhiMargin_;
0314       phiHigh += regionPhiMargin_;
0315
0316       theRegions->push_back(RectangularEtaPhiRegion(etaLow, etaHigh, phiLow, phiHigh));
0317     }
0318   }
0319 }
0320
0321 template <typename T1>
0322 void HLTRechitInRegionsProducer<T1>::getEtaPhiRegions(std::vector<RectangularEtaPhiRegion>* theRegions,
0323                                                       T1Collection const& theCandidateCollection,
0324                                                       edm::EventSetup const&,
0325                                                       bool const) {
0326   for (unsigned int candItr = 0; candItr < theCandidateCollection.size(); candItr++) {
0327     const T1& emItr = theCandidateCollection.at(candItr);
0328     if ((emItr.et() > l1LowerThr_) and (emItr.et() < l1UpperThr_)) {
0329       double etaLow = emItr.eta() - regionEtaMargin_;
0330       double etaHigh = emItr.eta() + regionEtaMargin_;
0331       double phiLow = emItr.phi() - regionPhiMargin_;
0332       double phiHigh = emItr.phi() + regionPhiMargin_;
0333
0334       theRegions->push_back(RectangularEtaPhiRegion(etaLow, etaHigh, phiLow, phiHigh));
0335     }
0336   }
0337 }
0338
0339 typedef HLTRechitInRegionsProducer<l1extra::L1EmParticle> EgammaHLTRechitInRegionsProducer;
0340 DEFINE_FWK_MODULE(EgammaHLTRechitInRegionsProducer);
0341
0342 typedef HLTRechitInRegionsProducer<reco::RecoChargedCandidate> MuonHLTRechitInRegionsProducer;
0343 DEFINE_FWK_MODULE(MuonHLTRechitInRegionsProducer);