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