EgammaClusterProducers/src/Multi5x5ClusterProducer.cc

0001 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0002 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0003 #include "DataFormats/CaloRecHit/interface/CaloID.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0006 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0007 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0008 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0009 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.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 "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0019 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0020 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0021 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0022 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0023 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "RecoEcal/EgammaClusterAlgos/interface/Multi5x5ClusterAlgo.h"
0026 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0027
0028 #include <ctime>
0029 #include <iostream>
0030 #include <memory>
0031 #include <vector>
0032
0033 class Multi5x5ClusterProducer : public edm::stream::EDProducer<> {
0034 public:
0035   Multi5x5ClusterProducer(const edm::ParameterSet& ps);
0036
0037   ~Multi5x5ClusterProducer() override;
0038
0039   void produce(edm::Event&, const edm::EventSetup&) override;
0040
0041 private:
0042   int nMaxPrintout_;  // max # of printouts
0043   int nEvt_;          // internal counter of events
0044
0045   // cluster which regions
0046   bool doBarrel_;
0047   bool doEndcap_;
0048
0049   edm::EDGetTokenT<EcalRecHitCollection> barrelHitToken_;
0050   edm::EDGetTokenT<EcalRecHitCollection> endcapHitToken_;
0051   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0052
0053   std::string barrelClusterCollection_;
0054   std::string endcapClusterCollection_;
0055
0056   PositionCalc posCalculator_;  // position calculation algorithm
0057   Multi5x5ClusterAlgo* island_p;
0058
0059   bool counterExceeded() const { return ((nEvt_ > nMaxPrintout_) || (nMaxPrintout_ < 0)); }
0060
0061   const EcalRecHitCollection* getCollection(edm::Event& evt, const edm::EDGetTokenT<EcalRecHitCollection>& token);
0062
0063   void clusterizeECALPart(edm::Event& evt,
0064                           const edm::EventSetup& es,
0065                           const edm::EDGetTokenT<EcalRecHitCollection>& token,
0066                           const std::string& clusterCollection,
0067                           const reco::CaloID::Detectors detector);
0068
0069   void outputValidationInfo(reco::CaloClusterPtrVector& clusterPtrVector);
0070 };
0071
0072 #include "FWCore/Framework/interface/MakerMacros.h"
0073 DEFINE_FWK_MODULE(Multi5x5ClusterProducer);
0074
0075 Multi5x5ClusterProducer::Multi5x5ClusterProducer(const edm::ParameterSet& ps) {
0076   // Parameters to identify the hit collections
0077   barrelHitToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("barrelHitTag"));
0078
0079   endcapHitToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("endcapHitTag"));
0080
0081   //EventSetup Token for CaloGeometry
0082   caloGeometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0083
0084   // should cluster algo be run in barrel and endcap?
0085   doEndcap_ = ps.getParameter<bool>("doEndcap");
0086   doBarrel_ = ps.getParameter<bool>("doBarrel");
0087
0088   // The names of the produced cluster collections
0089   barrelClusterCollection_ = ps.getParameter<std::string>("barrelClusterCollection");
0090   endcapClusterCollection_ = ps.getParameter<std::string>("endcapClusterCollection");
0091
0092   // Island algorithm parameters
0093   double barrelSeedThreshold = ps.getParameter<double>("IslandBarrelSeedThr");
0094   double endcapSeedThreshold = ps.getParameter<double>("IslandEndcapSeedThr");
0095
0096   const std::vector<std::string> flagnames = ps.getParameter<std::vector<std::string> >("RecHitFlagToBeExcluded");
0097
0098   const std::vector<int> v_chstatus = StringToEnumValue<EcalRecHit::Flags>(flagnames);
0099
0100   // Parameters for the position calculation:
0101   edm::ParameterSet posCalcParameters = ps.getParameter<edm::ParameterSet>("posCalcParameters");
0102   posCalculator_ = PositionCalc(posCalcParameters);
0103
0104   // Produces a collection of barrel and a collection of endcap clusters
0105   produces<reco::BasicClusterCollection>(endcapClusterCollection_);
0106   produces<reco::BasicClusterCollection>(barrelClusterCollection_);
0107
0108   bool reassignSeedCrysToClusterItSeeds = false;
0109   if (ps.exists("reassignSeedCrysToClusterItSeeds"))
0110     reassignSeedCrysToClusterItSeeds = ps.getParameter<bool>("reassignSeedCrysToClusterItSeeds");
0111
0112   island_p = new Multi5x5ClusterAlgo(
0113       barrelSeedThreshold, endcapSeedThreshold, v_chstatus, posCalculator_, reassignSeedCrysToClusterItSeeds);
0114
0115   nEvt_ = 0;
0116 }
0117
0118 Multi5x5ClusterProducer::~Multi5x5ClusterProducer() { delete island_p; }
0119
0120 void Multi5x5ClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0121   if (doEndcap_) {
0122     clusterizeECALPart(evt, es, endcapHitToken_, endcapClusterCollection_, reco::CaloID::DET_ECAL_ENDCAP);
0123   }
0124   if (doBarrel_) {
0125     clusterizeECALPart(evt, es, barrelHitToken_, barrelClusterCollection_, reco::CaloID::DET_ECAL_BARREL);
0126   }
0127
0128   nEvt_++;
0129 }
0130
0131 const EcalRecHitCollection* Multi5x5ClusterProducer::getCollection(
0132     edm::Event& evt, const edm::EDGetTokenT<EcalRecHitCollection>& token) {
0133   edm::Handle<EcalRecHitCollection> rhcHandle;
0134   evt.getByToken(token, rhcHandle);
0135   return rhcHandle.product();
0136 }
0137
0138 void Multi5x5ClusterProducer::clusterizeECALPart(edm::Event& evt,
0139                                                  const edm::EventSetup& es,
0140                                                  const edm::EDGetTokenT<EcalRecHitCollection>& token,
0141                                                  const std::string& clusterCollection,
0142                                                  const reco::CaloID::Detectors detector) {
0143   // get the hit collection from the event:
0144   const EcalRecHitCollection* hitCollection_p = getCollection(evt, token);
0145
0146   // get the geometry and topology from the event setup:
0147   edm::ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0148
0149   const CaloSubdetectorGeometry* geometry_p;
0150   std::unique_ptr<CaloSubdetectorTopology> topology_p;
0151
0152   if (detector == reco::CaloID::DET_ECAL_BARREL) {
0153     geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0154     topology_p = std::make_unique<EcalBarrelTopology>(*geoHandle);
0155   } else {
0156     geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0157     topology_p = std::make_unique<EcalEndcapTopology>(*geoHandle);
0158   }
0159
0160   const CaloSubdetectorGeometry* geometryES_p;
0161   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0162
0163   // Run the clusterization algorithm:
0164   reco::BasicClusterCollection clusters;
0165   clusters = island_p->makeClusters(hitCollection_p, geometry_p, topology_p.get(), geometryES_p, detector);
0166
0167   // create a unique_ptr to a BasicClusterCollection, copy the barrel clusters into it and put in the Event:
0168   auto clusters_p = std::make_unique<reco::BasicClusterCollection>();
0169   clusters_p->assign(clusters.begin(), clusters.end());
0170   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle;
0171   if (detector == reco::CaloID::DET_ECAL_BARREL)
0172     bccHandle = evt.put(std::move(clusters_p), barrelClusterCollection_);
0173   else
0174     bccHandle = evt.put(std::move(clusters_p), endcapClusterCollection_);
0175 }