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/stream/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::stream::EDProducer<> {
0051 public:
0052   //! ctor
0053   explicit InterestingDetIdCollectionProducer(const edm::ParameterSet&);
0054   void beginRun(edm::Run const&, const edm::EventSetup&) final;
0055   //! producer
0056   void produce(edm::Event&, const edm::EventSetup&) override;
0057
0058 private:
0059   // ----------member data ---------------------------
0060   edm::EDGetTokenT<EcalRecHitCollection> recHitsToken_;
0061   edm::EDGetTokenT<reco::BasicClusterCollection> basicClustersToken_;
0062   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0063   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> sevLVToken_;
0064   edm::ESGetToken<EcalNextToDeadChannel, EcalNextToDeadChannelRcd> nextToDeadToken_;
0065
0066   std::string interestingDetIdCollection_;
0067   int minimalEtaSize_;
0068   int minimalPhiSize_;
0069   const CaloTopology* caloTopology_;
0070
0071   int severityLevel_;
0072   const EcalSeverityLevelAlgo* severity_;
0073   bool keepNextToDead_;
0074   bool keepNextToBoundary_;
0075 };
0076
0077 #include "FWCore/Framework/interface/MakerMacros.h"
0078 DEFINE_FWK_MODULE(InterestingDetIdCollectionProducer);
0079
0080 InterestingDetIdCollectionProducer::InterestingDetIdCollectionProducer(const edm::ParameterSet& iConfig) {
0081   recHitsToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsLabel"));
0082   basicClustersToken_ =
0083       consumes<reco::BasicClusterCollection>(iConfig.getParameter<edm::InputTag>("basicClustersLabel"));
0084   caloTopologyToken_ = esConsumes<CaloTopology, CaloTopologyRecord, edm::Transition::BeginRun>();
0085   sevLVToken_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd, edm::Transition::BeginRun>();
0086   nextToDeadToken_ = esConsumes<EcalNextToDeadChannel, EcalNextToDeadChannelRcd>();
0087
0088   interestingDetIdCollection_ = iConfig.getParameter<std::string>("interestingDetIdCollection");
0089
0090   minimalEtaSize_ = iConfig.getParameter<int>("etaSize");
0091   minimalPhiSize_ = iConfig.getParameter<int>("phiSize");
0092   if (minimalPhiSize_ % 2 == 0 || minimalEtaSize_ % 2 == 0)
0093     edm::LogError("InterestingDetIdCollectionProducerError") << "Size of eta/phi should be odd numbers";
0094
0095   //register your products
0096   produces<DetIdCollection>(interestingDetIdCollection_);
0097
0098   severityLevel_ = iConfig.getParameter<int>("severityLevel");
0099   keepNextToDead_ = iConfig.getParameter<bool>("keepNextToDead");
0100   keepNextToBoundary_ = iConfig.getParameter<bool>("keepNextToBoundary");
0101 }
0102
0103 void InterestingDetIdCollectionProducer::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0104   edm::ESHandle<CaloTopology> theCaloTopology = iSetup.getHandle(caloTopologyToken_);
0105   caloTopology_ = &(*theCaloTopology);
0106
0107   edm::ESHandle<EcalSeverityLevelAlgo> sevLv = iSetup.getHandle(sevLVToken_);
0108   severity_ = sevLv.product();
0109 }
0110
0111 // ------------ method called to produce the data  ------------
0112 void InterestingDetIdCollectionProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0113   using namespace edm;
0114   using namespace std;
0115
0116   // take BasicClusters
0117   Handle<reco::BasicClusterCollection> pClusters;
0118   iEvent.getByToken(basicClustersToken_, pClusters);
0119
0120   // take EcalRecHits
0121   Handle<EcalRecHitCollection> recHitsHandle;
0122   iEvent.getByToken(recHitsToken_, recHitsHandle);
0123
0124   //Create empty output collections
0125   std::vector<DetId> indexToStore;
0126   indexToStore.reserve(1000);
0127
0128   reco::BasicClusterCollection::const_iterator clusIt;
0129
0130   std::vector<DetId> xtalsToStore;
0131   xtalsToStore.reserve(50);
0132   for (clusIt = pClusters->begin(); clusIt != pClusters->end(); clusIt++) {
0133     const std::vector<std::pair<DetId, float> >& clusterDetIds = clusIt->hitsAndFractions();
0134     for (const auto& detidpair : clusterDetIds) {
0135       indexToStore.push_back(detidpair.first);
0136     }
0137
0138     //below checks and additional code are only relevant for EB/EE, not for ES
0139     if (clusterDetIds.front().first.subdetId() == EcalBarrel || clusterDetIds.front().first.subdetId() == EcalEndcap) {
0140       std::vector<std::pair<DetId, float> >::const_iterator posCurrent;
0141
0142       float eMax = 0.;
0143       DetId eMaxId(0);
0144
0145       EcalRecHit testEcalRecHit;
0146
0147       for (posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); posCurrent++) {
0148         EcalRecHitCollection::const_iterator itt = recHitsHandle->find((*posCurrent).first);
0149         if ((!((*posCurrent).first.null())) && (itt != recHitsHandle->end()) && ((*itt).energy() > eMax)) {
0150           eMax = (*itt).energy();
0151           eMaxId = (*itt).id();
0152         }
0153       }
0154
0155       if (eMaxId.null())
0156         continue;
0157
0158       const CaloSubdetectorTopology* topology = caloTopology_->getSubdetectorTopology(eMaxId.det(), eMaxId.subdetId());
0159
0160       xtalsToStore = topology->getWindow(eMaxId, minimalEtaSize_, minimalPhiSize_);
0161
0162       for (const auto& detid : xtalsToStore) {
0163         indexToStore.push_back(detid);
0164       }
0165     }
0166   }
0167
0168   if (severityLevel_ >= 0 || keepNextToDead_ || keepNextToBoundary_) {
0169     for (EcalRecHitCollection::const_iterator it = recHitsHandle->begin(); it != recHitsHandle->end(); ++it) {
0170       // also add recHits of dead TT if the corresponding TP is saturated
0171       if (it->checkFlag(EcalRecHit::kTPSaturated)) {
0172         indexToStore.push_back(it->id());
0173       }
0174       // add hits for severities above a threshold
0175       if (severityLevel_ >= 0 && severity_->severityLevel(*it) >= severityLevel_) {
0176         indexToStore.push_back(it->id());
0177       }
0178       if (keepNextToDead_) {
0179         edm::ESHandle<EcalNextToDeadChannel> dch = iSetup.getHandle(nextToDeadToken_);
0180         // also keep channels next to dead ones
0181         if (EcalTools::isNextToDead(it->id(), *dch)) {
0182           indexToStore.push_back(it->id());
0183         }
0184       }
0185
0186       if (keepNextToBoundary_) {
0187         // keep channels around EB/EE boundary
0188         if (it->id().subdetId() == EcalBarrel) {
0189           EBDetId ebid(it->id());
0190           if (abs(ebid.ieta()) == 85)
0191             indexToStore.push_back(it->id());
0192         } else {
0193           if (EEDetId::isNextToRingBoundary(it->id()))
0194             indexToStore.push_back(it->id());
0195         }
0196       }
0197     }
0198   }
0199
0200   //unify the vector
0201   std::sort(indexToStore.begin(), indexToStore.end());
0202   std::unique(indexToStore.begin(), indexToStore.end());
0203
0204   iEvent.put(std::make_unique<DetIdCollection>(indexToStore), interestingDetIdCollection_);
0205 }