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