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