EgammaHLTProducers/plugins/EgammaHLTNxNClusterProducer.cc

0001 /**
0002
0003 Description: simple NxN ( 3x3 etc) clustering ,( for low energy photon reconstrution, currently used for pi0/eta HLT path)
0004
0005  Implementation:
0006      <Notes on implementation>
0007 */
0008
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/stream/EDProducer.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0017 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0018 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0019 #include "DataFormats/CaloRecHit/interface/CaloID.h"
0020
0021 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0022 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0023
0024 #include <vector>
0025 #include <memory>
0026 #include <ctime>
0027
0028 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0029 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0030 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0031 #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
0032 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0034
0035 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0036 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0037 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0038 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0039 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0040 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0041 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0042 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0043
0044 #include "TVector3.h"
0045
0046 class EgammaHLTNxNClusterProducer : public edm::stream::EDProducer<> {
0047 public:
0048   EgammaHLTNxNClusterProducer(const edm::ParameterSet &ps);
0049   ~EgammaHLTNxNClusterProducer() override;
0050
0051   void produce(edm::Event &, const edm::EventSetup &) override;
0052   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0053
0054 private:
0055   void makeNxNClusters(edm::Event &evt,
0056                        const edm::EventSetup &es,
0057                        const EcalRecHitCollection &hits,
0058                        const reco::CaloID::Detectors detector);
0059
0060   bool checkStatusOfEcalRecHit(const EcalChannelStatus &channelStatus, const EcalRecHit &rh);
0061
0062   //std::map<std::string,double> providedParameters;
0063
0064   const bool doBarrel_;
0065   const bool doEndcaps_;
0066   const edm::EDGetTokenT<EcalRecHitCollection> barrelHitProducer_;
0067   const edm::EDGetTokenT<EcalRecHitCollection> endcapHitProducer_;
0068   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0069   const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> ecalChannelStatusToken_;
0070
0071   const int clusEtaSize_;
0072   const int clusPhiSize_;
0073   const std::string barrelClusterCollection_;
0074   const std::string endcapClusterCollection_;
0075   const double clusSeedThr_;
0076   const double clusSeedThrEndCap_;
0077
0078   const bool useRecoFlag_;
0079   const int flagLevelRecHitsToUse_;
0080   const bool useDBStatus_;
0081   const int statusLevelRecHitsToUse_;
0082
0083   const int maxNumberofSeeds_;
0084   const int maxNumberofClusters_;
0085
0086   const int debug_;
0087
0088   PositionCalc posCalculator_;  // position calculation algorithm
0089 };
0090
0091 EgammaHLTNxNClusterProducer::EgammaHLTNxNClusterProducer(const edm::ParameterSet &ps)
0092     : doBarrel_(ps.getParameter<bool>("doBarrel")),
0093       doEndcaps_(ps.getParameter<bool>("doEndcaps")),
0094       barrelHitProducer_(consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("barrelHitProducer"))),
0095       endcapHitProducer_(consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("endcapHitProducer"))),
0096       caloGeometryToken_{esConsumes()},
0097       ecalChannelStatusToken_{esConsumes()},
0098       clusEtaSize_(ps.getParameter<int>("clusEtaSize")),
0099       clusPhiSize_(ps.getParameter<int>("clusPhiSize")),
0100       barrelClusterCollection_(ps.getParameter<std::string>("barrelClusterCollection")),
0101       endcapClusterCollection_(ps.getParameter<std::string>("endcapClusterCollection")),
0102       clusSeedThr_(ps.getParameter<double>("clusSeedThr")),
0103       clusSeedThrEndCap_(ps.getParameter<double>("clusSeedThrEndCap")),
0104       useRecoFlag_(ps.getParameter<bool>("useRecoFlag")),
0105       flagLevelRecHitsToUse_(ps.getParameter<int>("flagLevelRecHitsToUse")),
0106       useDBStatus_(ps.getParameter<bool>("useDBStatus")),
0107       statusLevelRecHitsToUse_(ps.getParameter<int>("statusLevelRecHitsToUse")),
0108       maxNumberofSeeds_(ps.getParameter<int>("maxNumberofSeeds")),
0109       maxNumberofClusters_(ps.getParameter<int>("maxNumberofClusters")),
0110       debug_(ps.getParameter<int>("debugLevel")),
0111       posCalculator_(PositionCalc(ps.getParameter<edm::ParameterSet>("posCalcParameters"))) {
0112   produces<reco::BasicClusterCollection>(barrelClusterCollection_);
0113   produces<reco::BasicClusterCollection>(endcapClusterCollection_);
0114 }
0115
0116 EgammaHLTNxNClusterProducer::~EgammaHLTNxNClusterProducer() {}
0117
0118 void EgammaHLTNxNClusterProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0119   edm::ParameterSetDescription desc;
0120   desc.add<bool>(("doBarrel"), true);
0121   desc.add<bool>(("doEndcaps"), true);
0122   desc.add<edm::InputTag>(("barrelHitProducer"), edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsEB"));
0123   desc.add<edm::InputTag>(("endcapHitProducer"), edm::InputTag("hltEcalRegionalPi0EtaRecHit", "EcalRecHitsEE"));
0124   desc.add<int>(("clusEtaSize"), 3);
0125   desc.add<int>(("clusPhiSize"), 3);
0126   desc.add<std::string>(("barrelClusterCollection"), "Simple3x3ClustersBarrel");
0127   desc.add<std::string>(("endcapClusterCollection"), "Simple3x3ClustersEndcap");
0128   desc.add<double>(("clusSeedThr"), 0.5);
0129   desc.add<double>(("clusSeedThrEndCap"), 1.0);
0130   desc.add<bool>(("useRecoFlag"), false);
0131   desc.add<int>(("flagLevelRecHitsToUse"), 1);
0132   desc.add<bool>(("useDBStatus"), true);
0133   desc.add<int>(("statusLevelRecHitsToUse"), 1);
0134
0135   edm::ParameterSetDescription posCalcPSET;
0136   posCalcPSET.add<double>("T0_barl", 7.4);
0137   posCalcPSET.add<double>("T0_endc", 3.1);
0138   posCalcPSET.add<double>("T0_endcPresh", 1.2);
0139   posCalcPSET.add<double>("W0", 4.2);
0140   posCalcPSET.add<double>("X0", 0.89);
0141   posCalcPSET.add<bool>("LogWeighted", true);
0142   desc.add<edm::ParameterSetDescription>("posCalcParameters", posCalcPSET);
0143
0144   desc.add<int>(("maxNumberofSeeds"), 1000);
0145   desc.add<int>(("maxNumberofClusters"), 200);
0146   desc.add<int>(("debugLevel"), 0);
0147   descriptions.add(("hltEgammaHLTNxNClusterProducer"), desc);
0148 }
0149
0150 void EgammaHLTNxNClusterProducer::produce(edm::Event &evt, const edm::EventSetup &eventSetup) {
0151   if (doBarrel_) {
0152     edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0153     evt.getByToken(barrelHitProducer_, barrelRecHitsHandle);
0154     if (!barrelRecHitsHandle.isValid()) {
0155       LogDebug("") << "EgammaHLTNxNClusterProducer Error! can't get product eb hit!" << std::endl;
0156     }
0157
0158     const EcalRecHitCollection *hits_eb = barrelRecHitsHandle.product();
0159     if (debug_ >= 2)
0160       LogDebug("") << "EgammaHLTNxNClusterProducer nEBrechits: " << evt.id().run() << " event " << evt.id().event()
0161                    << " " << hits_eb->size() << std::endl;
0162
0163     makeNxNClusters(evt, eventSetup, *hits_eb, reco::CaloID::DET_ECAL_BARREL);
0164   }
0165
0166   if (doEndcaps_) {
0167     edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0168     evt.getByToken(endcapHitProducer_, endcapRecHitsHandle);
0169     if (!endcapRecHitsHandle.isValid()) {
0170       LogDebug("") << "EgammaHLTNxNClusterProducer Error! can't get product ee hit!" << std::endl;
0171     }
0172
0173     const EcalRecHitCollection *hits_ee = endcapRecHitsHandle.product();
0174     if (debug_ >= 2)
0175       LogDebug("") << "EgammaHLTNxNClusterProducer nEErechits: " << evt.id().run() << " event " << evt.id().event()
0176                    << " " << hits_ee->size() << std::endl;
0177     makeNxNClusters(evt, eventSetup, *hits_ee, reco::CaloID::DET_ECAL_ENDCAP);
0178   }
0179 }
0180
0181 bool EgammaHLTNxNClusterProducer::checkStatusOfEcalRecHit(const EcalChannelStatus &channelStatus,
0182                                                           const EcalRecHit &rh) {
0183   if (useRecoFlag_) {  ///from recoFlag()
0184     int flag = rh.recoFlag();
0185     if (flagLevelRecHitsToUse_ == 0) {  ///good
0186       if (flag != 0)
0187         return false;
0188     } else if (flagLevelRecHitsToUse_ == 1) {  ///good || PoorCalib
0189       if (flag != 0 && flag != 4)
0190         return false;
0191     } else if (flagLevelRecHitsToUse_ == 2) {  ///good || PoorCalib || LeadingEdgeRecovered || kNeighboursRecovered,
0192       if (flag != 0 && flag != 4 && flag != 6 && flag != 7)
0193         return false;
0194     }
0195   }
0196   if (useDBStatus_) {  //// from DB
0197     int status = int(channelStatus[rh.id().rawId()].getStatusCode());
0198     if (status > statusLevelRecHitsToUse_)
0199       return false;
0200   }
0201
0202   return true;
0203 }
0204
0205 void EgammaHLTNxNClusterProducer::makeNxNClusters(edm::Event &evt,
0206                                                   const edm::EventSetup &eventSetup,
0207                                                   const EcalRecHitCollection &hits,
0208                                                   const reco::CaloID::Detectors detector) {
0209   ///get status from DB
0210   edm::ESHandle<EcalChannelStatus> csHandle;
0211   if (useDBStatus_) {
0212     csHandle = eventSetup.getHandle(ecalChannelStatusToken_);
0213   }
0214   const EcalChannelStatus &channelStatus = *csHandle;
0215
0216   std::vector<EcalRecHit> seeds;
0217
0218   double clusterSeedThreshold;
0219   if (detector == reco::CaloID::DET_ECAL_BARREL) {
0220     clusterSeedThreshold = clusSeedThr_;
0221   } else {
0222     clusterSeedThreshold = clusSeedThrEndCap_;
0223   }
0224
0225   for (auto const itt : hits) {
0226     double energy = itt.energy();
0227     if (!checkStatusOfEcalRecHit(channelStatus, itt))
0228       continue;
0229     if (energy > clusterSeedThreshold)
0230       seeds.push_back(itt);
0231
0232     if (int(seeds.size()) > maxNumberofSeeds_) {  //too many seeds, like beam splash events
0233       seeds.clear();                              //// empty seeds vector, don't do clustering anymore
0234       break;
0235     }
0236   }
0237
0238   // get the geometry and topology from the event setup:
0239   auto const &caloGeometry = eventSetup.getData(caloGeometryToken_);
0240
0241   const CaloSubdetectorGeometry *geometry_p;
0242   std::unique_ptr<CaloSubdetectorTopology> topology_p;
0243   if (detector == reco::CaloID::DET_ECAL_BARREL) {
0244     geometry_p = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0245     topology_p = std::make_unique<EcalBarrelTopology>(caloGeometry);
0246   } else {
0247     geometry_p = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0248     topology_p = std::make_unique<EcalEndcapTopology>(caloGeometry);
0249   }
0250
0251   const CaloSubdetectorGeometry *geometryES_p;
0252   geometryES_p = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0253
0254   std::vector<reco::BasicCluster> clusters;
0255   std::vector<DetId> usedXtals;
0256
0257   // sort seed according to Energy
0258   sort(seeds.begin(), seeds.end(), [](auto const &x, auto const &y) { return (x.energy() > y.energy()); });
0259
0260   for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
0261     DetId seed_id = itseed->id();
0262
0263     std::vector<DetId>::iterator itdet = find(usedXtals.begin(), usedXtals.end(), seed_id);
0264     if (itdet != usedXtals.end())
0265       continue;
0266
0267     std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
0268     std::vector<std::pair<DetId, float> > clus_used;
0269
0270     float clus_energy = 0;
0271
0272     for (auto const &detid : clus_v) {
0273       //not yet used
0274       std::vector<DetId>::iterator itdet = find(usedXtals.begin(), usedXtals.end(), detid);
0275       if (itdet != usedXtals.end())
0276         continue;
0277       //inside the collection
0278       EcalRecHitCollection::const_iterator hit = hits.find(detid);
0279       if (hit == hits.end())
0280         continue;
0281
0282       if (!checkStatusOfEcalRecHit(channelStatus, *hit))
0283         continue;
0284
0285       usedXtals.push_back(detid);
0286       clus_used.push_back(std::pair<DetId, float>(detid, 1.));
0287       clus_energy += hit->energy();
0288
0289     }  //// end of making one nxn simple cluster
0290
0291     if (clus_energy <= 0)
0292       continue;
0293
0294     math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used, &hits, geometry_p, geometryES_p);
0295
0296     if (debug_ >= 2)
0297       LogDebug("") << "nxn_cluster in run " << evt.id().run() << " event " << evt.id().event()
0298                    << " energy: " << clus_energy << " eta: " << clus_pos.Eta() << " phi: " << clus_pos.Phi()
0299                    << " nRecHits: " << clus_used.size() << std::endl;
0300
0301     clusters.push_back(reco::BasicCluster(
0302         clus_energy, clus_pos, reco::CaloID(detector), clus_used, reco::CaloCluster::island, seed_id));
0303     if (int(clusters.size()) > maxNumberofClusters_) {  ///if too much clusters made, then return 0 also
0304       clusters.clear();
0305       break;
0306     }
0307   }
0308
0309   //Create empty output collections
0310   auto clusters_p = std::make_unique<reco::BasicClusterCollection>();
0311   clusters_p->assign(clusters.begin(), clusters.end());
0312   if (detector == reco::CaloID::DET_ECAL_BARREL) {
0313     if (debug_ >= 1)
0314       LogDebug("") << "nxnclusterProducer: " << clusters_p->size() << " made in barrel" << std::endl;
0315     evt.put(std::move(clusters_p), barrelClusterCollection_);
0316   } else {
0317     if (debug_ >= 1)
0318       LogDebug("") << "nxnclusterProducer: " << clusters_p->size() << " made in endcap" << std::endl;
0319     evt.put(std::move(clusters_p), endcapClusterCollection_);
0320   }
0321 }
0322
0323 #include "FWCore/Framework/interface/MakerMacros.h"
0324 DEFINE_FWK_MODULE(EgammaHLTNxNClusterProducer);