EgammaClusterProducers/src/InterestingDetIdFromSuperClusterProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    InterestingDetIdFromSuperClusterProducer
0004 // Class:      InterestingDetIdFromSuperClusterProducer
0005 //
0006 /**\class InterestingDetIdFromSuperClusterProducer
0007 Adapted from InterestingDetIdCollectionProducer by J.Bendavid
0008
0009 Make a collection of detids to be kept tipically in a AOD rechit collection
0010
0011 The following classes of "interesting id" are considered
0012
0013     1.All rechits included in all subclusters, plus in a region around  the seed of each subcluster
0014       The size of the region is specified by
0015       minimalEtaSize_, minimalPhiSize_
0016
0017     2. if the severity of the hit is >= severityLevel_
0018        If severityLevel=0 this class is ignored
0019
0020     3. Channels next to dead ones,  keepNextToDead_ is true
0021     4. Channels next to the EB/EE transition if keepNextToBoundary_ is true
0022 */
0023
0024 #include "DataFormats/DetId/interface/DetIdCollection.h"
0025 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0026 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0027 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0028 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0029 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/Framework/interface/global/EDProducer.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "FWCore/Utilities/interface/ESGetToken.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0040 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0041 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0042 #include "RecoEcal/EgammaCoreTools/interface/EcalNextToDeadChannelRcd.h"
0043 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0044 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0045 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0046
0047 #include <memory>
0048
0049 class InterestingDetIdFromSuperClusterProducer : public edm::global::EDProducer<> {
0050 public:
0051   //! ctor
0052   explicit InterestingDetIdFromSuperClusterProducer(const edm::ParameterSet&);
0053   //! producer
0054   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0055
0056 private:
0057   // ----------member data ---------------------------
0058   edm::EDGetTokenT<EcalRecHitCollection> recHitsToken_;
0059   edm::EDGetTokenT<reco::SuperClusterCollection> superClustersToken_;
0060   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0061   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> severityLevelToken_;
0062   edm::ESGetToken<EcalNextToDeadChannel, EcalNextToDeadChannelRcd> nextToDeadToken_;
0063   std::string interestingDetIdCollection_;
0064   int minimalEtaSize_;
0065   int minimalPhiSize_;
0066
0067   int severityLevel_;
0068   bool keepNextToDead_;
0069   bool keepNextToBoundary_;
0070 };
0071
0072 #include "FWCore/Framework/interface/MakerMacros.h"
0073 DEFINE_FWK_MODULE(InterestingDetIdFromSuperClusterProducer);
0074
0075 InterestingDetIdFromSuperClusterProducer::InterestingDetIdFromSuperClusterProducer(const edm::ParameterSet& iConfig) {
0076   recHitsToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsLabel"));
0077   superClustersToken_ =
0078       consumes<reco::SuperClusterCollection>(iConfig.getParameter<edm::InputTag>("superClustersLabel"));
0079   caloTopologyToken_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0080   severityLevelToken_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0081   interestingDetIdCollection_ = iConfig.getParameter<std::string>("interestingDetIdCollection");
0082
0083   minimalEtaSize_ = iConfig.getParameter<int>("etaSize");
0084   minimalPhiSize_ = iConfig.getParameter<int>("phiSize");
0085   if (minimalPhiSize_ % 2 == 0 || minimalEtaSize_ % 2 == 0)
0086     edm::LogError("InterestingDetIdFromSuperClusterProducerError") << "Size of eta/phi should be odd numbers";
0087
0088   //register your products
0089   produces<DetIdCollection>(interestingDetIdCollection_);
0090
0091   severityLevel_ = iConfig.getParameter<int>("severityLevel");
0092   keepNextToDead_ = iConfig.getParameter<bool>("keepNextToDead");
0093   keepNextToBoundary_ = iConfig.getParameter<bool>("keepNextToBoundary");
0094   if (keepNextToDead_) {
0095     nextToDeadToken_ = esConsumes<EcalNextToDeadChannel, EcalNextToDeadChannelRcd>();
0096   }
0097 }
0098
0099 // ------------ method called to produce the data  ------------
0100 void InterestingDetIdFromSuperClusterProducer::produce(edm::StreamID,
0101                                                        edm::Event& iEvent,
0102                                                        const edm::EventSetup& iSetup) const {
0103   using namespace edm;
0104   using namespace std;
0105
0106   auto const& caloTopology = iSetup.getData(caloTopologyToken_);
0107
0108   auto const& severity = iSetup.getData(severityLevelToken_);
0109   // take BasicClusters
0110   Handle<reco::SuperClusterCollection> pClusters;
0111   iEvent.getByToken(superClustersToken_, pClusters);
0112
0113   // take EcalRecHits
0114   Handle<EcalRecHitCollection> recHitsHandle;
0115   iEvent.getByToken(recHitsToken_, recHitsHandle);
0116
0117   //Create empty output collections
0118   std::vector<DetId> indexToStore;
0119   indexToStore.reserve(1000);
0120
0121   reco::SuperClusterCollection::const_iterator sclusIt;
0122
0123   std::vector<DetId> xtalsToStore;
0124   xtalsToStore.reserve(50);
0125
0126   //loop over superclusters
0127   for (sclusIt = pClusters->begin(); sclusIt != pClusters->end(); sclusIt++) {
0128     //loop over subclusters
0129     for (reco::CaloCluster_iterator clusIt = sclusIt->clustersBegin(); clusIt != sclusIt->clustersEnd(); ++clusIt) {
0130       //PG barrel
0131
0132       float eMax = 0.;
0133       DetId eMaxId(0);
0134
0135       std::vector<std::pair<DetId, float> > clusterDetIds = (*clusIt)->hitsAndFractions();
0136       std::vector<std::pair<DetId, float> >::iterator posCurrent;
0137
0138       EcalRecHit testEcalRecHit;
0139
0140       for (posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); posCurrent++) {
0141         EcalRecHitCollection::const_iterator itt = recHitsHandle->find((*posCurrent).first);
0142         if ((!((*posCurrent).first.null())) && (itt != recHitsHandle->end()) && ((*itt).energy() > eMax)) {
0143           eMax = (*itt).energy();
0144           eMaxId = (*itt).id();
0145         }
0146       }
0147
0148       if (eMaxId.null())
0149         continue;
0150
0151       const CaloSubdetectorTopology* topology = caloTopology.getSubdetectorTopology(eMaxId.det(), eMaxId.subdetId());
0152
0153       xtalsToStore = topology->getWindow(eMaxId, minimalEtaSize_, minimalPhiSize_);
0154       std::vector<std::pair<DetId, float> > xtalsInClus = (*clusIt)->hitsAndFractions();
0155
0156       for (unsigned int ii = 0; ii < xtalsInClus.size(); ii++) {
0157         xtalsToStore.push_back(xtalsInClus[ii].first);
0158       }
0159
0160       indexToStore.insert(indexToStore.end(), xtalsToStore.begin(), xtalsToStore.end());
0161     }
0162   }
0163
0164   for (EcalRecHitCollection::const_iterator it = recHitsHandle->begin(); it != recHitsHandle->end(); ++it) {
0165     // also add recHits of dead TT if the corresponding TP is saturated
0166     if (it->checkFlag(EcalRecHit::kTPSaturated)) {
0167       indexToStore.push_back(it->id());
0168     }
0169     // add hits for severities above a threshold
0170     if (severityLevel_ >= 0 && severity.severityLevel(*it) >= severityLevel_) {
0171       indexToStore.push_back(it->id());
0172     }
0173     if (keepNextToDead_) {
0174       edm::ESHandle<EcalNextToDeadChannel> dch = iSetup.getHandle(nextToDeadToken_);
0175       // also keep channels next to dead ones
0176       if (EcalTools::isNextToDead(it->id(), *dch)) {
0177         indexToStore.push_back(it->id());
0178       }
0179     }
0180
0181     if (keepNextToBoundary_) {
0182       // keep channels around EB/EE boundary
0183       if (it->id().subdetId() == EcalBarrel) {
0184         EBDetId ebid(it->id());
0185         if (abs(ebid.ieta()) == 85)
0186           indexToStore.push_back(it->id());
0187       } else {
0188         if (EEDetId::isNextToRingBoundary(it->id()))
0189           indexToStore.push_back(it->id());
0190       }
0191     }
0192   }
0193
0194   //unify the vector
0195   std::sort(indexToStore.begin(), indexToStore.end());
0196   auto last = std::unique(indexToStore.begin(), indexToStore.end());
0197   indexToStore.erase(last, indexToStore.end());
0198
0199   auto detIdCollection = std::make_unique<DetIdCollection>(indexToStore);
0200
0201   iEvent.put(std::move(detIdCollection), interestingDetIdCollection_);
0202 }