EgammaClusterProducers/src/CosmicClusterProducer.cc

0001 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0004 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0005 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.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/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/ESGetToken.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0021 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0022 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0023 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0024 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0025 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0026 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0027 #include "RecoEcal/EgammaClusterAlgos/interface/CosmicClusterAlgo.h"
0028 #include "RecoEcal/EgammaCoreTools/interface/ClusterShapeAlgo.h"
0029 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0030
0031 #include <ctime>
0032 #include <iostream>
0033 #include <memory>
0034 #include <vector>
0035
0036 class CosmicClusterProducer : public edm::stream::EDProducer<> {
0037 public:
0038   CosmicClusterProducer(const edm::ParameterSet& ps);
0039
0040   ~CosmicClusterProducer() override;
0041
0042   void produce(edm::Event&, const edm::EventSetup&) override;
0043
0044 private:
0045   int nMaxPrintout_;  // max # of printouts
0046   int nEvt_;          // internal counter of events
0047
0048   CosmicClusterAlgo::VerbosityLevel verbosity;
0049
0050   edm::EDGetTokenT<EcalRecHitCollection> ebHitsToken_;
0051   edm::EDGetTokenT<EcalRecHitCollection> eeHitsToken_;
0052
0053   edm::EDGetTokenT<EcalUncalibratedRecHitCollection> ebUHitsToken_;
0054   edm::EDGetTokenT<EcalUncalibratedRecHitCollection> eeUHitsToken_;
0055   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0056
0057   std::string barrelClusterCollection_;
0058   std::string endcapClusterCollection_;
0059
0060   std::string clustershapecollectionEB_;
0061   std::string clustershapecollectionEE_;
0062
0063   //BasicClusterShape AssociationMap
0064   std::string barrelClusterShapeAssociation_;
0065   std::string endcapClusterShapeAssociation_;
0066
0067   PositionCalc posCalculator_;  // position calculation algorithm
0068   ClusterShapeAlgo shapeAlgo_;  // cluster shape algorithm
0069   CosmicClusterAlgo* island_p;
0070
0071   bool counterExceeded() const { return ((nEvt_ > nMaxPrintout_) || (nMaxPrintout_ < 0)); }
0072
0073   void clusterizeECALPart(edm::Event& evt,
0074                           const edm::EventSetup& es,
0075                           const edm::EDGetTokenT<EcalRecHitCollection>& hitsToken,
0076                           const edm::EDGetTokenT<EcalUncalibratedRecHitCollection>& uhitsToken,
0077                           const std::string& clusterCollection,
0078                           const std::string& clusterShapeAssociation,
0079                           const CosmicClusterAlgo::EcalPart& ecalPart);
0080
0081   void outputValidationInfo(reco::CaloClusterPtrVector& clusterPtrVector);
0082 };
0083
0084 #include "FWCore/Framework/interface/MakerMacros.h"
0085 DEFINE_FWK_MODULE(CosmicClusterProducer);
0086
0087 CosmicClusterProducer::CosmicClusterProducer(const edm::ParameterSet& ps) {
0088   // The verbosity level
0089   std::string verbosityString = ps.getParameter<std::string>("VerbosityLevel");
0090   if (verbosityString == "DEBUG")
0091     verbosity = CosmicClusterAlgo::pDEBUG;
0092   else if (verbosityString == "WARNING")
0093     verbosity = CosmicClusterAlgo::pWARNING;
0094   else if (verbosityString == "INFO")
0095     verbosity = CosmicClusterAlgo::pINFO;
0096   else
0097     verbosity = CosmicClusterAlgo::pERROR;
0098
0099   // Parameters to identify the hit collections
0100   ebHitsToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("barrelHits"));
0101   eeHitsToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("endcapHits"));
0102
0103   ebUHitsToken_ = consumes<EcalUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("barrelUncalibHits"));
0104   eeUHitsToken_ = consumes<EcalUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("endcapUncalibHits"));
0105
0106   caloGeometryToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0107
0108   // The names of the produced cluster collections
0109   barrelClusterCollection_ = ps.getParameter<std::string>("barrelClusterCollection");
0110   endcapClusterCollection_ = ps.getParameter<std::string>("endcapClusterCollection");
0111
0112   // Island algorithm parameters
0113   double barrelSeedThreshold = ps.getParameter<double>("BarrelSeedThr");
0114   double barrelSingleThreshold = ps.getParameter<double>("BarrelSingleThr");
0115   double barrelSecondThreshold = ps.getParameter<double>("BarrelSecondThr");
0116   double barrelSupThreshold = ps.getParameter<double>("BarrelSupThr");
0117   double endcapSeedThreshold = ps.getParameter<double>("EndcapSeedThr");
0118   double endcapSingleThreshold = ps.getParameter<double>("EndcapSingleThr");
0119   double endcapSecondThreshold = ps.getParameter<double>("EndcapSecondThr");
0120   double endcapSupThreshold = ps.getParameter<double>("EndcapSupThr");
0121
0122   // Parameters for the position calculation:
0123   edm::ParameterSet posCalcParameters = ps.getParameter<edm::ParameterSet>("posCalcParameters");
0124
0125   posCalculator_ = PositionCalc(posCalcParameters);
0126   shapeAlgo_ = ClusterShapeAlgo(posCalcParameters);
0127
0128   clustershapecollectionEB_ = ps.getParameter<std::string>("clustershapecollectionEB");
0129   clustershapecollectionEE_ = ps.getParameter<std::string>("clustershapecollectionEE");
0130
0131   //AssociationMap
0132   barrelClusterShapeAssociation_ = ps.getParameter<std::string>("barrelShapeAssociation");
0133   endcapClusterShapeAssociation_ = ps.getParameter<std::string>("endcapShapeAssociation");
0134
0135   // Produces a collection of barrel and a collection of endcap clusters
0136
0137   produces<reco::ClusterShapeCollection>(clustershapecollectionEE_);
0138   produces<reco::BasicClusterCollection>(endcapClusterCollection_);
0139   produces<reco::ClusterShapeCollection>(clustershapecollectionEB_);
0140   produces<reco::BasicClusterCollection>(barrelClusterCollection_);
0141   produces<reco::BasicClusterShapeAssociationCollection>(barrelClusterShapeAssociation_);
0142   produces<reco::BasicClusterShapeAssociationCollection>(endcapClusterShapeAssociation_);
0143
0144   island_p = new CosmicClusterAlgo(barrelSeedThreshold,
0145                                    barrelSingleThreshold,
0146                                    barrelSecondThreshold,
0147                                    barrelSupThreshold,
0148                                    endcapSeedThreshold,
0149                                    endcapSingleThreshold,
0150                                    endcapSecondThreshold,
0151                                    endcapSupThreshold,
0152                                    posCalculator_,
0153                                    verbosity);
0154
0155   nEvt_ = 0;
0156 }
0157
0158 CosmicClusterProducer::~CosmicClusterProducer() { delete island_p; }
0159
0160 void CosmicClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0161   clusterizeECALPart(evt,
0162                      es,
0163                      eeHitsToken_,
0164                      eeUHitsToken_,
0165                      endcapClusterCollection_,
0166                      endcapClusterShapeAssociation_,
0167                      CosmicClusterAlgo::endcap);
0168   clusterizeECALPart(evt,
0169                      es,
0170                      eeHitsToken_,
0171                      eeUHitsToken_,
0172                      barrelClusterCollection_,
0173                      barrelClusterShapeAssociation_,
0174                      CosmicClusterAlgo::barrel);
0175   nEvt_++;
0176 }
0177
0178 void CosmicClusterProducer::clusterizeECALPart(edm::Event& evt,
0179                                                const edm::EventSetup& es,
0180                                                const edm::EDGetTokenT<EcalRecHitCollection>& hitsToken,
0181                                                const edm::EDGetTokenT<EcalUncalibratedRecHitCollection>& uhitsToken,
0182                                                const std::string& clusterCollection,
0183                                                const std::string& clusterShapeAssociation,
0184                                                const CosmicClusterAlgo::EcalPart& ecalPart) {
0185   // get the hit collection from the event:
0186
0187   edm::Handle<EcalRecHitCollection> hits_h;
0188   edm::Handle<EcalUncalibratedRecHitCollection> uhits_h;
0189
0190   evt.getByToken(hitsToken, hits_h);
0191   evt.getByToken(uhitsToken, uhits_h);
0192
0193   const EcalRecHitCollection* hitCollection_p = hits_h.product();
0194   const EcalUncalibratedRecHitCollection* uhitCollection_p = uhits_h.product();
0195
0196   // get the geometry and topology from the event setup:
0197   edm::ESHandle<CaloGeometry> geoHandle = es.getHandle(caloGeometryToken_);
0198
0199   const CaloSubdetectorGeometry* geometry_p;
0200   std::unique_ptr<CaloSubdetectorTopology> topology_p;
0201
0202   std::string clustershapetag;
0203   if (ecalPart == CosmicClusterAlgo::barrel) {
0204     geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0205     topology_p = std::make_unique<EcalBarrelTopology>(*geoHandle);
0206   } else {
0207     geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0208     topology_p = std::make_unique<EcalEndcapTopology>(*geoHandle);
0209   }
0210
0211   const CaloSubdetectorGeometry* geometryES_p;
0212   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0213
0214   // Run the clusterization algorithm:
0215   reco::BasicClusterCollection clusters;
0216   clusters =
0217       island_p->makeClusters(hitCollection_p, uhitCollection_p, geometry_p, topology_p.get(), geometryES_p, ecalPart);
0218
0219   //Create associated ClusterShape objects.
0220   std::vector<reco::ClusterShape> ClusVec;
0221
0222   for (int erg = 0; erg < int(clusters.size()); ++erg) {
0223     reco::ClusterShape TestShape = shapeAlgo_.Calculate(clusters[erg], hitCollection_p, geometry_p, topology_p.get());
0224     ClusVec.push_back(TestShape);
0225   }
0226
0227   //Put clustershapes in event, but retain a Handle on them.
0228   auto clustersshapes_p = std::make_unique<reco::ClusterShapeCollection>();
0229   clustersshapes_p->assign(ClusVec.begin(), ClusVec.end());
0230   edm::OrphanHandle<reco::ClusterShapeCollection> clusHandle;
0231   if (ecalPart == CosmicClusterAlgo::barrel)
0232     clusHandle = evt.put(std::move(clustersshapes_p), clustershapecollectionEB_);
0233   else
0234     clusHandle = evt.put(std::move(clustersshapes_p), clustershapecollectionEE_);
0235
0236   // create a unique_ptr to a BasicClusterCollection, copy the barrel clusters into it and put in the Event:
0237   auto clusters_p = std::make_unique<reco::BasicClusterCollection>();
0238   clusters_p->assign(clusters.begin(), clusters.end());
0239   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle;
0240
0241   if (ecalPart == CosmicClusterAlgo::barrel)
0242     bccHandle = evt.put(std::move(clusters_p), barrelClusterCollection_);
0243   else
0244     bccHandle = evt.put(std::move(clusters_p), endcapClusterCollection_);
0245
0246   // BasicClusterShapeAssociationMap
0247   auto shapeAssocs_p = std::make_unique<reco::BasicClusterShapeAssociationCollection>(bccHandle, clusHandle);
0248
0249   for (unsigned int i = 0; i < clusHandle->size(); i++) {
0250     shapeAssocs_p->insert(edm::Ref<reco::BasicClusterCollection>(bccHandle, i),
0251                           edm::Ref<reco::ClusterShapeCollection>(clusHandle, i));
0252   }
0253   evt.put(std::move(shapeAssocs_p), clusterShapeAssociation);
0254 }