EgammaHLTProducers/plugins/HLTCaloObjInRegionsProducer.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/ConsumesCollector.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0015 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0016 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
0017 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0018 #include "DataFormats/L1Trigger/interface/P2GTCandidate.h"
0019 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0020 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
0021 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0022 #include "DataFormats/DetId/interface/DetId.h"
0023
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
0029 #include "DataFormats/L1Trigger/interface/EGamma.h"
0030 #include "DataFormats/L1Trigger/interface/Jet.h"
0031 #include "DataFormats/L1Trigger/interface/Muon.h"
0032
0033 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0034 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0035
0036 /**************************************************************
0037 / purpose: enable filtering of calo objects in eta/phi or deltaR
0038 /          regions around generic objects
0039 /
0040 / operation : accepts all objects with
0041 /             (dEta <dEtaMax  && dPhi < dPhiMax) || dR < dRMax
0042 /             so the OR of a rectangular region and cone region
0043 ****************************************************************/
0044
0045 //this is a struct which contains all the eta/phi regions
0046 //from which to filter the calo objs
0047 class EtaPhiRegion {
0048 private:
0049   float centreEta_;
0050   float centrePhi_;
0051   float maxDeltaR2_;
0052   float maxDEta_;
0053   float maxDPhi_;
0054
0055 public:
0056   EtaPhiRegion(float iEta, float iPhi, float iDR, float iDEta, float iDPhi)
0057       : centreEta_(iEta), centrePhi_(iPhi), maxDeltaR2_(iDR * iDR), maxDEta_(iDEta), maxDPhi_(iDPhi) {}
0058   ~EtaPhiRegion() {}
0059   bool operator()(float eta, float phi) const {
0060     return reco::deltaR2(eta, phi, centreEta_, centrePhi_) < maxDeltaR2_ ||
0061            (std::abs(eta - centreEta_) < maxDEta_ && std::abs(reco::deltaPhi(phi, centrePhi_)) < maxDPhi_);
0062   }
0063 };
0064
0065 class EtaPhiRegionDataBase {
0066 public:
0067   EtaPhiRegionDataBase() {}
0068   virtual ~EtaPhiRegionDataBase() = default;
0069   virtual void getEtaPhiRegions(const edm::Event&, std::vector<EtaPhiRegion>&) const = 0;
0070 };
0071
0072 //this class stores the tokens to access the objects around which we wish to filter
0073 //it makes a vector of EtaPhiRegions which are then used to filter the CaloObjs
0074 template <typename T1>
0075 class EtaPhiRegionData : public EtaPhiRegionDataBase {
0076 private:
0077   float minEt_;
0078   float maxEt_;
0079   float maxDeltaR_;
0080   float maxDEta_;
0081   float maxDPhi_;
0082   edm::EDGetTokenT<T1> token_;
0083
0084 public:
0085   EtaPhiRegionData(const edm::ParameterSet& para, edm::ConsumesCollector& consumesColl)
0086       : minEt_(para.getParameter<double>("minEt")),
0087         maxEt_(para.getParameter<double>("maxEt")),
0088         maxDeltaR_(para.getParameter<double>("maxDeltaR")),
0089         maxDEta_(para.getParameter<double>("maxDEta")),
0090         maxDPhi_(para.getParameter<double>("maxDPhi")),
0091         token_(consumesColl.consumes<T1>(para.getParameter<edm::InputTag>("inputColl"))) {}
0092
0093   void getEtaPhiRegions(const edm::Event&, std::vector<EtaPhiRegion>&) const override;
0094 };
0095
0096 template <typename CaloObjType, typename CaloObjCollType = edm::SortedCollection<CaloObjType>>
0097 class HLTCaloObjInRegionsProducer : public edm::stream::EDProducer<> {
0098 public:
0099   HLTCaloObjInRegionsProducer(const edm::ParameterSet& ps);
0100   ~HLTCaloObjInRegionsProducer() override {}
0101
0102   void produce(edm::Event&, const edm::EventSetup&) override;
0103   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0104
0105 private:
0106   EtaPhiRegionDataBase* createEtaPhiRegionData(const std::string&,
0107                                                const edm::ParameterSet&,
0108                                                edm::ConsumesCollector&&);  //calling function owns this
0109   static std::unique_ptr<CaloObjCollType> makeFilteredColl(const edm::Handle<CaloObjCollType>& inputColl,
0110                                                            CaloGeometry const& caloGeomHandle,
0111                                                            const std::vector<EtaPhiRegion>& regions);
0112   static bool validIDForGeom(const DetId& id);
0113   std::vector<std::string> outputProductNames_;
0114   std::vector<edm::InputTag> inputCollTags_;
0115   std::vector<edm::EDGetTokenT<CaloObjCollType>> inputTokens_;
0116   std::vector<std::unique_ptr<EtaPhiRegionDataBase>> etaPhiRegionData_;
0117   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0118 };
0119
0120 template <typename CaloObjType, typename CaloObjCollType>
0121 HLTCaloObjInRegionsProducer<CaloObjType, CaloObjCollType>::HLTCaloObjInRegionsProducer(const edm::ParameterSet& para)
0122     : caloGeometryToken_{esConsumes()} {
0123   const std::vector<edm::ParameterSet> etaPhiRegions =
0124       para.getParameter<std::vector<edm::ParameterSet>>("etaPhiRegions");
0125   for (auto& pset : etaPhiRegions) {
0126     const std::string type = pset.getParameter<std::string>("type");
0127     etaPhiRegionData_.emplace_back(createEtaPhiRegionData(
0128         type,
0129         pset,
0130         consumesCollector()));  //meh I was going to use a factory but it was going to be overly complex for my needs
0131   }
0132
0133   outputProductNames_ = para.getParameter<std::vector<std::string>>("outputProductNames");
0134   inputCollTags_ = para.getParameter<std::vector<edm::InputTag>>("inputCollTags");
0135   if (outputProductNames_.size() != inputCollTags_.size()) {
0136     throw cms::Exception("InvalidConfiguration")
0137         << " error outputProductNames and inputCollTags must be the same size, they are " << outputProductNames_.size()
0138         << " vs " << inputCollTags_.size();
0139   }
0140   for (unsigned int collNr = 0; collNr < inputCollTags_.size(); collNr++) {
0141     inputTokens_.push_back(consumes<CaloObjCollType>(inputCollTags_[collNr]));
0142     produces<CaloObjCollType>(outputProductNames_[collNr]);
0143   }
0144 }
0145
0146 template <typename CaloObjType, typename CaloObjCollType>
0147 void HLTCaloObjInRegionsProducer<CaloObjType, CaloObjCollType>::fillDescriptions(
0148     edm::ConfigurationDescriptions& descriptions) {
0149   edm::ParameterSetDescription desc;
0150   std::vector<std::string> outputProductNames;
0151   outputProductNames.push_back("EcalRegionalRecHitsEB");
0152   desc.add<std::vector<std::string>>("outputProductNames", outputProductNames);
0153   std::vector<edm::InputTag> inputColls;
0154   inputColls.push_back(edm::InputTag("hltHcalDigis"));
0155   desc.add<std::vector<edm::InputTag>>("inputCollTags", inputColls);
0156   std::vector<edm::ParameterSet> etaPhiRegions;
0157
0158   edm::ParameterSet ecalCandPSet;
0159   ecalCandPSet.addParameter<std::string>("type", "RecoEcalCandidate");
0160   ecalCandPSet.addParameter<double>("minEt", -1);
0161   ecalCandPSet.addParameter<double>("maxEt", -1);
0162   ecalCandPSet.addParameter<double>("maxDeltaR", 0.5);
0163   ecalCandPSet.addParameter<double>("maxDEta", 0.);
0164   ecalCandPSet.addParameter<double>("maxDPhi", 0.);
0165   ecalCandPSet.addParameter<edm::InputTag>("inputColl", edm::InputTag("hltEgammaCandidates"));
0166   etaPhiRegions.push_back(ecalCandPSet);
0167
0168   edm::ParameterSetDescription etaPhiRegionDesc;
0169   etaPhiRegionDesc.add<std::string>("type");
0170   etaPhiRegionDesc.add<double>("minEt");
0171   etaPhiRegionDesc.add<double>("maxEt");
0172   etaPhiRegionDesc.add<double>("maxDeltaR");
0173   etaPhiRegionDesc.add<double>("maxDEta");
0174   etaPhiRegionDesc.add<double>("maxDPhi");
0175   etaPhiRegionDesc.add<edm::InputTag>("inputColl");
0176   desc.addVPSet("etaPhiRegions", etaPhiRegionDesc, etaPhiRegions);
0177
0178   descriptions.add(defaultModuleLabel<HLTCaloObjInRegionsProducer<CaloObjType, CaloObjCollType>>(), desc);
0179 }
0180
0181 template <typename CaloObjType, typename CaloObjCollType>
0182 void HLTCaloObjInRegionsProducer<CaloObjType, CaloObjCollType>::produce(edm::Event& event,
0183                                                                         const edm::EventSetup& setup) {
0184   // get the collection geometry:
0185   auto const& caloGeom = setup.getData(caloGeometryToken_);
0186
0187   std::vector<EtaPhiRegion> regions;
0188   std::for_each(etaPhiRegionData_.begin(),
0189                 etaPhiRegionData_.end(),
0190                 [&event, &regions](const std::unique_ptr<EtaPhiRegionDataBase>& input) {
0191                   input->getEtaPhiRegions(event, regions);
0192                 });
0193
0194   for (size_t inputCollNr = 0; inputCollNr < inputTokens_.size(); inputCollNr++) {
0195     edm::Handle<CaloObjCollType> inputColl;
0196     event.getByToken(inputTokens_[inputCollNr], inputColl);
0197
0198     if (!(inputColl.isValid())) {
0199       edm::LogError("ProductNotFound") << "could not get a handle on the " << typeid(CaloObjCollType).name()
0200                                        << " named " << inputCollTags_[inputCollNr].encode() << std::endl;
0201       continue;
0202     }
0203     auto outputColl = makeFilteredColl(inputColl, caloGeom, regions);
0204     event.put(std::move(outputColl), outputProductNames_[inputCollNr]);
0205   }
0206 }
0207
0208 template <typename CaloObjType, typename CaloObjCollType>
0209 std::unique_ptr<CaloObjCollType> HLTCaloObjInRegionsProducer<CaloObjType, CaloObjCollType>::makeFilteredColl(
0210     const edm::Handle<CaloObjCollType>& inputColl,
0211     CaloGeometry const& caloGeomHandle,
0212     const std::vector<EtaPhiRegion>& regions) {
0213   auto outputColl = std::make_unique<CaloObjCollType>();
0214   if (!inputColl->empty()) {
0215     const CaloSubdetectorGeometry* subDetGeom = caloGeomHandle.getSubdetectorGeometry(inputColl->begin()->id());
0216     if (!regions.empty()) {
0217       for (auto const& obj : *inputColl) {
0218         auto objGeom = subDetGeom->getGeometry(obj.id());
0219         if (objGeom == nullptr) {
0220           //wondering what to do here
0221           //something is very very wrong
0222           //given HLT should never crash or throw, decided to log an error
0223           //update: so turns out HCAL can pass through calibration channels in QIE11 so for that module, its an expected behaviour
0224           //so we check if the ID is valid
0225           if (validIDForGeom(obj.id())) {
0226             edm::LogError("HLTCaloObjInRegionsProducer")
0227                 << "for an object of type " << typeid(CaloObjType).name() << " the geometry returned null for id "
0228                 << DetId(obj.id()).rawId() << " with initial ID " << DetId(inputColl->begin()->id()).rawId()
0229                 << " in HLTCaloObjsInRegion, this shouldnt be possible and something has gone wrong, auto accepting "
0230                    "hit";
0231           }
0232           outputColl->push_back(obj);
0233           continue;
0234         }
0235         float eta = objGeom->getPosition().eta();
0236         float phi = objGeom->getPosition().phi();
0237
0238         for (const auto& region : regions) {
0239           if (region(eta, phi)) {
0240             outputColl->push_back(obj);
0241             break;
0242           }
0243         }
0244       }
0245     }  //end check of empty regions
0246   }  //end check of empty rec-hits
0247   return outputColl;
0248 }
0249
0250 //tells us if an ID should have a valid geometry
0251 //it assumes that all IDs do except those specifically mentioned
0252 //HCAL for example have laser calibs in the digi collection so
0253 //so we have to ensure that HCAL is HB,HE or HO
0254 template <typename CaloObjType, typename CaloObjCollType>
0255 bool HLTCaloObjInRegionsProducer<CaloObjType, CaloObjCollType>::validIDForGeom(const DetId& id) {
0256   if (id.det() == DetId::Hcal) {
0257     if (id.subdetId() == HcalSubdetector::HcalEmpty || id.subdetId() == HcalSubdetector::HcalOther) {
0258       return false;
0259     }
0260   }
0261   return true;
0262 }
0263
0264 template <typename CaloObjType, typename CaloObjCollType>
0265 EtaPhiRegionDataBase* HLTCaloObjInRegionsProducer<CaloObjType, CaloObjCollType>::createEtaPhiRegionData(
0266     const std::string& type, const edm::ParameterSet& para, edm::ConsumesCollector&& consumesColl) {
0267   if (type == "L1P2GTCandidate") {
0268     return new EtaPhiRegionData<l1t::P2GTCandidateCollection>(para, consumesColl);
0269   } else if (type == "L1EGamma") {
0270     return new EtaPhiRegionData<l1t::EGammaBxCollection>(para, consumesColl);
0271   } else if (type == "L1Jet") {
0272     return new EtaPhiRegionData<l1t::JetBxCollection>(para, consumesColl);
0273   } else if (type == "L1Muon") {
0274     return new EtaPhiRegionData<l1t::MuonBxCollection>(para, consumesColl);
0275   } else if (type == "L1Tau") {
0276     return new EtaPhiRegionData<l1t::TauBxCollection>(para, consumesColl);
0277   } else if (type == "RecoEcalCandidate") {
0278     return new EtaPhiRegionData<reco::RecoEcalCandidateCollection>(para, consumesColl);
0279   } else if (type == "RecoChargedCandidate") {
0280     return new EtaPhiRegionData<reco::RecoChargedCandidateCollection>(para, consumesColl);
0281   } else if (type == "Electron") {
0282     return new EtaPhiRegionData<reco::Electron>(para, consumesColl);
0283   } else {
0284     //this is a major issue and could lead to rather subtle efficiency losses, so if its incorrectly configured, we're aborting the job!
0285     throw cms::Exception("InvalidConfig")
0286         << " type " << type
0287         << " is not recognised, this means the rec-hit you think you are keeping may not be and you should fix this "
0288            "error as it can lead to hard to find efficiency loses"
0289         << std::endl;
0290   }
0291 }
0292
0293 template <typename CandCollType>
0294 void EtaPhiRegionData<CandCollType>::getEtaPhiRegions(const edm::Event& event,
0295                                                       std::vector<EtaPhiRegion>& regions) const {
0296   edm::Handle<CandCollType> cands;
0297   event.getByToken(token_, cands);
0298
0299   for (auto const& cand : *cands) {
0300     if (cand.et() >= minEt_ && (maxEt_ < 0 || cand.et() < maxEt_)) {
0301       regions.push_back(EtaPhiRegion(cand.eta(), cand.phi(), maxDeltaR_, maxDEta_, maxDPhi_));
0302     }
0303   }
0304 }
0305
0306 #include "FWCore/Framework/interface/MakerMacros.h"
0307
0308 #include "DataFormats/HcalDigi/interface/HBHEDataFrame.h"
0309 using HLTHcalDigisInRegionsProducer = HLTCaloObjInRegionsProducer<HBHEDataFrame>;
0310 DEFINE_FWK_MODULE(HLTHcalDigisInRegionsProducer);
0311
0312 #include "DataFormats/HcalDigi/interface/QIE11DataFrame.h"
0313 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0314 using HLTHcalQIE11DigisInRegionsProducer = HLTCaloObjInRegionsProducer<QIE11DataFrame, QIE11DigiCollection>;
0315 DEFINE_FWK_MODULE(HLTHcalQIE11DigisInRegionsProducer);
0316
0317 #include "DataFormats/HcalDigi/interface/QIE10DataFrame.h"
0318 using HLTHcalQIE10DigisInRegionsProducer = HLTCaloObjInRegionsProducer<QIE10DataFrame, QIE10DigiCollection>;
0319 DEFINE_FWK_MODULE(HLTHcalQIE10DigisInRegionsProducer);
0320
0321 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0322 using HLTEcalEBDigisInRegionsProducer = HLTCaloObjInRegionsProducer<EBDataFrame, EBDigiCollection>;
0323 DEFINE_FWK_MODULE(HLTEcalEBDigisInRegionsProducer);
0324 using HLTEcalEEDigisInRegionsProducer = HLTCaloObjInRegionsProducer<EEDataFrame, EEDigiCollection>;
0325 DEFINE_FWK_MODULE(HLTEcalEEDigisInRegionsProducer);
0326 using HLTEcalESDigisInRegionsProducer = HLTCaloObjInRegionsProducer<ESDataFrame, ESDigiCollection>;
0327 DEFINE_FWK_MODULE(HLTEcalESDigisInRegionsProducer);
0328
0329 //these two classes are intended to ultimately replace the EcalRecHit and EcalUncalibratedRecHit
0330 //instances of HLTRecHitInAllL1RegionsProducer, particulary as we're free of legacy / stage-1 L1 now
0331 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0332 using HLTEcalRecHitsInRegionsProducer = HLTCaloObjInRegionsProducer<EcalRecHit>;
0333 DEFINE_FWK_MODULE(HLTEcalRecHitsInRegionsProducer);
0334 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0335 using HLTEcalUnCalibRecHitsInRegionsProducer = HLTCaloObjInRegionsProducer<EcalUncalibratedRecHit>;
0336 DEFINE_FWK_MODULE(HLTEcalUnCalibRecHitsInRegionsProducer);
0337
0338 // HGCAL Digis
0339 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0340 using HLTHGCalDigisInRegionsProducer = HLTCaloObjInRegionsProducer<HGCalDataFrame, HGCalDigiCollection>;
0341 DEFINE_FWK_MODULE(HLTHGCalDigisInRegionsProducer);
0342
0343 // HGCAL RecHits
0344 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0345 using HLTHGCalRecHitsInRegionsProducer = HLTCaloObjInRegionsProducer<HGCRecHit, std::vector<HGCRecHit>>;
0346 DEFINE_FWK_MODULE(HLTHGCalRecHitsInRegionsProducer);
0347 #include "DataFormats/HGCRecHit/interface/HGCUncalibratedRecHit.h"
0348 using HLTHGCalUncalibratedRecHitsInRegionsProducer = HLTCaloObjInRegionsProducer<HGCUncalibratedRecHit>;
0349 DEFINE_FWK_MODULE(HLTHGCalUncalibratedRecHitsInRegionsProducer);