EgammaHLTProducers/plugins/HLTRecHitInAllL1RegionsProducer.cc

0001 #include <memory>
0002
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/TypeDemangler.h"
0013
0014 // Reco candidates (might not need)
0015 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0016 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0017
0018 // Geometry and topology
0019 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0021 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0022 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0023 #include "DataFormats/Math/interface/RectangularEtaPhiRegion.h"
0024
0025 // Level 1 Trigger
0026 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
0027 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0028 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
0029 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
0030 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0031 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
0032
0033 #include "DataFormats/L1Trigger/interface/EGamma.h"
0034 #include "DataFormats/L1Trigger/interface/Jet.h"
0035 #include "DataFormats/L1Trigger/interface/Muon.h"
0036
0037 #include "CondFormats/L1TObjects/interface/L1CaloGeometry.h"
0038 #include "CondFormats/DataRecord/interface/L1CaloGeometryRecord.h"
0039
0040 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0041 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0042
0043 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0044
0045 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0046
0047 //this is a re-write of HLTRechitInRegionsProducer to be able to handle arbitary L1 collections as inputs
0048 //in the process, some of the cruft was cleaned up but it mantains almost all the old behaviour
0049 //think the only difference now is that it wont throw if its not ECALBarrel, ECALEndcap or ECAL PS rec-hit type
0050 class L1RegionDataBase {
0051 public:
0052   virtual ~L1RegionDataBase() {}
0053   virtual void getEtaPhiRegions(const edm::Event&,
0054                                 const edm::EventSetup&,
0055                                 std::vector<RectangularEtaPhiRegion>&) const = 0;
0056 };
0057
0058 template <typename T1>
0059 class L1RegionData : public L1RegionDataBase {
0060 private:
0061   double const minEt_;
0062   double const maxEt_;
0063   double const regionEtaMargin_;
0064   double const regionPhiMargin_;
0065   edm::EDGetTokenT<T1> const token_;
0066   edm::ESGetToken<L1CaloGeometry, L1CaloGeometryRecord> l1CaloGeometryToken_;
0067
0068   void eventSetupConsumes(edm::ConsumesCollector& consumesColl);
0069
0070 public:
0071   L1RegionData(const edm::ParameterSet& para, edm::ConsumesCollector& consumesColl)
0072       : minEt_(para.getParameter<double>("minEt")),
0073         maxEt_(para.getParameter<double>("maxEt")),
0074         regionEtaMargin_(para.getParameter<double>("regionEtaMargin")),
0075         regionPhiMargin_(para.getParameter<double>("regionPhiMargin")),
0076         token_(consumesColl.consumes<T1>(para.getParameter<edm::InputTag>("inputColl"))) {
0077     eventSetupConsumes(consumesColl);
0078   }
0079
0080   void getEtaPhiRegions(const edm::Event&,
0081                         const edm::EventSetup&,
0082                         std::vector<RectangularEtaPhiRegion>&) const override;
0083   template <typename T2>
0084   bool isEmpty(const T2& coll) const {
0085     return coll.empty();
0086   }
0087   template <typename T2>
0088   static typename T2::const_iterator beginIt(const T2& coll) {
0089     return coll.begin();
0090   }
0091   template <typename T2>
0092   static typename T2::const_iterator endIt(const T2& coll) {
0093     return coll.end();
0094   }
0095   template <typename T2>
0096   bool isEmpty(const BXVector<T2>& coll) const {
0097     return (coll.size() == 0 or coll.isEmpty(0));
0098   }
0099   template <typename T2>
0100   static typename BXVector<T2>::const_iterator beginIt(const BXVector<T2>& coll) {
0101     return coll.begin(0);
0102   }
0103   template <typename T2>
0104   static typename BXVector<T2>::const_iterator endIt(const BXVector<T2>& coll) {
0105     return coll.end(0);
0106   }
0107 };
0108
0109 template <typename RecHitType>
0110 class HLTRecHitInAllL1RegionsProducer : public edm::stream::EDProducer<> {
0111   using RecHitCollectionType = edm::SortedCollection<RecHitType>;
0112
0113 public:
0114   HLTRecHitInAllL1RegionsProducer(const edm::ParameterSet& ps);
0115   ~HLTRecHitInAllL1RegionsProducer() override {}
0116
0117   void produce(edm::Event&, const edm::EventSetup&) override;
0118   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0119
0120 private:
0121   L1RegionDataBase* createL1RegionData(const std::string&,
0122                                        const edm::ParameterSet&,
0123                                        edm::ConsumesCollector&&);  //calling function owns this
0124
0125   std::vector<std::unique_ptr<L1RegionDataBase>> l1RegionData_;
0126
0127   std::vector<edm::InputTag> recHitLabels_;
0128   std::vector<std::string> productLabels_;
0129
0130   std::vector<edm::EDGetTokenT<RecHitCollectionType>> recHitTokens_;
0131
0132   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0133 };
0134
0135 template <typename RecHitType>
0136 HLTRecHitInAllL1RegionsProducer<RecHitType>::HLTRecHitInAllL1RegionsProducer(const edm::ParameterSet& para)
0137     : caloGeometryToken_{esConsumes()} {
0138   const std::vector<edm::ParameterSet> l1InputRegions =
0139       para.getParameter<std::vector<edm::ParameterSet>>("l1InputRegions");
0140   for (auto& pset : l1InputRegions) {
0141     const std::string type = pset.getParameter<std::string>("type");
0142     // meh I was going to use a factory but it was going to be overly complex for my needs
0143     l1RegionData_.emplace_back(createL1RegionData(type, pset, consumesCollector()));
0144   }
0145   recHitLabels_ = para.getParameter<std::vector<edm::InputTag>>("recHitLabels");
0146   productLabels_ = para.getParameter<std::vector<std::string>>("productLabels");
0147
0148   for (unsigned int collNr = 0; collNr < recHitLabels_.size(); collNr++) {
0149     recHitTokens_.push_back(consumes<RecHitCollectionType>(recHitLabels_[collNr]));
0150     produces<RecHitCollectionType>(productLabels_[collNr]);
0151   }
0152 }
0153
0154 template <typename RecHitType>
0155 void HLTRecHitInAllL1RegionsProducer<RecHitType>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0156   edm::ParameterSetDescription desc;
0157   std::vector<std::string> productTags;
0158   productTags.push_back("EcalRegionalRecHitsEB");
0159   productTags.push_back("EcalRegionalRecHitsEE");
0160   desc.add<std::vector<std::string>>("productLabels", productTags);
0161   std::vector<edm::InputTag> recHitLabels;
0162   recHitLabels.push_back(edm::InputTag("hltEcalRegionalEgammaRecHit:EcalRecHitsEB"));
0163   recHitLabels.push_back(edm::InputTag("hltEcalRegionalEgammaRecHit:EcalRecHitsEE"));
0164   recHitLabels.push_back(edm::InputTag("hltESRegionalEgammaRecHit:EcalRecHitsES"));
0165   desc.add<std::vector<edm::InputTag>>("recHitLabels", recHitLabels);
0166   std::vector<edm::ParameterSet> l1InputRegions;
0167
0168   edm::ParameterSet emIsoPSet;
0169   emIsoPSet.addParameter<std::string>("type", "L1EmParticle");
0170   emIsoPSet.addParameter<double>("minEt", 5);
0171   emIsoPSet.addParameter<double>("maxEt", 999);
0172   emIsoPSet.addParameter<double>("regionEtaMargin", 0.14);
0173   emIsoPSet.addParameter<double>("regionPhiMargin", 0.4);
0174   emIsoPSet.addParameter<edm::InputTag>("inputColl", edm::InputTag("hltL1extraParticles:NonIsolated"));
0175   l1InputRegions.push_back(emIsoPSet);
0176   edm::ParameterSet emNonIsoPSet;
0177   emNonIsoPSet.addParameter<std::string>("type", "L1EmParticle");
0178   emNonIsoPSet.addParameter<double>("minEt", 5);
0179   emNonIsoPSet.addParameter<double>("maxEt", 999);
0180   emNonIsoPSet.addParameter<double>("regionEtaMargin", 0.14);
0181   emNonIsoPSet.addParameter<double>("regionPhiMargin", 0.4);
0182   emNonIsoPSet.addParameter<edm::InputTag>("inputColl", edm::InputTag("hltL1extraParticles:Isolated"));
0183   l1InputRegions.push_back(emNonIsoPSet);
0184
0185   // Why no Central Jets here? They are present in the python config, e.g. OnLine_HLT_GRun.py
0186   // SHarper: because these are the default parameters designed to reproduce the original (no jets) behaviour
0187   //
0188   edm::ParameterSet egPSet;
0189   egPSet.addParameter<std::string>("type", "EGamma");
0190   egPSet.addParameter<double>("minEt", 5);
0191   egPSet.addParameter<double>("maxEt", 999);
0192   egPSet.addParameter<double>("regionEtaMargin", 0.4);
0193   egPSet.addParameter<double>("regionPhiMargin", 0.5);
0194   egPSet.addParameter<edm::InputTag>("inputColl", edm::InputTag("hltCaloStage2Digis"));
0195   l1InputRegions.push_back(egPSet);
0196
0197   edm::ParameterSet jetPSet;
0198   jetPSet.addParameter<std::string>("type", "EGamma");
0199   jetPSet.addParameter<double>("minEt", 200);
0200   jetPSet.addParameter<double>("maxEt", 999);
0201   jetPSet.addParameter<double>("regionEtaMargin", 0.4);
0202   jetPSet.addParameter<double>("regionPhiMargin", 0.5);
0203   jetPSet.addParameter<edm::InputTag>("inputColl", edm::InputTag("hltCaloStage2Digis"));
0204   l1InputRegions.push_back(jetPSet);
0205
0206   edm::ParameterSetDescription l1InputRegionDesc;
0207   l1InputRegionDesc.add<std::string>("type");
0208   l1InputRegionDesc.add<double>("minEt");
0209   l1InputRegionDesc.add<double>("maxEt");
0210   l1InputRegionDesc.add<double>("regionEtaMargin");
0211   l1InputRegionDesc.add<double>("regionPhiMargin");
0212   l1InputRegionDesc.add<edm::InputTag>("inputColl");
0213   desc.addVPSet("l1InputRegions", l1InputRegionDesc, l1InputRegions);
0214
0215   descriptions.add(defaultModuleLabel<HLTRecHitInAllL1RegionsProducer<RecHitType>>(), desc);
0216 }
0217
0218 template <typename RecHitType>
0219 void HLTRecHitInAllL1RegionsProducer<RecHitType>::produce(edm::Event& event, const edm::EventSetup& setup) {
0220   // get the collection geometry:
0221   auto const& caloGeom = setup.getData(caloGeometryToken_);
0222
0223   std::vector<RectangularEtaPhiRegion> regions;
0224   std::for_each(l1RegionData_.begin(),
0225                 l1RegionData_.end(),
0226                 [&event, &setup, &regions](const std::unique_ptr<L1RegionDataBase>& input) {
0227                   input->getEtaPhiRegions(event, setup, regions);
0228                 });
0229
0230   for (size_t recHitCollNr = 0; recHitCollNr < recHitTokens_.size(); recHitCollNr++) {
0231     edm::Handle<RecHitCollectionType> recHits;
0232     event.getByToken(recHitTokens_[recHitCollNr], recHits);
0233
0234     if (!(recHits.isValid())) {
0235       edm::LogError("ProductNotFound") << "could not get a handle on the " << typeid(RecHitCollectionType).name()
0236                                        << " named " << recHitLabels_[recHitCollNr].encode() << std::endl;
0237       continue;
0238     }
0239
0240     auto filteredRecHits = std::make_unique<RecHitCollectionType>();
0241
0242     if (!recHits->empty()) {
0243       const CaloSubdetectorGeometry* subDetGeom = caloGeom.getSubdetectorGeometry(recHits->front().id());
0244       if (!regions.empty()) {
0245         for (const RecHitType& recHit : *recHits) {
0246           auto this_cell = subDetGeom->getGeometry(recHit.id());
0247           for (const auto& region : regions) {
0248             if (region.inRegion(this_cell->etaPos(), this_cell->phiPos())) {
0249               filteredRecHits->push_back(recHit);
0250               break;
0251             }
0252           }
0253         }
0254       }  //end check of empty regions
0255     }    //end check of empty rec-hits
0256     //   std::cout <<"putting fileter coll in "<<filteredRecHits->size()<<std::endl;
0257     event.put(std::move(filteredRecHits), productLabels_[recHitCollNr]);
0258   }  //end loop over all rec hit collections
0259 }
0260
0261 template <typename RecHitType>
0262 L1RegionDataBase* HLTRecHitInAllL1RegionsProducer<RecHitType>::createL1RegionData(
0263     const std::string& type, const edm::ParameterSet& para, edm::ConsumesCollector&& consumesColl) {
0264   if (type == "L1EmParticle") {
0265     return new L1RegionData<l1extra::L1EmParticleCollection>(para, consumesColl);
0266   } else if (type == "L1JetParticle") {
0267     return new L1RegionData<l1extra::L1JetParticleCollection>(para, consumesColl);
0268   } else if (type == "L1MuonParticle") {
0269     return new L1RegionData<l1extra::L1MuonParticleCollection>(para, consumesColl);
0270   } else if (type == "EGamma") {
0271     return new L1RegionData<l1t::EGammaBxCollection>(para, consumesColl);
0272   } else if (type == "Jet") {
0273     return new L1RegionData<l1t::JetBxCollection>(para, consumesColl);
0274   } else if (type == "Muon") {
0275     return new L1RegionData<l1t::MuonBxCollection>(para, consumesColl);
0276   } else if (type == "Tau") {
0277     return new L1RegionData<l1t::TauBxCollection>(para, consumesColl);
0278   } else {
0279     //this is a major issue and could lead to rather subtle efficiency losses, so if its incorrectly configured, we're aborting the job!
0280     throw cms::Exception("InvalidConfig")
0281         << " type " << type
0282         << " is not recognised, this means the rec-hit you think you are keeping may not be and you should fix this "
0283            "error as it can lead to hard to find efficiency loses"
0284         << std::endl;
0285   }
0286 }
0287
0288 template <typename L1CollType>
0289 void L1RegionData<L1CollType>::eventSetupConsumes(edm::ConsumesCollector&) {}
0290
0291 template <typename L1CollType>
0292 void L1RegionData<L1CollType>::getEtaPhiRegions(const edm::Event& event,
0293                                                 const edm::EventSetup&,
0294                                                 std::vector<RectangularEtaPhiRegion>& regions) const {
0295   edm::Handle<L1CollType> l1Cands;
0296   event.getByToken(token_, l1Cands);
0297
0298   if (isEmpty(*l1Cands)) {
0299     LogDebug("HLTRecHitInAllL1RegionsProducerL1RegionData")
0300         << "The input collection of L1T candidates is empty (L1CollType = \""
0301         << edm::typeDemangle(typeid(L1CollType).name()) << "\"). No regions selected.";
0302     return;
0303   }
0304
0305   for (auto l1CandIt = beginIt(*l1Cands); l1CandIt != endIt(*l1Cands); ++l1CandIt) {
0306     if (l1CandIt->et() >= minEt_ && l1CandIt->et() < maxEt_) {
0307       double etaLow = l1CandIt->eta() - regionEtaMargin_;
0308       double etaHigh = l1CandIt->eta() + regionEtaMargin_;
0309       double phiLow = l1CandIt->phi() - regionPhiMargin_;
0310       double phiHigh = l1CandIt->phi() + regionPhiMargin_;
0311
0312       regions.push_back(RectangularEtaPhiRegion(etaLow, etaHigh, phiLow, phiHigh));
0313     }
0314   }
0315 }
0316
0317 template <>
0318 void L1RegionData<l1extra::L1JetParticleCollection>::eventSetupConsumes(edm::ConsumesCollector& consumesColl) {
0319   l1CaloGeometryToken_ = consumesColl.esConsumes();
0320 }
0321
0322 template <>
0323 void L1RegionData<l1extra::L1JetParticleCollection>::getEtaPhiRegions(
0324     const edm::Event& event, const edm::EventSetup& setup, std::vector<RectangularEtaPhiRegion>& regions) const {
0325   edm::Handle<l1extra::L1JetParticleCollection> l1Cands;
0326   event.getByToken(token_, l1Cands);
0327
0328   auto const& l1CaloGeom = setup.getData(l1CaloGeometryToken_);
0329
0330   for (const auto& l1Cand : *l1Cands) {
0331     if (l1Cand.et() >= minEt_ && l1Cand.et() < maxEt_) {
0332       // Access the GCT hardware object corresponding to the L1Extra EM object.
0333       int etaIndex = l1Cand.gctJetCand()->etaIndex();
0334       int phiIndex = l1Cand.gctJetCand()->phiIndex();
0335
0336       // Use the L1CaloGeometry to find the eta, phi bin boundaries.
0337       double etaLow = l1CaloGeom.etaBinLowEdge(etaIndex);
0338       double etaHigh = l1CaloGeom.etaBinHighEdge(etaIndex);
0339       double phiLow = l1CaloGeom.emJetPhiBinLowEdge(phiIndex);
0340       double phiHigh = l1CaloGeom.emJetPhiBinHighEdge(phiIndex);
0341
0342       etaLow -= regionEtaMargin_;
0343       etaHigh += regionEtaMargin_;
0344       phiLow -= regionPhiMargin_;
0345       phiHigh += regionPhiMargin_;
0346
0347       regions.push_back(RectangularEtaPhiRegion(etaLow, etaHigh, phiLow, phiHigh));
0348     }
0349   }
0350 }
0351
0352 template <>
0353 void L1RegionData<l1extra::L1EmParticleCollection>::eventSetupConsumes(edm::ConsumesCollector& consumesColl) {
0354   l1CaloGeometryToken_ = consumesColl.esConsumes();
0355 }
0356
0357 template <>
0358 void L1RegionData<l1extra::L1EmParticleCollection>::getEtaPhiRegions(
0359     const edm::Event& event, const edm::EventSetup& setup, std::vector<RectangularEtaPhiRegion>& regions) const {
0360   edm::Handle<l1extra::L1EmParticleCollection> l1Cands;
0361   event.getByToken(token_, l1Cands);
0362
0363   auto const& l1CaloGeom = setup.getData(l1CaloGeometryToken_);
0364
0365   for (const auto& l1Cand : *l1Cands) {
0366     if (l1Cand.et() >= minEt_ && l1Cand.et() < maxEt_) {
0367       // Access the GCT hardware object corresponding to the L1Extra EM object.
0368       int etaIndex = l1Cand.gctEmCand()->etaIndex();
0369       int phiIndex = l1Cand.gctEmCand()->phiIndex();
0370
0371       // Use the L1CaloGeometry to find the eta, phi bin boundaries.
0372       double etaLow = l1CaloGeom.etaBinLowEdge(etaIndex);
0373       double etaHigh = l1CaloGeom.etaBinHighEdge(etaIndex);
0374       double phiLow = l1CaloGeom.emJetPhiBinLowEdge(phiIndex);
0375       double phiHigh = l1CaloGeom.emJetPhiBinHighEdge(phiIndex);
0376
0377       etaLow -= regionEtaMargin_;
0378       etaHigh += regionEtaMargin_;
0379       phiLow -= regionPhiMargin_;
0380       phiHigh += regionPhiMargin_;
0381
0382       regions.push_back(RectangularEtaPhiRegion(etaLow, etaHigh, phiLow, phiHigh));
0383     }
0384   }
0385 }
0386
0387 typedef HLTRecHitInAllL1RegionsProducer<EcalRecHit> HLTEcalRecHitInAllL1RegionsProducer;
0388 DEFINE_FWK_MODULE(HLTEcalRecHitInAllL1RegionsProducer);
0389
0390 typedef HLTRecHitInAllL1RegionsProducer<EcalUncalibratedRecHit> HLTEcalUncalibratedRecHitInAllL1RegionsProducer;
0391 DEFINE_FWK_MODULE(HLTEcalUncalibratedRecHitInAllL1RegionsProducer);