EgammaHLTProducers/plugins/EgammaHLTExtraProducer.cc

0001
0002 #include "FWCore/Framework/interface/global/EDProducer.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0010
0011 #include "DataFormats/HLTReco/interface/EgammaObject.h"
0012 #include "DataFormats/HLTReco/interface/EgammaObjectFwd.h"
0013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0014 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0015 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0016 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
0017 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0018 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
0019 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0020 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0021 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0022 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0023 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0024 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0025 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0026 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0027 #include "DataFormats/TrackReco/interface/Track.h"
0028 #include "DataFormats/Math/interface/deltaR.h"
0029
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0032 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0033
0034 #include <sstream>
0035 #include <vector>
0036
0037 namespace {
0038   //changes double to string for product name
0039   //ie "." is replaced with "p" and for -ve vals, string is M instead so -28 is M28
0040   //has a fixed precision of precision although it removes trailing zeros and the .
0041   std::string convertToProdNameStr(double val, int precision = 3) {
0042     std::ostringstream valOStr;
0043     valOStr << std::fixed << std::setprecision(precision) << val;
0044     std::string valStr = valOStr.str();
0045     while (valStr.size() > 1 && valStr.back() == '0') {
0046       valStr.pop_back();
0047     }
0048     if (valStr.size() > 1 && valStr.back() == '.') {
0049       valStr.pop_back();
0050     }
0051     auto decPoint = valStr.find('.');
0052     if (decPoint != std::string::npos) {
0053       valStr.replace(decPoint, 1, "p");
0054     }
0055     if (val < 0)
0056       valStr.replace(0, 1, "M");
0057     return valStr;
0058   }
0059
0060   template <typename T>
0061   std::vector<std::unique_ptr<int>> countRecHits(const T& recHitHandle, const std::vector<double>& thresholds) {
0062     std::vector<std::unique_ptr<int>> counts(thresholds.size());
0063     for (auto& count : counts)
0064       count = std::make_unique<int>(0);
0065     if (recHitHandle.isValid()) {
0066       for (const auto& recHit : *recHitHandle) {
0067         for (size_t thresNr = 0; thresNr < thresholds.size(); thresNr++) {
0068           if (recHit.energy() >= thresholds[thresNr]) {
0069             (*counts[thresNr])++;
0070           }
0071         }
0072       }
0073     }
0074     return counts;
0075   }
0076 }  // namespace
0077
0078 class EgammaHLTExtraProducer : public edm::global::EDProducer<> {
0079 public:
0080   explicit EgammaHLTExtraProducer(const edm::ParameterSet& pset);
0081   ~EgammaHLTExtraProducer() override {}
0082
0083   void produce(edm::StreamID streamID, edm::Event& event, const edm::EventSetup& eventSetup) const override;
0084   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0085
0086 private:
0087   static void setVars(trigger::EgammaObject& egTrigObj,
0088                       const reco::RecoEcalCandidateRef& ecalCandRef,
0089                       const std::vector<edm::Handle<reco::RecoEcalCandidateIsolationMap>>& valueMapHandles);
0090   static reco::GsfTrackRefVector matchingGsfTrks(const reco::SuperClusterRef& scRef,
0091                                                  const edm::Handle<reco::GsfTrackCollection>& gsfTrksHandle);
0092   static void setGsfTracks(trigger::EgammaObject& egTrigObj,
0093                            const edm::Handle<reco::GsfTrackCollection>& gsfTrksHandle);
0094   static void setSeeds(trigger::EgammaObject& egTrigObj, edm::Handle<reco::ElectronSeedCollection>& eleSeedsHandle);
0095
0096   //these three filter functions are overly similar but with annoying differences
0097   //eg rechits needs to access geometry, trk dr is also w.r.t the track eta/phi
0098   //still could collapse into a single function
0099   template <typename RecHitCollection>
0100   std::unique_ptr<RecHitCollection> filterRecHits(
0101       const std::vector<std::unique_ptr<trigger::EgammaObjectCollection>>& egTrigObjs,
0102       const edm::Handle<RecHitCollection>& recHits,
0103       const CaloGeometry& geom,
0104       float maxDR2 = 0.4 * 0.4) const;
0105
0106   std::unique_ptr<reco::TrackCollection> filterTrks(
0107       const std::vector<std::unique_ptr<trigger::EgammaObjectCollection>>& egTrigObjs,
0108       const edm::Handle<reco::TrackCollection>& trks,
0109       float maxDR2 = 0.4 * 0.4) const;
0110
0111   std::unique_ptr<reco::PFClusterCollection> filterPFClusIso(
0112       const std::vector<std::unique_ptr<trigger::EgammaObjectCollection>>& egTrigObjs,
0113       const edm::Handle<reco::PFClusterCollection>& pfClus,
0114       float maxDR2 = 0.4 * 0.4) const;
0115
0116   struct Tokens {
0117     //these are the tokens which comprise the E/gamma candidate (eg cand,gsf track, pixel seeds)
0118     struct EgObjTokens {
0119       edm::EDGetTokenT<reco::RecoEcalCandidateCollection> ecalCands;
0120       edm::EDGetTokenT<reco::GsfTrackCollection> gsfTracks;
0121       edm::EDGetTokenT<reco::ElectronSeedCollection> pixelSeeds;
0122     };
0123     std::vector<std::pair<EgObjTokens, std::string>> egCands;
0124
0125     std::vector<std::pair<edm::EDGetTokenT<EcalRecHitCollection>, std::string>> ecal;
0126     std::vector<std::pair<edm::EDGetTokenT<HBHERecHitCollection>, std::string>> hcal;
0127     std::vector<std::pair<edm::EDGetTokenT<reco::TrackCollection>, std::string>> trks;
0128     std::vector<std::pair<edm::EDGetTokenT<reco::PFClusterCollection>, std::string>> pfClusIso;
0129
0130     template <typename T>
0131     static void setToken(edm::EDGetTokenT<T>& token,
0132                          edm::ConsumesCollector& cc,
0133                          const edm::ParameterSet& pset,
0134                          const std::string& tagname) {
0135       token = cc.consumes<T>(pset.getParameter<edm::InputTag>(tagname));
0136     }
0137     template <typename T>
0138     static void setToken(std::vector<edm::EDGetTokenT<T>>& tokens,
0139                          edm::ConsumesCollector& cc,
0140                          const edm::ParameterSet& pset,
0141                          const std::string& tagname) {
0142       auto inputTags = pset.getParameter<std::vector<edm::InputTag>>(tagname);
0143       tokens.resize(inputTags.size());
0144       for (size_t tagNr = 0; tagNr < inputTags.size(); tagNr++) {
0145         tokens[tagNr] = cc.consumes<T>(inputTags[tagNr]);
0146       }
0147     }
0148     template <typename T>
0149     static void setToken(std::vector<std::pair<edm::EDGetTokenT<T>, std::string>>& tokens,
0150                          edm::ConsumesCollector& cc,
0151                          const edm::ParameterSet& pset,
0152                          const std::string& tagname) {
0153       const auto& collectionPSets = pset.getParameter<std::vector<edm::ParameterSet>>(tagname);
0154       for (const auto& collPSet : collectionPSets) {
0155         edm::EDGetTokenT<T> token = cc.consumes<T>(collPSet.getParameter<edm::InputTag>("src"));
0156         std::string label = collPSet.getParameter<std::string>("label");
0157         tokens.emplace_back(token, std::move(label));
0158       }
0159     }
0160
0161     static void setToken(std::vector<std::pair<EgObjTokens, std::string>>& tokens,
0162                          edm::ConsumesCollector& cc,
0163                          const edm::ParameterSet& pset,
0164                          const std::string& tagname) {
0165       const auto& collectionPSets = pset.getParameter<std::vector<edm::ParameterSet>>(tagname);
0166       for (const auto& collPSet : collectionPSets) {
0167         EgObjTokens objTokens;
0168         setToken(objTokens.ecalCands, cc, collPSet, "ecalCands");
0169         setToken(objTokens.gsfTracks, cc, collPSet, "gsfTracks");
0170         setToken(objTokens.pixelSeeds, cc, collPSet, "pixelSeeds");
0171         std::string label = collPSet.getParameter<std::string>("label");
0172         tokens.emplace_back(objTokens, std::move(label));
0173       }
0174     }
0175     Tokens(const edm::ParameterSet& pset, edm::ConsumesCollector&& cc);
0176   };
0177
0178   const Tokens tokens_;
0179
0180   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> const caloGeomToken_;
0181
0182   float minPtToSaveHits_;
0183   bool saveHitsPlusPi_;
0184   bool saveHitsPlusHalfPi_;
0185   std::vector<double> recHitCountThresholds_;
0186 };
0187
0188 EgammaHLTExtraProducer::Tokens::Tokens(const edm::ParameterSet& pset, edm::ConsumesCollector&& cc) {
0189   setToken(egCands, cc, pset, "egCands");
0190   setToken(ecal, cc, pset, "ecal");
0191   setToken(hcal, cc, pset, "hcal");
0192   setToken(trks, cc, pset, "trks");
0193   setToken(pfClusIso, cc, pset, "pfClusIso");
0194 }
0195
0196 EgammaHLTExtraProducer::EgammaHLTExtraProducer(const edm::ParameterSet& pset)
0197     : tokens_(pset, consumesCollector()),
0198       caloGeomToken_{esConsumes()},
0199       minPtToSaveHits_(pset.getParameter<double>("minPtToSaveHits")),
0200       saveHitsPlusPi_(pset.getParameter<bool>("saveHitsPlusPi")),
0201       saveHitsPlusHalfPi_(pset.getParameter<bool>("saveHitsPlusHalfPi")),
0202       recHitCountThresholds_(pset.getParameter<std::vector<double>>("recHitCountThresholds")) {
0203   consumesMany<reco::RecoEcalCandidateIsolationMap>();
0204
0205   for (auto& tokenLabel : tokens_.egCands) {
0206     produces<trigger::EgammaObjectCollection>(tokenLabel.second);
0207   }
0208   for (auto& tokenLabel : tokens_.ecal) {
0209     produces<EcalRecHitCollection>(tokenLabel.second);
0210     for (const auto& thres : recHitCountThresholds_) {
0211       produces<int>("countEcalRecHits" + tokenLabel.second + "Thres" + convertToProdNameStr(thres) + "GeV");
0212     }
0213   }
0214   for (auto& tokenLabel : tokens_.hcal) {
0215     produces<HBHERecHitCollection>(tokenLabel.second);
0216     for (const auto& thres : recHitCountThresholds_) {
0217       produces<int>("countHcalRecHits" + tokenLabel.second + "Thres" + convertToProdNameStr(thres) + "GeV");
0218     }
0219   }
0220   for (auto& tokenLabel : tokens_.trks) {
0221     produces<reco::TrackCollection>(tokenLabel.second);
0222   }
0223   for (auto& tokenLabel : tokens_.pfClusIso) {
0224     produces<reco::PFClusterCollection>(tokenLabel.second);
0225   }
0226 }
0227
0228 void EgammaHLTExtraProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0229   edm::ParameterSetDescription desc;
0230   desc.add<double>("minPtToSaveHits", 0.);
0231   desc.add<bool>("saveHitsPlusPi", false);
0232   desc.add<bool>("saveHitsPlusHalfPi", true);
0233   desc.add<std::vector<double>>("recHitCountThresholds", std::vector{0., 0.5, 1.0, 1.5, 2.0});
0234
0235   edm::ParameterSetDescription egCandsDesc;
0236   egCandsDesc.add<edm::InputTag>("ecalCands", edm::InputTag(""));
0237   egCandsDesc.add<edm::InputTag>("pixelSeeds", edm::InputTag(""));
0238   egCandsDesc.add<edm::InputTag>("gsfTracks", edm::InputTag(""));
0239   egCandsDesc.add<std::string>("label", "");
0240   std::vector<edm::ParameterSet> egCandsDefaults(1);
0241   egCandsDefaults[0].addParameter("ecalCands", edm::InputTag("hltEgammaCandidates"));
0242   egCandsDefaults[0].addParameter("pixelSeeds", edm::InputTag("hltEgammaElectronPixelSeeds"));
0243   egCandsDefaults[0].addParameter("gsfTracks", edm::InputTag("hltEgammaGsfTracks"));
0244   egCandsDefaults[0].addParameter("label", std::string(""));
0245
0246   edm::ParameterSetDescription tokenLabelDesc;
0247   tokenLabelDesc.add<edm::InputTag>("src", edm::InputTag(""));
0248   tokenLabelDesc.add<std::string>("label", "");
0249   std::vector<edm::ParameterSet> ecalDefaults(2);
0250   ecalDefaults[0].addParameter("src", edm::InputTag("hltEcalRecHit", "EcalRecHitsEB"));
0251   ecalDefaults[0].addParameter("label", std::string("EcalRecHitsEB"));
0252   ecalDefaults[1].addParameter("src", edm::InputTag("hltEcalRecHit", "EcalRecHitsEE"));
0253   ecalDefaults[1].addParameter("label", std::string("EcalRecHitsEE"));
0254   std::vector<edm::ParameterSet> hcalDefaults(1);
0255   hcalDefaults[0].addParameter("src", edm::InputTag("hltHbhereco"));
0256   hcalDefaults[0].addParameter("label", std::string(""));
0257   std::vector<edm::ParameterSet> trksDefaults(1);
0258   trksDefaults[0].addParameter("src", edm::InputTag("generalTracks"));
0259   trksDefaults[0].addParameter("label", std::string(""));
0260   std::vector<edm::ParameterSet> pfClusIsoDefaults(3);
0261   pfClusIsoDefaults[0].addParameter("src", edm::InputTag("hltParticleFlowClusterECALL1Seeded"));
0262   pfClusIsoDefaults[0].addParameter("label", std::string("Ecal"));
0263   pfClusIsoDefaults[1].addParameter("src", edm::InputTag("hltParticleFlowClusterECALUnseeded"));
0264   pfClusIsoDefaults[1].addParameter("label", std::string("EcalUnseeded"));
0265   pfClusIsoDefaults[2].addParameter("src", edm::InputTag("hltParticleFlowClusterHCAL"));
0266   pfClusIsoDefaults[2].addParameter("label", std::string("Hcal"));
0267
0268   desc.addVPSet("egCands", egCandsDesc, egCandsDefaults);
0269   desc.addVPSet("ecal", tokenLabelDesc, ecalDefaults);
0270   desc.addVPSet("hcal", tokenLabelDesc, hcalDefaults);
0271   desc.addVPSet("trks", tokenLabelDesc, trksDefaults);
0272   desc.addVPSet("pfClusIso", tokenLabelDesc, pfClusIsoDefaults);
0273
0274   descriptions.add(("hltEgammaHLTExtraProducer"), desc);
0275 }
0276
0277 void EgammaHLTExtraProducer::produce(edm::StreamID streamID,
0278                                      edm::Event& event,
0279                                      const edm::EventSetup& eventSetup) const {
0280   std::vector<edm::Handle<reco::RecoEcalCandidateIsolationMap>> valueMapHandles;
0281   event.getManyByType(valueMapHandles);
0282
0283   std::vector<std::unique_ptr<trigger::EgammaObjectCollection>> egTrigObjColls;
0284   for (const auto& egCandsToken : tokens_.egCands) {
0285     auto ecalCandsHandle = event.getHandle(egCandsToken.first.ecalCands);
0286     auto gsfTrksHandle = event.getHandle(egCandsToken.first.gsfTracks);
0287     auto pixelSeedsHandle = event.getHandle(egCandsToken.first.pixelSeeds);
0288
0289     auto egTrigObjs = std::make_unique<trigger::EgammaObjectCollection>();
0290     for (size_t candNr = 0; ecalCandsHandle.isValid() && candNr < ecalCandsHandle->size(); candNr++) {
0291       reco::RecoEcalCandidateRef candRef(ecalCandsHandle, candNr);
0292       egTrigObjs->push_back(*candRef);
0293       auto& egTrigObj = egTrigObjs->back();
0294       setVars(egTrigObj, candRef, valueMapHandles);
0295       setGsfTracks(egTrigObj, gsfTrksHandle);
0296       setSeeds(egTrigObj, pixelSeedsHandle);
0297     }
0298     egTrigObjColls.emplace_back(std::move(egTrigObjs));
0299   }
0300
0301   auto const caloGeomHandle = eventSetup.getHandle(caloGeomToken_);
0302
0303   auto filterAndStoreRecHits = [caloGeomHandle, &event, this](const auto& egTrigObjs, const auto& tokenLabels) {
0304     for (const auto& tokenLabel : tokenLabels) {
0305       auto handle = event.getHandle(tokenLabel.first);
0306       auto recHits = filterRecHits(egTrigObjs, handle, *caloGeomHandle);
0307       event.put(std::move(recHits), tokenLabel.second);
0308     }
0309   };
0310
0311   auto storeCountRecHits = [&event](const auto& tokenLabels, const auto& thresholds, const std::string& prefixLabel) {
0312     for (const auto& tokenLabel : tokenLabels) {
0313       auto handle = event.getHandle(tokenLabel.first);
0314       auto count = countRecHits(handle, thresholds);
0315       for (size_t thresNr = 0; thresNr < thresholds.size(); thresNr++) {
0316         const auto& thres = thresholds[thresNr];
0317         event.put(std::move(count[thresNr]),
0318                   prefixLabel + tokenLabel.second + "Thres" + convertToProdNameStr(thres) + "GeV");
0319       }
0320     }
0321   };
0322
0323   auto filterAndStore = [&event, this](const auto& egTrigObjs, const auto& tokenLabels, auto filterFunc) {
0324     for (const auto& tokenLabel : tokenLabels) {
0325       auto handle = event.getHandle(tokenLabel.first);
0326       auto filtered = (this->*filterFunc)(egTrigObjs, handle, 0.4 * 0.4);
0327       event.put(std::move(filtered), tokenLabel.second);
0328     }
0329   };
0330
0331   filterAndStoreRecHits(egTrigObjColls, tokens_.ecal);
0332   filterAndStoreRecHits(egTrigObjColls, tokens_.hcal);
0333   storeCountRecHits(tokens_.ecal, recHitCountThresholds_, "countEcalRecHits");
0334   storeCountRecHits(tokens_.hcal, recHitCountThresholds_, "countHcalRecHits");
0335   filterAndStore(egTrigObjColls, tokens_.pfClusIso, &EgammaHLTExtraProducer::filterPFClusIso);
0336   filterAndStore(egTrigObjColls, tokens_.trks, &EgammaHLTExtraProducer::filterTrks);
0337
0338   for (size_t collNr = 0; collNr < egTrigObjColls.size(); collNr++) {
0339     event.put(std::move(egTrigObjColls[collNr]), tokens_.egCands[collNr].second);
0340   }
0341 }
0342
0343 void EgammaHLTExtraProducer::setVars(
0344     trigger::EgammaObject& egTrigObj,
0345     const reco::RecoEcalCandidateRef& ecalCandRef,
0346     const std::vector<edm::Handle<reco::RecoEcalCandidateIsolationMap>>& valueMapHandles) {
0347   std::vector<std::pair<std::string, float>> vars;
0348   for (auto& valueMapHandle : valueMapHandles) {
0349     auto mapIt = valueMapHandle->find(ecalCandRef);
0350     if (mapIt != valueMapHandle->end()) {
0351       std::string name = valueMapHandle.provenance()->moduleLabel();
0352       if (!valueMapHandle.provenance()->productInstanceName().empty()) {
0353         name += "_" + valueMapHandle.provenance()->productInstanceName();
0354       }
0355       vars.emplace_back(std::move(name), mapIt->val);
0356     }
0357   }
0358   egTrigObj.setVars(std::move(vars));
0359 }
0360
0361 reco::GsfTrackRefVector EgammaHLTExtraProducer::matchingGsfTrks(
0362     const reco::SuperClusterRef& scRef, const edm::Handle<reco::GsfTrackCollection>& gsfTrksHandle) {
0363   if (!gsfTrksHandle.isValid()) {
0364     return reco::GsfTrackRefVector();
0365   }
0366
0367   reco::GsfTrackRefVector gsfTrkRefs(gsfTrksHandle.id());
0368   for (size_t trkNr = 0; gsfTrksHandle.isValid() && trkNr < gsfTrksHandle->size(); trkNr++) {
0369     reco::GsfTrackRef trkRef(gsfTrksHandle, trkNr);
0370     edm::RefToBase<TrajectorySeed> seed = trkRef->extra()->seedRef();
0371     reco::ElectronSeedRef eleSeed = seed.castTo<reco::ElectronSeedRef>();
0372     edm::RefToBase<reco::CaloCluster> caloCluster = eleSeed->caloCluster();
0373     reco::SuperClusterRef scRefFromTrk = caloCluster.castTo<reco::SuperClusterRef>();
0374     if (scRefFromTrk == scRef) {
0375       gsfTrkRefs.push_back(trkRef);
0376     }
0377   }
0378   return gsfTrkRefs;
0379 }
0380
0381 void EgammaHLTExtraProducer::setGsfTracks(trigger::EgammaObject& egTrigObj,
0382                                           const edm::Handle<reco::GsfTrackCollection>& gsfTrksHandle) {
0383   egTrigObj.setGsfTracks(matchingGsfTrks(egTrigObj.superCluster(), gsfTrksHandle));
0384 }
0385
0386 void EgammaHLTExtraProducer::setSeeds(trigger::EgammaObject& egTrigObj,
0387                                       edm::Handle<reco::ElectronSeedCollection>& eleSeedsHandle) {
0388   if (!eleSeedsHandle.isValid()) {
0389     egTrigObj.setSeeds(reco::ElectronSeedRefVector());
0390   } else {
0391     reco::ElectronSeedRefVector trigObjSeeds(eleSeedsHandle.id());
0392
0393     for (size_t seedNr = 0; eleSeedsHandle.isValid() && seedNr < eleSeedsHandle->size(); seedNr++) {
0394       reco::ElectronSeedRef eleSeed(eleSeedsHandle, seedNr);
0395       edm::RefToBase<reco::CaloCluster> caloCluster = eleSeed->caloCluster();
0396       reco::SuperClusterRef scRefFromSeed = caloCluster.castTo<reco::SuperClusterRef>();
0397
0398       if (scRefFromSeed == egTrigObj.superCluster()) {
0399         trigObjSeeds.push_back(eleSeed);
0400       }
0401     }
0402     egTrigObj.setSeeds(std::move(trigObjSeeds));
0403   }
0404 }
0405
0406 template <typename RecHitCollection>
0407 std::unique_ptr<RecHitCollection> EgammaHLTExtraProducer::filterRecHits(
0408     const std::vector<std::unique_ptr<trigger::EgammaObjectCollection>>& egTrigObjColls,
0409     const edm::Handle<RecHitCollection>& recHits,
0410     const CaloGeometry& geom,
0411     float maxDR2) const {
0412   auto filteredHits = std::make_unique<RecHitCollection>();
0413   if (!recHits.isValid())
0414     return filteredHits;
0415
0416   std::vector<std::pair<float, float>> etaPhis;
0417   for (const auto& egTrigObjs : egTrigObjColls) {
0418     for (const auto& egTrigObj : *egTrigObjs) {
0419       if (egTrigObj.pt() >= minPtToSaveHits_) {
0420         etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi()});
0421         if (saveHitsPlusPi_)
0422           etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi() + 3.14159});
0423         if (saveHitsPlusHalfPi_)
0424           etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi() + 3.14159 / 2.});
0425       }
0426     }
0427   }
0428   auto deltaR2Match = [&etaPhis, &maxDR2](const GlobalPoint& pos) {
0429     float eta = pos.eta();
0430     float phi = pos.phi();
0431     for (auto& etaPhi : etaPhis) {
0432       if (reco::deltaR2(eta, phi, etaPhi.first, etaPhi.second) < maxDR2)
0433         return true;
0434     }
0435     return false;
0436   };
0437
0438   for (auto& hit : *recHits) {
0439     const CaloSubdetectorGeometry* subDetGeom = geom.getSubdetectorGeometry(hit.id());
0440     if (subDetGeom) {
0441       auto cellGeom = subDetGeom->getGeometry(hit.id());
0442       if (deltaR2Match(cellGeom->getPosition()))
0443         filteredHits->push_back(hit);
0444     } else {
0445       throw cms::Exception("GeomError") << "could not get geometry for det id " << hit.id().rawId();
0446     }
0447   }
0448   return filteredHits;
0449 }
0450
0451 std::unique_ptr<reco::TrackCollection> EgammaHLTExtraProducer::filterTrks(
0452     const std::vector<std::unique_ptr<trigger::EgammaObjectCollection>>& egTrigObjColls,
0453     const edm::Handle<reco::TrackCollection>& trks,
0454     float maxDR2) const {
0455   auto filteredTrks = std::make_unique<reco::TrackCollection>();
0456   if (!trks.isValid())
0457     return filteredTrks;
0458
0459   //so because each egamma object can have multiple eta/phi pairs
0460   //easier to just make a temp vector and then copy that in with the +pi and  +pi/2
0461   std::vector<std::pair<float, float>> etaPhisTmp;
0462   for (const auto& egTrigObjs : egTrigObjColls) {
0463     for (const auto& egTrigObj : *egTrigObjs) {
0464       if (egTrigObj.pt() >= minPtToSaveHits_) {
0465         etaPhisTmp.push_back({egTrigObj.eta(), egTrigObj.phi()});
0466         //also save the eta /phi of all gsf tracks with the object
0467         for (const auto& gsfTrk : egTrigObj.gsfTracks()) {
0468           etaPhisTmp.push_back({gsfTrk->eta(), gsfTrk->phi()});
0469         }
0470       }
0471     }
0472   }
0473   std::vector<std::pair<float, float>> etaPhis;
0474   for (const auto& etaPhi : etaPhisTmp) {
0475     etaPhis.push_back(etaPhi);
0476     if (saveHitsPlusPi_)
0477       etaPhis.push_back({etaPhi.first, etaPhi.second + 3.14159});
0478     if (saveHitsPlusHalfPi_)
0479       etaPhis.push_back({etaPhi.first, etaPhi.second + 3.14159 / 2.});
0480   }
0481
0482   auto deltaR2Match = [&etaPhis, &maxDR2](float eta, float phi) {
0483     for (auto& etaPhi : etaPhis) {
0484       if (reco::deltaR2(eta, phi, etaPhi.first, etaPhi.second) < maxDR2)
0485         return true;
0486     }
0487     return false;
0488   };
0489
0490   for (auto& trk : *trks) {
0491     if (deltaR2Match(trk.eta(), trk.phi()))
0492       filteredTrks->push_back(trk);
0493   }
0494   return filteredTrks;
0495 }
0496
0497 std::unique_ptr<reco::PFClusterCollection> EgammaHLTExtraProducer::filterPFClusIso(
0498     const std::vector<std::unique_ptr<trigger::EgammaObjectCollection>>& egTrigObjColls,
0499     const edm::Handle<reco::PFClusterCollection>& pfClus,
0500     float maxDR2) const {
0501   auto filteredPFClus = std::make_unique<reco::PFClusterCollection>();
0502   if (!pfClus.isValid())
0503     return filteredPFClus;
0504
0505   std::vector<std::pair<float, float>> etaPhis;
0506   for (const auto& egTrigObjs : egTrigObjColls) {
0507     for (const auto& egTrigObj : *egTrigObjs) {
0508       if (egTrigObj.pt() >= minPtToSaveHits_) {
0509         etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi()});
0510         if (saveHitsPlusPi_)
0511           etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi() + 3.14159});
0512         if (saveHitsPlusHalfPi_)
0513           etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi() + 3.14159 / 2.});
0514       }
0515     }
0516   }
0517   auto deltaR2Match = [&etaPhis, &maxDR2](float eta, float phi) {
0518     for (auto& etaPhi : etaPhis) {
0519       if (reco::deltaR2(eta, phi, etaPhi.first, etaPhi.second) < maxDR2)
0520         return true;
0521     }
0522     return false;
0523   };
0524
0525   for (auto& clus : *pfClus) {
0526     if (deltaR2Match(clus.eta(), clus.phi()))
0527       filteredPFClus->push_back(clus);
0528   }
0529   return filteredPFClus;
0530 }
0531
0532 #include "FWCore/Framework/interface/MakerMacros.h"
0533 DEFINE_FWK_MODULE(EgammaHLTExtraProducer);