EgammaClusterProducers/src/ReducedESRecHitCollectionProducer.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "DataFormats/DetId/interface/DetId.h"
0003 #include "DataFormats/DetId/interface/DetIdCollection.h"
0004 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0005 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0006 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0007 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0008 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/stream/EDProducer.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/ESGetToken.h"
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 #include "FWCore/Utilities/interface/transform.h"
0019 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0020 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0021 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0022 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0023 #include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
0024
0025 #include <iostream>
0026 #include <map>
0027 #include <set>
0028 #include <string>
0029 #include <vector>
0030
0031 class ReducedESRecHitCollectionProducer : public edm::stream::EDProducer<> {
0032 public:
0033   ReducedESRecHitCollectionProducer(const edm::ParameterSet& pset);
0034   ~ReducedESRecHitCollectionProducer() override;
0035   void beginRun(edm::Run const&, const edm::EventSetup&) final;
0036   void produce(edm::Event& e, const edm::EventSetup& c) override;
0037   void collectIds(const ESDetId strip1, const ESDetId strip2, const int& row = 0);
0038
0039 private:
0040   const EcalPreshowerGeometry* geometry_p;
0041   std::unique_ptr<CaloSubdetectorTopology> topology_p;
0042
0043   double scEtThresh_;
0044
0045   edm::EDGetTokenT<ESRecHitCollection> InputRecHitES_;
0046   edm::EDGetTokenT<reco::SuperClusterCollection> InputSuperClusterEE_;
0047   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0048   std::string OutputLabelES_;
0049   std::vector<edm::EDGetTokenT<DetIdCollection>> interestingDetIdCollections_;
0050   std::vector<edm::EDGetTokenT<DetIdCollection>>
0051       interestingDetIdCollectionsNotToClean_;  //theres a hard coded cut on rec-hit quality which some collections would prefer not to have...
0052
0053   std::set<DetId> collectedIds_;
0054 };
0055
0056 #include "FWCore/Framework/interface/MakerMacros.h"
0057 DEFINE_FWK_MODULE(ReducedESRecHitCollectionProducer);
0058
0059 using namespace edm;
0060 using namespace std;
0061 using namespace reco;
0062
0063 ReducedESRecHitCollectionProducer::ReducedESRecHitCollectionProducer(const edm::ParameterSet& ps)
0064     : geometry_p(nullptr) {
0065   scEtThresh_ = ps.getParameter<double>("scEtThreshold");
0066
0067   InputRecHitES_ = consumes<ESRecHitCollection>(ps.getParameter<edm::InputTag>("EcalRecHitCollectionES"));
0068   InputSuperClusterEE_ =
0069       consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("EndcapSuperClusterCollection"));
0070   caloGeometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>();
0071
0072   OutputLabelES_ = ps.getParameter<std::string>("OutputLabel_ES");
0073
0074   interestingDetIdCollections_ =
0075       edm::vector_transform(ps.getParameter<std::vector<edm::InputTag>>("interestingDetIds"),
0076                             [this](edm::InputTag const& tag) { return consumes<DetIdCollection>(tag); });
0077
0078   interestingDetIdCollectionsNotToClean_ =
0079       edm::vector_transform(ps.getParameter<std::vector<edm::InputTag>>("interestingDetIdsNotToClean"),
0080                             [this](edm::InputTag const& tag) { return consumes<DetIdCollection>(tag); });
0081
0082   produces<EcalRecHitCollection>(OutputLabelES_);
0083 }
0084
0085 ReducedESRecHitCollectionProducer::~ReducedESRecHitCollectionProducer() = default;
0086
0087 void ReducedESRecHitCollectionProducer::beginRun(edm::Run const&, const edm::EventSetup& iSetup) {
0088   ESHandle<CaloGeometry> geoHandle = iSetup.getHandle(caloGeometryToken_);
0089   geometry_p =
0090       dynamic_cast<const EcalPreshowerGeometry*>(geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower));
0091   if (!geometry_p) {
0092     edm::LogError("WrongGeometry") << "could not cast the subdet geometry to preshower geometry";
0093   }
0094
0095   if (geometry_p)
0096     topology_p = std::make_unique<EcalPreshowerTopology>();
0097 }
0098
0099 void ReducedESRecHitCollectionProducer::produce(edm::Event& e, const edm::EventSetup& iSetup) {
0100   edm::Handle<ESRecHitCollection> ESRecHits_;
0101   e.getByToken(InputRecHitES_, ESRecHits_);
0102
0103   auto output = std::make_unique<EcalRecHitCollection>();
0104
0105   edm::Handle<reco::SuperClusterCollection> pEndcapSuperClusters;
0106   e.getByToken(InputSuperClusterEE_, pEndcapSuperClusters);
0107   {
0108     const reco::SuperClusterCollection* eeSuperClusters = pEndcapSuperClusters.product();
0109
0110     for (reco::SuperClusterCollection::const_iterator isc = eeSuperClusters->begin(); isc != eeSuperClusters->end();
0111          ++isc) {
0112       if (isc->energy() < scEtThresh_)
0113         continue;
0114       if (fabs(isc->eta()) < 1.65 || fabs(isc->eta()) > 2.6)
0115         continue;
0116       //cout<<"SC energy : "<<isc->energy()<<" "<<isc->eta()<<endl;
0117
0118       //Int_t nBC = 0;
0119       reco::CaloCluster_iterator ibc = isc->clustersBegin();
0120       for (; ibc != isc->clustersEnd(); ++ibc) {
0121         //cout<<"BC : "<<nBC<<endl;
0122
0123         const GlobalPoint point((*ibc)->x(), (*ibc)->y(), (*ibc)->z());
0124
0125         ESDetId esId1 = geometry_p->getClosestCellInPlane(point, 1);
0126         ESDetId esId2 = geometry_p->getClosestCellInPlane(point, 2);
0127
0128         collectIds(esId1, esId2, 0);
0129         collectIds(esId1, esId2, 1);
0130         collectIds(esId1, esId2, -1);
0131
0132         //nBC++;
0133       }
0134     }
0135   }
0136
0137   edm::Handle<DetIdCollection> detId;
0138   for (unsigned int t = 0; t < interestingDetIdCollections_.size(); ++t) {
0139     e.getByToken(interestingDetIdCollections_[t], detId);
0140     if (!detId.isValid()) {
0141       Labels labels;
0142       labelsForToken(interestingDetIdCollections_[t], labels);
0143       edm::LogError("MissingInput") << "no reason to skip detid from : (" << labels.module << ", "
0144                                     << labels.productInstance << ", " << labels.process << ")" << std::endl;
0145       continue;
0146     }
0147     collectedIds_.insert(detId->begin(), detId->end());
0148   }
0149
0150   //screw it, cant think of a better solution, not the best but lets run over all the rec hits, remove the ones failing cleaning
0151   //and then merge in the collection not to be cleaned
0152   //mainly as I suspect its more efficient to find an object in the DetIdSet rather than the rec-hit in the rec-hit collecition
0153   //with only a det id
0154   //if its a CPU issues then revisit
0155   for (const auto& hit : *ESRecHits_) {
0156     if (hit.recoFlag() == 1 || hit.recoFlag() == 14 ||
0157         (hit.recoFlag() <= 10 && hit.recoFlag() >= 5)) {  //right we might need to erase it from the collection
0158       auto idIt = collectedIds_.find(hit.id());
0159       if (idIt != collectedIds_.end())
0160         collectedIds_.erase(idIt);
0161     }
0162   }
0163
0164   for (const auto& token : interestingDetIdCollectionsNotToClean_) {
0165     e.getByToken(token, detId);
0166     if (!detId.isValid()) {  //meh might as well keep the warning
0167       Labels labels;
0168       labelsForToken(token, labels);
0169       edm::LogError("MissingInput") << "no reason to skip detid from : (" << labels.module << ", "
0170                                     << labels.productInstance << ", " << labels.process << ")" << std::endl;
0171       continue;
0172     }
0173     collectedIds_.insert(detId->begin(), detId->end());
0174   }
0175
0176   output->reserve(collectedIds_.size());
0177   EcalRecHitCollection::const_iterator it;
0178   for (it = ESRecHits_->begin(); it != ESRecHits_->end(); ++it) {
0179     if (collectedIds_.find(it->id()) != collectedIds_.end()) {
0180       output->push_back(*it);
0181     }
0182   }
0183   collectedIds_.clear();
0184
0185   e.put(std::move(output), OutputLabelES_);
0186 }
0187
0188 void ReducedESRecHitCollectionProducer::collectIds(const ESDetId esDetId1, const ESDetId esDetId2, const int& row) {
0189   //cout<<row<<endl;
0190
0191   map<DetId, const EcalRecHit*>::iterator it;
0192   map<DetId, int>::iterator itu;
0193   ESDetId next;
0194   ESDetId strip1;
0195   ESDetId strip2;
0196
0197   strip1 = esDetId1;
0198   strip2 = esDetId2;
0199
0200   EcalPreshowerNavigator theESNav1(strip1, topology_p.get());
0201   theESNav1.setHome(strip1);
0202
0203   EcalPreshowerNavigator theESNav2(strip2, topology_p.get());
0204   theESNav2.setHome(strip2);
0205
0206   if (row == 1) {
0207     if (strip1 != ESDetId(0))
0208       strip1 = theESNav1.north();
0209     if (strip2 != ESDetId(0))
0210       strip2 = theESNav2.east();
0211   } else if (row == -1) {
0212     if (strip1 != ESDetId(0))
0213       strip1 = theESNav1.south();
0214     if (strip2 != ESDetId(0))
0215       strip2 = theESNav2.west();
0216   }
0217
0218   // Plane 1
0219   if (strip1 == ESDetId(0)) {
0220   } else {
0221     collectedIds_.insert(strip1);
0222     //cout<<"center : "<<strip1<<endl;
0223     // east road
0224     for (int i = 0; i < 15; ++i) {
0225       next = theESNav1.east();
0226       //cout<<"east : "<<i<<" "<<next<<endl;
0227       if (next != ESDetId(0)) {
0228         collectedIds_.insert(next);
0229       } else {
0230         break;
0231       }
0232     }
0233
0234     // west road
0235     theESNav1.setHome(strip1);
0236     theESNav1.home();
0237     for (int i = 0; i < 15; ++i) {
0238       next = theESNav1.west();
0239       //cout<<"west : "<<i<<" "<<next<<endl;
0240       if (next != ESDetId(0)) {
0241         collectedIds_.insert(next);
0242       } else {
0243         break;
0244       }
0245     }
0246   }
0247
0248   if (strip2 == ESDetId(0)) {
0249   } else {
0250     collectedIds_.insert(strip2);
0251     //cout<<"center : "<<strip2<<endl;
0252     // north road
0253     for (int i = 0; i < 15; ++i) {
0254       next = theESNav2.north();
0255       //cout<<"north : "<<i<<" "<<next<<endl;
0256       if (next != ESDetId(0)) {
0257         collectedIds_.insert(next);
0258       } else {
0259         break;
0260       }
0261     }
0262
0263     // south road
0264     theESNav2.setHome(strip2);
0265     theESNav2.home();
0266     for (int i = 0; i < 15; ++i) {
0267       next = theESNav2.south();
0268       //cout<<"south : "<<i<<" "<<next<<endl;
0269       if (next != ESDetId(0)) {
0270         collectedIds_.insert(next);
0271       } else {
0272         break;
0273       }
0274     }
0275   }
0276 }