EgammaClusterProducers/src/IslandClusterProducer.cc

0001 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0002 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0005 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0006 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0007 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0008 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0009 #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
0010 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0011 #include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/Utilities/interface/ESGetToken.h"
0020 #include "FWCore/Utilities/interface/Exception.h"
0021 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0022 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0023 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0024 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0025 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0026 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0028 #include "RecoEcal/EgammaClusterAlgos/interface/IslandClusterAlgo.h"
0029 #include "RecoEcal/EgammaCoreTools/interface/ClusterShapeAlgo.h"
0030 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0031
0032 #include <ctime>
0033 #include <iostream>
0034 #include <memory>
0035 #include <vector>
0036
0037 class IslandClusterProducer : public edm::stream::EDProducer<> {
0038 public:
0039   IslandClusterProducer(const edm::ParameterSet& ps);
0040
0041   ~IslandClusterProducer() override;
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043   void produce(edm::Event&, const edm::EventSetup&) override;
0044
0045 private:
0046   int nMaxPrintout_;  // max # of printouts
0047   int nEvt_;          // internal counter of events
0048
0049   IslandClusterAlgo::VerbosityLevel verbosity;
0050
0051   edm::EDGetTokenT<EcalRecHitCollection> barrelRecHits_;
0052   edm::EDGetTokenT<EcalRecHitCollection> endcapRecHits_;
0053   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0054
0055   std::string barrelClusterCollection_;
0056   std::string endcapClusterCollection_;
0057
0058   std::string clustershapecollectionEB_;
0059   std::string clustershapecollectionEE_;
0060
0061   //BasicClusterShape AssociationMap
0062   std::string barrelClusterShapeAssociation_;
0063   std::string endcapClusterShapeAssociation_;
0064
0065   PositionCalc posCalculator_;  // position calculation algorithm
0066   ClusterShapeAlgo shapeAlgo_;  // cluster shape algorithm
0067   IslandClusterAlgo* island_p;
0068
0069   bool counterExceeded() const { return ((nEvt_ > nMaxPrintout_) || (nMaxPrintout_ < 0)); }
0070
0071   const EcalRecHitCollection* getCollection(edm::Event& evt, const edm::EDGetTokenT<EcalRecHitCollection>& token);
0072
0073   void clusterizeECALPart(edm::Event& evt,
0074                           const edm::EventSetup& es,
0075                           const edm::EDGetTokenT<EcalRecHitCollection>& token,
0076                           const std::string& clusterCollection,
0077                           const std::string& clusterShapeAssociation,
0078                           const IslandClusterAlgo::EcalPart& ecalPart);
0079
0080   void outputValidationInfo(reco::CaloClusterPtrVector& clusterPtrVector);
0081 };
0082
0083 #include "FWCore/Framework/interface/MakerMacros.h"
0084 DEFINE_FWK_MODULE(IslandClusterProducer);
0085
0086 IslandClusterProducer::IslandClusterProducer(const edm::ParameterSet& ps) {
0087   // The verbosity level
0088   std::string verbosityString = ps.getParameter<std::string>("VerbosityLevel");
0089   if (verbosityString == "DEBUG")
0090     verbosity = IslandClusterAlgo::pDEBUG;
0091   else if (verbosityString == "WARNING")
0092     verbosity = IslandClusterAlgo::pWARNING;
0093   else if (verbosityString == "INFO")
0094     verbosity = IslandClusterAlgo::pINFO;
0095   else
0096     verbosity = IslandClusterAlgo::pERROR;
0097
0098   // Parameters to identify the hit collections
0099   barrelRecHits_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("barrelHits"));
0100   endcapRecHits_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("endcapHits"));
0101
0102   //EventSetup Token for CaloGeometry
0103   caloGeometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0104
0105   // The names of the produced cluster collections
0106   barrelClusterCollection_ = ps.getParameter<std::string>("barrelClusterCollection");
0107   endcapClusterCollection_ = ps.getParameter<std::string>("endcapClusterCollection");
0108
0109   // Island algorithm parameters
0110   double barrelSeedThreshold = ps.getParameter<double>("IslandBarrelSeedThr");
0111   double endcapSeedThreshold = ps.getParameter<double>("IslandEndcapSeedThr");
0112
0113   // Parameters for the position calculation:
0114   edm::ParameterSet posCalcParameters = ps.getParameter<edm::ParameterSet>("posCalcParameters");
0115   posCalculator_ = PositionCalc(posCalcParameters);
0116   shapeAlgo_ = ClusterShapeAlgo(posCalcParameters);
0117
0118   clustershapecollectionEB_ = ps.getParameter<std::string>("clustershapecollectionEB");
0119   clustershapecollectionEE_ = ps.getParameter<std::string>("clustershapecollectionEE");
0120
0121   //AssociationMap
0122   barrelClusterShapeAssociation_ = ps.getParameter<std::string>("barrelShapeAssociation");
0123   endcapClusterShapeAssociation_ = ps.getParameter<std::string>("endcapShapeAssociation");
0124
0125   const std::vector<std::string> seedflagnamesEB =
0126       ps.getParameter<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEB");
0127   const std::vector<int> seedflagsexclEB = StringToEnumValue<EcalRecHit::Flags>(seedflagnamesEB);
0128
0129   const std::vector<std::string> seedflagnamesEE =
0130       ps.getParameter<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEE");
0131   const std::vector<int> seedflagsexclEE = StringToEnumValue<EcalRecHit::Flags>(seedflagnamesEE);
0132
0133   const std::vector<std::string> flagnamesEB = ps.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEB");
0134   const std::vector<int> flagsexclEB = StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
0135
0136   const std::vector<std::string> flagnamesEE = ps.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEE");
0137   const std::vector<int> flagsexclEE = StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
0138
0139   // Produces a collection of barrel and a collection of endcap clusters
0140
0141   produces<reco::ClusterShapeCollection>(clustershapecollectionEE_);
0142   produces<reco::BasicClusterCollection>(endcapClusterCollection_);
0143   produces<reco::ClusterShapeCollection>(clustershapecollectionEB_);
0144   produces<reco::BasicClusterCollection>(barrelClusterCollection_);
0145   produces<reco::BasicClusterShapeAssociationCollection>(barrelClusterShapeAssociation_);
0146   produces<reco::BasicClusterShapeAssociationCollection>(endcapClusterShapeAssociation_);
0147
0148   island_p = new IslandClusterAlgo(barrelSeedThreshold,
0149                                    endcapSeedThreshold,
0150                                    posCalculator_,
0151                                    seedflagsexclEB,
0152                                    seedflagsexclEE,
0153                                    flagsexclEB,
0154                                    flagsexclEE,
0155                                    verbosity);
0156
0157   nEvt_ = 0;
0158 }
0159
0160 IslandClusterProducer::~IslandClusterProducer() { delete island_p; }
0161
0162 void IslandClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0163   edm::ParameterSetDescription desc;
0164   desc.add<std::string>("VerbosityLevel", "ERROR");
0165   desc.add<edm::InputTag>("barrelHits", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0166   desc.add<edm::InputTag>("endcapHits", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0167   desc.add<std::string>("barrelClusterCollection", "islandBarrelBasicClusters");
0168   desc.add<std::string>("endcapClusterCollection", "islandEndcapBasicClusters");
0169   desc.add<double>("IslandBarrelSeedThr", 0.5);
0170   desc.add<double>("IslandEndcapSeedThr", 0.18);
0171
0172   edm::ParameterSetDescription posCalcParameters;
0173   posCalcParameters.add<bool>("LogWeighted", true);
0174   posCalcParameters.add<double>("T0_barl", 7.4);
0175   posCalcParameters.add<double>("T0_endc", 3.1);
0176   posCalcParameters.add<double>("T0_endcPresh", 1.2);
0177   posCalcParameters.add<double>("W0", 4.2);
0178   posCalcParameters.add<double>("X0", 0.89);
0179   desc.add<edm::ParameterSetDescription>("posCalcParameters", posCalcParameters);
0180
0181   desc.add<std::string>("clustershapecollectionEE", "islandEndcapShape");
0182   desc.add<std::string>("clustershapecollectionEB", "islandBarrelShape");
0183   desc.add<std::string>("barrelShapeAssociation", "islandBarrelShapeAssoc");
0184   desc.add<std::string>("endcapShapeAssociation", "islandEndcapShapeAssoc");
0185   desc.add<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEB", {});
0186   desc.add<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEE", {});
0187   desc.add<std::vector<std::string>>("RecHitFlagToBeExcludedEB", {});
0188   desc.add<std::vector<std::string>>("RecHitFlagToBeExcludedEE", {});
0189   descriptions.add("IslandClusterProducer", desc);
0190 }
0191
0192 void IslandClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0193   clusterizeECALPart(
0194       evt, es, endcapRecHits_, endcapClusterCollection_, endcapClusterShapeAssociation_, IslandClusterAlgo::endcap);
0195   clusterizeECALPart(
0196       evt, es, barrelRecHits_, barrelClusterCollection_, barrelClusterShapeAssociation_, IslandClusterAlgo::barrel);
0197   nEvt_++;
0198 }
0199
0200 const EcalRecHitCollection* IslandClusterProducer::getCollection(edm::Event& evt,
0201                                                                  const edm::EDGetTokenT<EcalRecHitCollection>& token) {
0202   edm::Handle<EcalRecHitCollection> rhcHandle;
0203   evt.getByToken(token, rhcHandle);
0204   return rhcHandle.product();
0205 }
0206
0207 void IslandClusterProducer::clusterizeECALPart(edm::Event& evt,
0208                                                const edm::EventSetup& es,
0209                                                const edm::EDGetTokenT<EcalRecHitCollection>& token,
0210                                                const std::string& clusterCollection,
0211                                                const std::string& clusterShapeAssociation,
0212                                                const IslandClusterAlgo::EcalPart& ecalPart) {
0213   // get the hit collection from the event:
0214   const EcalRecHitCollection* hitCollection_p = getCollection(evt, token);
0215
0216   // get the geometry and topology from the event setup:
0217   edm::ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0218
0219   const CaloSubdetectorGeometry* geometry_p;
0220   std::unique_ptr<CaloSubdetectorTopology> topology_p;
0221
0222   std::string clustershapetag;
0223   if (ecalPart == IslandClusterAlgo::barrel) {
0224     geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0225     topology_p = std::make_unique<EcalBarrelTopology>(*geoHandle);
0226   } else {
0227     geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0228     topology_p = std::make_unique<EcalEndcapTopology>(*geoHandle);
0229   }
0230
0231   const CaloSubdetectorGeometry* geometryES_p;
0232   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0233
0234   // Run the clusterization algorithm:
0235   reco::BasicClusterCollection clusters;
0236   clusters = island_p->makeClusters(hitCollection_p, geometry_p, topology_p.get(), geometryES_p, ecalPart);
0237
0238   //Create associated ClusterShape objects.
0239   std::vector<reco::ClusterShape> ClusVec;
0240   for (int erg = 0; erg < int(clusters.size()); ++erg) {
0241     reco::ClusterShape TestShape = shapeAlgo_.Calculate(clusters[erg], hitCollection_p, geometry_p, topology_p.get());
0242     ClusVec.push_back(TestShape);
0243   }
0244
0245   //Put clustershapes in event, but retain a Handle on them.
0246   auto clustersshapes_p = std::make_unique<reco::ClusterShapeCollection>();
0247   clustersshapes_p->assign(ClusVec.begin(), ClusVec.end());
0248   edm::OrphanHandle<reco::ClusterShapeCollection> clusHandle;
0249   if (ecalPart == IslandClusterAlgo::barrel)
0250     clusHandle = evt.put(std::move(clustersshapes_p), clustershapecollectionEB_);
0251   else
0252     clusHandle = evt.put(std::move(clustersshapes_p), clustershapecollectionEE_);
0253
0254   // create a unique_ptr to a BasicClusterCollection, copy the barrel clusters into it and put in the Event:
0255   auto clusters_p = std::make_unique<reco::BasicClusterCollection>();
0256   clusters_p->assign(clusters.begin(), clusters.end());
0257   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle;
0258   if (ecalPart == IslandClusterAlgo::barrel)
0259     bccHandle = evt.put(std::move(clusters_p), barrelClusterCollection_);
0260   else
0261     bccHandle = evt.put(std::move(clusters_p), endcapClusterCollection_);
0262
0263   // BasicClusterShapeAssociationMap
0264   auto shapeAssocs_p = std::make_unique<reco::BasicClusterShapeAssociationCollection>(bccHandle, clusHandle);
0265   for (unsigned int i = 0; i < clusHandle->size(); i++) {
0266     shapeAssocs_p->insert(edm::Ref<reco::BasicClusterCollection>(bccHandle, i),
0267                           edm::Ref<reco::ClusterShapeCollection>(clusHandle, i));
0268   }
0269   evt.put(std::move(shapeAssocs_p), clusterShapeAssociation);
0270 }