EgammaClusterProducers/src/InterestingDetIdCollectionProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    InterestingDetIdCollectionProducer
0004 // Class:      InterestingDetIdCollectionProducer
0005 //
0006 /**\class InterestingDetIdCollectionProducer
0007 Original author: Paolo Meridiani PH/CMG
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.in a region around  the seed of the cluster collection specified
0014       by paramter basicClusters. 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/BasicCluster.h"
0029 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.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/Frameworkfwd.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/Framework/interface/global/EDProducer.h"
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 #include "FWCore/Utilities/interface/ESGetToken.h"
0039 #include "FWCore/Utilities/interface/InputTag.h"
0040 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0041 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0042 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0043 #include "RecoEcal/EgammaCoreTools/interface/EcalNextToDeadChannelRcd.h"
0044 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0045 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0046 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0047
0048 #include <memory>
0049
0050 class InterestingDetIdCollectionProducer : public edm::global::EDProducer<> {
0051 public:
0052   //! ctor
0053   explicit InterestingDetIdCollectionProducer(const edm::ParameterSet&);
0054   //! producer
0055   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0056
0057 private:
0058   // ----------member data ---------------------------
0059   edm::EDGetTokenT<EcalRecHitCollection> recHitsToken_;
0060   edm::EDGetTokenT<reco::BasicClusterCollection> basicClustersToken_;
0061   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0062   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> sevLVToken_;
0063   edm::ESGetToken<EcalNextToDeadChannel, EcalNextToDeadChannelRcd> nextToDeadToken_;
0064
0065   std::string interestingDetIdCollection_;
0066   int minimalEtaSize_;
0067   int minimalPhiSize_;
0068
0069   int severityLevel_;
0070   bool keepNextToDead_;
0071   bool keepNextToBoundary_;
0072 };
0073
0074 #include "FWCore/Framework/interface/MakerMacros.h"
0075 DEFINE_FWK_MODULE(InterestingDetIdCollectionProducer);
0076
0077 InterestingDetIdCollectionProducer::InterestingDetIdCollectionProducer(const edm::ParameterSet& iConfig) {
0078   recHitsToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsLabel"));
0079   basicClustersToken_ =
0080       consumes<reco::BasicClusterCollection>(iConfig.getParameter<edm::InputTag>("basicClustersLabel"));
0081   caloTopologyToken_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0082   sevLVToken_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0083   nextToDeadToken_ = esConsumes<EcalNextToDeadChannel, EcalNextToDeadChannelRcd>();
0084
0085   interestingDetIdCollection_ = iConfig.getParameter<std::string>("interestingDetIdCollection");
0086
0087   minimalEtaSize_ = iConfig.getParameter<int>("etaSize");
0088   minimalPhiSize_ = iConfig.getParameter<int>("phiSize");
0089   if (minimalPhiSize_ % 2 == 0 || minimalEtaSize_ % 2 == 0)
0090     edm::LogError("InterestingDetIdCollectionProducerError") << "Size of eta/phi should be odd numbers";
0091
0092   //register your products
0093   produces<DetIdCollection>(interestingDetIdCollection_);
0094
0095   severityLevel_ = iConfig.getParameter<int>("severityLevel");
0096   keepNextToDead_ = iConfig.getParameter<bool>("keepNextToDead");
0097   keepNextToBoundary_ = iConfig.getParameter<bool>("keepNextToBoundary");
0098 }
0099
0100 // ------------ method called to produce the data  ------------
0101 void InterestingDetIdCollectionProducer::produce(edm::StreamID,
0102                                                  edm::Event& iEvent,
0103                                                  const edm::EventSetup& iSetup) const {
0104   using namespace edm;
0105   using namespace std;
0106
0107   auto const& caloTopology = iSetup.getData(caloTopologyToken_);
0108
0109   auto const& severity = iSetup.getData(sevLVToken_);
0110
0111   // take BasicClusters
0112   Handle<reco::BasicClusterCollection> pClusters;
0113   iEvent.getByToken(basicClustersToken_, pClusters);
0114
0115   // take EcalRecHits
0116   Handle<EcalRecHitCollection> recHitsHandle;
0117   iEvent.getByToken(recHitsToken_, recHitsHandle);
0118
0119   //Create empty output collections
0120   std::vector<DetId> indexToStore;
0121   indexToStore.reserve(1000);
0122
0123   reco::BasicClusterCollection::const_iterator clusIt;
0124
0125   std::vector<DetId> xtalsToStore;
0126   xtalsToStore.reserve(50);
0127   for (clusIt = pClusters->begin(); clusIt != pClusters->end(); clusIt++) {
0128     const std::vector<std::pair<DetId, float> >& clusterDetIds = clusIt->hitsAndFractions();
0129     for (const auto& detidpair : clusterDetIds) {
0130       indexToStore.push_back(detidpair.first);
0131     }
0132
0133     //below checks and additional code are only relevant for EB/EE, not for ES
0134     if (clusterDetIds.front().first.subdetId() == EcalBarrel || clusterDetIds.front().first.subdetId() == EcalEndcap) {
0135       std::vector<std::pair<DetId, float> >::const_iterator posCurrent;
0136
0137       float eMax = 0.;
0138       DetId eMaxId(0);
0139
0140       EcalRecHit testEcalRecHit;
0141
0142       for (posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); posCurrent++) {
0143         EcalRecHitCollection::const_iterator itt = recHitsHandle->find((*posCurrent).first);
0144         if ((!((*posCurrent).first.null())) && (itt != recHitsHandle->end()) && ((*itt).energy() > eMax)) {
0145           eMax = (*itt).energy();
0146           eMaxId = (*itt).id();
0147         }
0148       }
0149
0150       if (eMaxId.null())
0151         continue;
0152
0153       const CaloSubdetectorTopology* topology = caloTopology.getSubdetectorTopology(eMaxId.det(), eMaxId.subdetId());
0154
0155       xtalsToStore = topology->getWindow(eMaxId, minimalEtaSize_, minimalPhiSize_);
0156
0157       for (const auto& detid : xtalsToStore) {
0158         indexToStore.push_back(detid);
0159       }
0160     }
0161   }
0162
0163   if (severityLevel_ >= 0 || keepNextToDead_ || keepNextToBoundary_) {
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
0195   //unify the vector
0196   std::sort(indexToStore.begin(), indexToStore.end());
0197   std::unique(indexToStore.begin(), indexToStore.end());
0198
0199   iEvent.put(std::make_unique<DetIdCollection>(indexToStore), interestingDetIdCollection_);
0200 }