EgammaHLTProducers/plugins/EgammaHLTPhase2ExtraProducer.cc

0001 #include "FWCore/Framework/interface/global/EDProducer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0008
0009 #include "DataFormats/HLTReco/interface/EgammaObject.h"
0010 #include "DataFormats/HLTReco/interface/EgammaObjectFwd.h"
0011 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0012 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0014 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0015 #include "DataFormats/L1TrackTrigger/interface/L1Track.h"
0016 #include "DataFormats/Math/interface/deltaR.h"
0017 #include "DataFormats/Common/interface/RefToPtr.h"
0018 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0019 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0020 #include "DataFormats/Common/interface/ValueMap.h"
0021
0022 #include "SimDataFormats/Associations/interface/L1TrackTruthPair.h"
0023 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0024 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0025 #include "SimTracker/TrackTriggerAssociation/interface/TTTrackAssociationMap.h"
0026
0027 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0028 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0029 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0030
0031 #include <vector>
0032 #include <unordered_map>
0033
0034 //class complements EgammaHLTExtraProducer and adds all the phase-II specific E/g HLT debug information to the event
0035 //this allows phase-II to be factorised from the standard class rather than having to extend EgammaHLTExtraProducer to deal with it
0036 //although to be fair, given all the phase-II code is now in the release, the need for this factorisation is not as great
0037 //and ultimately it could be merged into EgammaHLTExtraProducer
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 EgammaHLTPhase2ExtraProducer : public edm::global::EDProducer<> {
0081 public:
0082   explicit EgammaHLTPhase2ExtraProducer(const edm::ParameterSet& pset);
0083
0084   void produce(edm::StreamID streamID, edm::Event& event, const edm::EventSetup& eventSetup) const override;
0085   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0086
0087 private:
0088   template <typename CollType, typename RefType>
0089   std::unique_ptr<CollType> filterObjs(const trigger::EgammaObjectCollection& egTrigObjs,
0090                                        const edm::Handle<CollType>& objs,
0091                                        std::vector<RefType>& orgRefs,
0092                                        float maxDR2 = 0.4 * 0.4) const;
0093
0094   //these three filter functions are overly similar but with annoying differences
0095   //eg rechits needs to access geometry, trk dr is also w.r.t the track eta/phi
0096   //still could collapse into a single function
0097   template <typename RecHitCollection>
0098   std::unique_ptr<RecHitCollection> filterRecHits(const trigger::EgammaObjectCollection& egTrigObjs,
0099                                                   const edm::Handle<RecHitCollection>& recHits,
0100                                                   const CaloGeometry& geom,
0101                                                   float maxDR2 = 0.4 * 0.4) const;
0102
0103   struct Tokens {
0104     edm::EDGetTokenT<trigger::EgammaObjectCollection> egTrigObjs;
0105     edm::EDGetTokenT<L1TrackCollection> l1Trks;
0106     edm::EDGetTokenT<TrackingParticleCollection> trkParts;
0107     edm::EDGetTokenT<TTTrackAssociationMap<Ref_Phase2TrackerDigi_>> l1TrkToTrkPartMap;
0108     edm::EDGetTokenT<reco::CaloClusterCollection> hgcalLayerClusters;
0109     edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> hgcalLayerClustersTime;
0110     std::vector<std::pair<edm::EDGetTokenT<HGCRecHitCollection>, std::string>> hgcal;
0111
0112     template <typename T>
0113     static void setToken(edm::EDGetTokenT<T>& token,
0114                          edm::ConsumesCollector& cc,
0115                          const edm::ParameterSet& pset,
0116                          const std::string& tagname) {
0117       token = cc.consumes<T>(pset.getParameter<edm::InputTag>(tagname));
0118     }
0119     template <typename T>
0120     static void setToken(std::vector<edm::EDGetTokenT<T>>& tokens,
0121                          edm::ConsumesCollector& cc,
0122                          const edm::ParameterSet& pset,
0123                          const std::string& tagname) {
0124       auto inputTags = pset.getParameter<std::vector<edm::InputTag>>(tagname);
0125       tokens.resize(inputTags.size());
0126       for (size_t tagNr = 0; tagNr < inputTags.size(); tagNr++) {
0127         tokens[tagNr] = cc.consumes<T>(inputTags[tagNr]);
0128       }
0129     }
0130     template <typename T>
0131     static void setToken(std::vector<std::pair<edm::EDGetTokenT<T>, std::string>>& tokens,
0132                          edm::ConsumesCollector& cc,
0133                          const edm::ParameterSet& pset,
0134                          const std::string& tagname) {
0135       const auto& collectionPSets = pset.getParameter<std::vector<edm::ParameterSet>>(tagname);
0136       for (const auto& collPSet : collectionPSets) {
0137         edm::EDGetTokenT<T> token = cc.consumes<T>(collPSet.getParameter<edm::InputTag>("src"));
0138         std::string label = collPSet.getParameter<std::string>("label");
0139         tokens.emplace_back(std::make_pair(token, label));
0140       }
0141     }
0142     Tokens(const edm::ParameterSet& pset, edm::ConsumesCollector&& cc);
0143   };
0144   template <typename T, typename H>
0145   static std::unique_ptr<edm::ValueMap<T>> makeValueMap(const H& handle, const std::vector<T>& values) {
0146     auto valueMap = std::make_unique<edm::ValueMap<T>>();
0147     typename edm::ValueMap<T>::Filler filler(*valueMap);
0148     filler.insert(handle, values.begin(), values.end());
0149     filler.fill();
0150     return valueMap;
0151   }
0152
0153   const Tokens tokens_;
0154
0155   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> const caloGeomToken_;
0156
0157   float minPtToSaveHits_;
0158   bool saveHitsPlusPi_;
0159   bool saveHitsPlusHalfPi_;
0160   std::vector<double> recHitCountThresholds_;
0161 };
0162
0163 EgammaHLTPhase2ExtraProducer::Tokens::Tokens(const edm::ParameterSet& pset, edm::ConsumesCollector&& cc) {
0164   setToken(egTrigObjs, cc, pset, "egTrigObjs");
0165   setToken(l1Trks, cc, pset, "l1Trks");
0166   setToken(trkParts, cc, pset, "trkParts");
0167   setToken(l1TrkToTrkPartMap, cc, pset, "l1TrkToTrkPartMap");
0168   setToken(hgcalLayerClusters, cc, pset, "hgcalLayerClusters");
0169   setToken(hgcalLayerClustersTime, cc, pset, "hgcalLayerClustersTime");
0170   setToken(hgcal, cc, pset, "hgcal");
0171 }
0172
0173 EgammaHLTPhase2ExtraProducer::EgammaHLTPhase2ExtraProducer(const edm::ParameterSet& pset)
0174     : tokens_(pset, consumesCollector()),
0175       caloGeomToken_{esConsumes()},
0176       minPtToSaveHits_(pset.getParameter<double>("minPtToSaveHits")),
0177       saveHitsPlusPi_(pset.getParameter<bool>("saveHitsPlusPi")),
0178       saveHitsPlusHalfPi_(pset.getParameter<bool>("saveHitsPlusHalfPi")),
0179       recHitCountThresholds_(pset.getParameter<std::vector<double>>("recHitCountThresholds")) {
0180   produces<L1TrackCollection>();
0181   produces<L1TrackTruthPairCollection>();
0182   produces<TrackingParticleCollection>();
0183   produces<reco::CaloClusterCollection>("hgcalLayerClusters");
0184   produces<edm::ValueMap<std::pair<float, float>>>("hgcalLayerClustersTime");
0185   for (auto& tokenLabel : tokens_.hgcal) {
0186     produces<HGCRecHitCollection>(tokenLabel.second);
0187     for (const auto& thres : recHitCountThresholds_) {
0188       produces<int>("countHgcalRecHits" + tokenLabel.second + "Thres" + convertToProdNameStr(thres) + "GeV");
0189     }
0190   }
0191 }
0192
0193 void EgammaHLTPhase2ExtraProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0194   edm::ParameterSetDescription desc;
0195   desc.add<edm::InputTag>("egTrigObjs", edm::InputTag("hltEgammaHLTExtra"));
0196   desc.add<edm::InputTag>("l1Trks", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks"));
0197   desc.add<edm::InputTag>("trkParts", edm::InputTag("mix", "MergedTrackTruth"));
0198   desc.add<edm::InputTag>("l1TrkToTrkPartMap", edm::InputTag("TTTrackAssociatorFromPixelDigis", "Level1TTTracks"));
0199   desc.add<edm::InputTag>("hgcalLayerClusters", edm::InputTag("hgcalLayerClusters"));
0200   desc.add<edm::InputTag>("hgcalLayerClustersTime", edm::InputTag("hgcalLayerClusters", "timeLayerCluster"));
0201   desc.add<double>("minPtToSaveHits", 0.);
0202   desc.add<bool>("saveHitsPlusPi", true);
0203   desc.add<bool>("saveHitsPlusHalfPi", true);
0204   desc.add<std::vector<double>>("recHitCountThresholds", std::vector{0., 0.5, 1.0, 1.5, 2.0});
0205   std::vector<edm::ParameterSet> ecalDefaults(2);
0206   edm::ParameterSetDescription tokenLabelDesc;
0207   tokenLabelDesc.add<edm::InputTag>("src", edm::InputTag(""));
0208   tokenLabelDesc.add<std::string>("label", "");
0209   std::vector<edm::ParameterSet> hgcalDefaults(3);
0210   hgcalDefaults[0].addParameter("src", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0211   hgcalDefaults[0].addParameter("label", std::string("HGCEERecHits"));
0212   hgcalDefaults[1].addParameter("src", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0213   hgcalDefaults[1].addParameter("label", std::string("HGCHEFRecHits"));
0214   hgcalDefaults[2].addParameter("src", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0215   hgcalDefaults[2].addParameter("label", std::string("HGCHEBRecHits"));
0216   desc.addVPSet("hgcal", tokenLabelDesc, hgcalDefaults);
0217   descriptions.add(("hltEgammaHLTPhase2ExtraProducer"), desc);
0218 }
0219
0220 void EgammaHLTPhase2ExtraProducer::produce(edm::StreamID streamID,
0221                                            edm::Event& event,
0222                                            const edm::EventSetup& eventSetup) const {
0223   auto egTrigObjs = event.getHandle(tokens_.egTrigObjs);
0224
0225   auto trkParts = event.getHandle(tokens_.trkParts);
0226   auto l1trks = event.getHandle(tokens_.l1Trks);
0227   auto l1TrkToTrkPartMap = event.getHandle(tokens_.l1TrkToTrkPartMap);
0228
0229   auto const caloGeomHandle = eventSetup.getHandle(caloGeomToken_);
0230
0231   for (const auto& tokenLabel : tokens_.hgcal) {
0232     auto handle = event.getHandle(tokenLabel.first);
0233     auto recHits = filterRecHits(*egTrigObjs, handle, *caloGeomHandle);
0234     event.put(std::move(recHits), tokenLabel.second);
0235   }
0236
0237   auto storeCountRecHits = [&event](const auto& tokenLabels, const auto& thresholds, const std::string& prefixLabel) {
0238     for (const auto& tokenLabel : tokenLabels) {
0239       auto handle = event.getHandle(tokenLabel.first);
0240       auto count = countRecHits(handle, thresholds);
0241       for (size_t thresNr = 0; thresNr < thresholds.size(); thresNr++) {
0242         const auto& thres = thresholds[thresNr];
0243         event.put(std::move(count[thresNr]),
0244                   prefixLabel + tokenLabel.second + "Thres" + convertToProdNameStr(thres) + "GeV");
0245       }
0246     }
0247   };
0248   storeCountRecHits(tokens_.hgcal, recHitCountThresholds_, "countHgcalRecHits");
0249
0250   auto hgcalLayerClusters = event.getHandle(tokens_.hgcalLayerClusters);
0251   auto hgcalLayerClustersTime = event.getHandle(tokens_.hgcalLayerClustersTime);
0252   std::vector<edm::Ref<reco::CaloClusterCollection>> orgHGCalLayerClusterRefs;
0253   auto hgcalLayerClustersFiltered = filterObjs(*egTrigObjs, hgcalLayerClusters, orgHGCalLayerClusterRefs);
0254   std::vector<std::pair<float, float>> timesFiltered;
0255   timesFiltered.reserve(orgHGCalLayerClusterRefs.size());
0256   for (auto& clusRef : orgHGCalLayerClusterRefs) {
0257     timesFiltered.push_back((*hgcalLayerClustersTime)[clusRef]);
0258   }
0259   auto hgcalLayerClustersFilteredHandle = event.put(std::move(hgcalLayerClustersFiltered), "hgcalLayerClusters");
0260   event.put(makeValueMap(hgcalLayerClustersFilteredHandle, timesFiltered), "hgcalLayerClustersTime");
0261
0262   std::vector<edm::Ref<L1TrackCollection>> orgL1TrkRefs;
0263   auto l1TrksFiltered = filterObjs(*egTrigObjs, l1trks, orgL1TrkRefs);
0264   std::vector<edm::Ref<TrackingParticleCollection>> orgTPRefs;
0265   auto trkPartsFiltered = filterObjs(*egTrigObjs, trkParts, orgTPRefs);
0266
0267   //meh should make this edm::Ref<T>::key_type
0268   std::unordered_map<size_t, size_t> orgTPIndxToNewIndx;
0269   for (size_t refNr = 0; refNr < orgTPRefs.size(); refNr++) {
0270     orgTPIndxToNewIndx.insert(std::make_pair(orgTPRefs[refNr].key(), refNr));
0271   }
0272
0273   edm::OrphanHandle<L1TrackCollection> l1TrksFiltHandle = event.put(std::move(l1TrksFiltered));
0274   edm::OrphanHandle<TrackingParticleCollection> trkPartsFiltHandle = event.put(std::move(trkPartsFiltered));
0275
0276   auto l1TrkExtraColl = std::make_unique<L1TrackTruthPairCollection>();
0277
0278   for (size_t l1TrkNr = 0; l1TrkNr < orgL1TrkRefs.size(); l1TrkNr++) {
0279     auto orgTrkRef = orgL1TrkRefs[l1TrkNr];
0280     auto orgTrkPtr = edm::refToPtr(orgTrkRef);
0281     int flags = 0;
0282     if (l1TrkToTrkPartMap->isGenuine(orgTrkPtr))
0283       flags |= L1TrackTruthPair::StatusFlags::IsGenuine;
0284     if (l1TrkToTrkPartMap->isLooselyGenuine(orgTrkPtr))
0285       flags |= L1TrackTruthPair::StatusFlags::IsLooselyGenuine;
0286     if (l1TrkToTrkPartMap->isCombinatoric(orgTrkPtr))
0287       flags |= L1TrackTruthPair::StatusFlags::IsCombinatoric;
0288     if (l1TrkToTrkPartMap->isUnknown(orgTrkPtr))
0289       flags |= L1TrackTruthPair::StatusFlags::IsUnknown;
0290
0291     auto orgTPRef = l1TrkToTrkPartMap->findTrackingParticlePtr(orgTrkPtr);
0292     auto getNewTPRef = [&orgTPIndxToNewIndx, &orgTPRef, &trkPartsFiltHandle]() {
0293       auto newIndexPair = orgTPIndxToNewIndx.find(orgTPRef.key());
0294       if (newIndexPair != orgTPIndxToNewIndx.end()) {
0295         return edm::Ref<TrackingParticleCollection>(trkPartsFiltHandle, newIndexPair->second);
0296       } else
0297         return edm::Ref<TrackingParticleCollection>(trkPartsFiltHandle.id());
0298     };
0299     auto newTPRef = getNewTPRef();
0300     edm::Ref<L1TrackCollection> newL1TrkRef(l1TrksFiltHandle, l1TrkNr);
0301
0302     L1TrackTruthPair l1TrkExtra(newL1TrkRef, newTPRef, flags);
0303     l1TrkExtraColl->push_back(l1TrkExtra);
0304   }
0305   event.put(std::move(l1TrkExtraColl));
0306 }
0307
0308 template <typename CollType, typename RefType>
0309 std::unique_ptr<CollType> EgammaHLTPhase2ExtraProducer::filterObjs(const trigger::EgammaObjectCollection& egTrigObjs,
0310                                                                    const edm::Handle<CollType>& objs,
0311                                                                    std::vector<RefType>& orgRefs,
0312                                                                    float maxDR2) const {
0313   auto filteredObjs = std::make_unique<CollType>();
0314   orgRefs.clear();
0315   if (!objs.isValid())
0316     return filteredObjs;
0317
0318   //so because each egamma object can have multiple eta/phi pairs
0319   //easier to just make a temp vector and then copy that in with the +pi and  +pi/2
0320   std::vector<std::pair<float, float>> etaPhisTmp;
0321   for (const auto& egTrigObj : egTrigObjs) {
0322     if (egTrigObj.pt() >= minPtToSaveHits_) {
0323       etaPhisTmp.push_back({egTrigObj.eta(), egTrigObj.phi()});
0324       //also save the eta /phi of all gsf tracks with the object
0325       for (const auto& gsfTrk : egTrigObj.gsfTracks()) {
0326         etaPhisTmp.push_back({gsfTrk->eta(), gsfTrk->phi()});
0327       }
0328     }
0329   }
0330   std::vector<std::pair<float, float>> etaPhis;
0331   for (const auto& etaPhi : etaPhisTmp) {
0332     etaPhis.push_back(etaPhi);
0333     if (saveHitsPlusPi_)
0334       etaPhis.push_back({etaPhi.first, etaPhi.second + 3.14159});
0335     if (saveHitsPlusHalfPi_)
0336       etaPhis.push_back({etaPhi.first, etaPhi.second + 3.14159 / 2.});
0337   }
0338
0339   auto deltaR2Match = [&etaPhis, &maxDR2](float eta, float phi) {
0340     for (auto& etaPhi : etaPhis) {
0341       if (reco::deltaR2(eta, phi, etaPhi.first, etaPhi.second) < maxDR2)
0342         return true;
0343     }
0344     return false;
0345   };
0346
0347   for (size_t objNr = 0; objNr < objs->size(); objNr++) {
0348     RefType ref(objs, objNr);
0349     if (deltaR2Match(ref->eta(), ref->phi())) {
0350       filteredObjs->push_back(*ref);
0351       orgRefs.push_back(ref);
0352     }
0353   }
0354   return filteredObjs;
0355 }
0356
0357 template <typename RecHitCollection>
0358 std::unique_ptr<RecHitCollection> EgammaHLTPhase2ExtraProducer::filterRecHits(
0359     const trigger::EgammaObjectCollection& egTrigObjs,
0360     const edm::Handle<RecHitCollection>& recHits,
0361     const CaloGeometry& geom,
0362     float maxDR2) const {
0363   auto filteredHits = std::make_unique<RecHitCollection>();
0364   if (!recHits.isValid())
0365     return filteredHits;
0366
0367   std::vector<std::pair<float, float>> etaPhis;
0368   for (const auto& egTrigObj : egTrigObjs) {
0369     if (egTrigObj.pt() >= minPtToSaveHits_) {
0370       etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi()});
0371       if (saveHitsPlusPi_)
0372         etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi() + 3.14159});
0373       if (saveHitsPlusHalfPi_)
0374         etaPhis.push_back({egTrigObj.eta(), egTrigObj.phi() + 3.14159 / 2.});
0375     }
0376   }
0377   auto deltaR2Match = [&etaPhis, &maxDR2](const GlobalPoint& pos) {
0378     float eta = pos.eta();
0379     float phi = pos.phi();
0380     for (auto& etaPhi : etaPhis) {
0381       if (reco::deltaR2(eta, phi, etaPhi.first, etaPhi.second) < maxDR2)
0382         return true;
0383     }
0384     return false;
0385   };
0386
0387   for (auto& hit : *recHits) {
0388     const CaloSubdetectorGeometry* subDetGeom = geom.getSubdetectorGeometry(hit.id());
0389     if (subDetGeom) {
0390       auto cellGeom = subDetGeom->getGeometry(hit.id());
0391       if (deltaR2Match(cellGeom->getPosition()))
0392         filteredHits->push_back(hit);
0393     } else {
0394       throw cms::Exception("GeomError") << "could not get geometry for det id " << hit.id().rawId();
0395     }
0396   }
0397   return filteredHits;
0398 }
0399
0400 #include "FWCore/Framework/interface/MakerMacros.h"
0401 DEFINE_FWK_MODULE(EgammaHLTPhase2ExtraProducer);