EgammaIsolationAlgos/plugins/EgammaIsoESDetIdCollectionProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EgammaIsoESDetIdCollectionProducer
0004 // Class:      EgammaIsoESDetIdCollectionProducer
0005 //
0006 /**\class EgammaIsoESDetIdCollectionProducer
0007
0008 author: Sam Harper (inspired by InterestingDetIdProducer)
0009
0010 Make a collection of detids to be kept in a AOD rechit collection
0011 These are all the ES DetIds of ES PFClusters associated to all PF clusters within dR of ele/pho/sc
0012 The aim is to save enough preshower info in the AOD to remake the PF clusters near an ele/pho/sc
0013 */
0014
0015 #include "DataFormats/DetId/interface/DetIdCollection.h"
0016 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0017 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0018 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0019 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0020 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0021 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0022 #include "DataFormats/Math/interface/deltaR.h"
0023 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0024 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/Framework/interface/global/EDProducer.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030
0031 class EgammaIsoESDetIdCollectionProducer : public edm::global::EDProducer<> {
0032 public:
0033   //! ctor
0034   explicit EgammaIsoESDetIdCollectionProducer(const edm::ParameterSet&);
0035   //! producer
0036   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0037
0038 private:
0039   void addDetIds(const reco::SuperCluster& superClus,
0040                  reco::PFClusterCollection clusters,
0041                  const reco::PFCluster::EEtoPSAssociation& eeClusToESMap,
0042                  std::vector<DetId>& detIdsToStore) const;
0043
0044   // ----------member data ---------------------------
0045   edm::EDGetTokenT<reco::PFCluster::EEtoPSAssociation> eeClusToESMapToken_;
0046   edm::EDGetTokenT<reco::PFClusterCollection> ecalPFClustersToken_;
0047   edm::EDGetTokenT<reco::SuperClusterCollection> superClustersToken_;
0048   edm::EDGetTokenT<reco::GsfElectronCollection> elesToken_;
0049   edm::EDGetTokenT<reco::PhotonCollection> phosToken_;
0050
0051   std::string interestingDetIdCollection_;
0052
0053   float minSCEt_;
0054   float minEleEt_;
0055   float minPhoEt_;
0056
0057   float maxDR_;
0058 };
0059
0060 #include "FWCore/Framework/interface/MakerMacros.h"
0061 DEFINE_FWK_MODULE(EgammaIsoESDetIdCollectionProducer);
0062
0063 EgammaIsoESDetIdCollectionProducer::EgammaIsoESDetIdCollectionProducer(const edm::ParameterSet& iConfig)
0064     : eeClusToESMapToken_{consumes(iConfig.getParameter<edm::InputTag>("eeClusToESMapLabel"))},
0065       ecalPFClustersToken_{consumes(iConfig.getParameter<edm::InputTag>("ecalPFClustersLabel"))},
0066       superClustersToken_{consumes(iConfig.getParameter<edm::InputTag>("superClustersLabel"))},
0067       elesToken_{consumes(iConfig.getParameter<edm::InputTag>("elesLabel"))},
0068       phosToken_{consumes(iConfig.getParameter<edm::InputTag>("phosLabel"))} {
0069   minSCEt_ = iConfig.getParameter<double>("minSCEt");
0070   minEleEt_ = iConfig.getParameter<double>("minEleEt");
0071   minPhoEt_ = iConfig.getParameter<double>("minPhoEt");
0072
0073   interestingDetIdCollection_ = iConfig.getParameter<std::string>("interestingDetIdCollection");
0074
0075   maxDR_ = iConfig.getParameter<double>("maxDR");
0076
0077   //register your products
0078   produces<DetIdCollection>(interestingDetIdCollection_);
0079 }
0080
0081 // ------------ method called to produce the data  ------------
0082 void EgammaIsoESDetIdCollectionProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0083   auto superClusters = iEvent.getHandle(superClustersToken_);
0084   auto eles = iEvent.getHandle(elesToken_);
0085   auto phos = iEvent.getHandle(phosToken_);
0086   auto eeClusToESMap = iEvent.getHandle(eeClusToESMapToken_);
0087   auto ecalPFClusters = iEvent.getHandle(ecalPFClustersToken_);
0088
0089   //Create empty output collections
0090   std::vector<DetId> indexToStore;
0091   indexToStore.reserve(100);
0092
0093   if (eles.isValid() && eeClusToESMap.isValid() && ecalPFClusters.isValid()) {
0094     for (auto& ele : *eles) {
0095       float scEt = ele.superCluster()->energy() * std::sin(ele.superCluster()->position().theta());
0096       if (scEt > minEleEt_ || ele.et() > minEleEt_)
0097         addDetIds(*ele.superCluster(), *ecalPFClusters, *eeClusToESMap, indexToStore);
0098     }
0099   }
0100   if (phos.isValid() && eeClusToESMap.isValid() && ecalPFClusters.isValid()) {
0101     for (auto& pho : *phos) {
0102       float scEt = pho.superCluster()->energy() * std::sin(pho.superCluster()->position().theta());
0103       if (scEt > minPhoEt_ || pho.et() > minPhoEt_)
0104         addDetIds(*pho.superCluster(), *ecalPFClusters, *eeClusToESMap, indexToStore);
0105     }
0106   }
0107   if (superClusters.isValid() && eeClusToESMap.isValid() && ecalPFClusters.isValid()) {
0108     for (auto& sc : *superClusters) {
0109       float scEt = sc.energy() * std::sin(sc.position().theta());
0110       if (scEt > minSCEt_)
0111         addDetIds(sc, *ecalPFClusters, *eeClusToESMap, indexToStore);
0112     }
0113   }
0114
0115   //unify the vector
0116   std::sort(indexToStore.begin(), indexToStore.end());
0117   std::unique(indexToStore.begin(), indexToStore.end());
0118
0119   auto detIdCollection = std::make_unique<DetIdCollection>(indexToStore);
0120
0121   iEvent.put(std::move(detIdCollection), interestingDetIdCollection_);
0122 }
0123
0124 //find all clusters within dR2<maxDR2_ of supercluster and then save det id's of hits of all es clusters matched to said ecal clusters
0125 void EgammaIsoESDetIdCollectionProducer::addDetIds(const reco::SuperCluster& superClus,
0126                                                    reco::PFClusterCollection clusters,
0127                                                    const reco::PFCluster::EEtoPSAssociation& eeClusToESMap,
0128                                                    std::vector<DetId>& detIdsToStore) const {
0129   const float scEta = superClus.eta();
0130   //  if(std::abs(scEta)+maxDR_<1.5) return; //not possible to have a endcap cluster, let alone one with preshower (eta>1.65) so exit without checking further
0131   const float scPhi = superClus.phi();
0132
0133   const float maxDR2 = maxDR_ * maxDR_;
0134
0135   for (size_t clusNr = 0; clusNr < clusters.size(); clusNr++) {
0136     const reco::PFCluster& clus = clusters[clusNr];
0137     if (clus.layer() == PFLayer::ECAL_ENDCAP && reco::deltaR2(scEta, scPhi, clus.eta(), clus.phi()) < maxDR2) {
0138       auto keyVal = std::make_pair(clusNr, edm::Ptr<reco::PFCluster>());
0139       const auto esClusters = std::equal_range(
0140           eeClusToESMap.begin(),
0141           eeClusToESMap.end(),
0142           keyVal,
0143           [](const reco::PFCluster::EEtoPSAssociation::value_type& rhs,  //roll on c++14, auto & lambda 4 evar!
0144              const reco::PFCluster::EEtoPSAssociation::value_type& lhs) -> bool { return rhs.first < lhs.first; });
0145       //   std::cout <<"cluster "<<clus.eta()<<"  had "<<std::distance(esClusters.first,esClusters.second)<<" es clusters"<<std::endl;
0146       for (auto esIt = esClusters.first; esIt != esClusters.second; ++esIt) {
0147         //   std::cout <<"es clus "<<esIt->second->hitsAndFractions().size()<<std::endl;
0148         for (const auto& hitAndFrac : esIt->second->hitsAndFractions()) {
0149           detIdsToStore.push_back(hitAndFrac.first);
0150         }
0151       }
0152
0153     }  //end of endcap & dR check
0154   }    //end of cluster loop
0155 }