File indexing completed on 2021-06-08 05:47:13
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
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
0052 explicit InterestingDetIdFromSuperClusterProducer(const edm::ParameterSet&);
0053 void beginRun(edm::Run const&, const edm::EventSetup&) final;
0054
0055 void produce(edm::Event&, const edm::EventSetup&) override;
0056
0057 private:
0058
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
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
0111 void InterestingDetIdFromSuperClusterProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0112 using namespace edm;
0113 using namespace std;
0114
0115
0116 Handle<reco::SuperClusterCollection> pClusters;
0117 iEvent.getByToken(superClustersToken_, pClusters);
0118
0119
0120 Handle<EcalRecHitCollection> recHitsHandle;
0121 iEvent.getByToken(recHitsToken_, recHitsHandle);
0122
0123
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
0133 for (sclusIt = pClusters->begin(); sclusIt != pClusters->end(); sclusIt++) {
0134
0135 for (reco::CaloCluster_iterator clusIt = sclusIt->clustersBegin(); clusIt != sclusIt->clustersEnd(); ++clusIt) {
0136
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
0172 if (it->checkFlag(EcalRecHit::kTPSaturated)) {
0173 indexToStore.push_back(it->id());
0174 }
0175
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
0182 if (EcalTools::isNextToDead(it->id(), *dch)) {
0183 indexToStore.push_back(it->id());
0184 }
0185 }
0186
0187 if (keepNextToBoundary_) {
0188
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
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 }