EgammaClusterProducers/src/HybridClusterProducer.cc

0001 #include "CommonTools/Utils/interface/StringToEnumValue.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/EgammaReco/interface/BasicCluster.h"
0007 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0008 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0009 #include "DataFormats/Math/interface/RectangularEtaPhiRegion.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 "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0018 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0019 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0020 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0021 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0022 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0023 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0024 #include "RecoEcal/EgammaClusterAlgos/interface/HybridClusterAlgo.h"
0025 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0026 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0027 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0028
0029 #include <iostream>
0030 #include <memory>
0031 #include <vector>
0032
0033 class HybridClusterProducer : public edm::stream::EDProducer<> {
0034 public:
0035   HybridClusterProducer(const edm::ParameterSet& ps);
0036
0037   ~HybridClusterProducer() override;
0038
0039   void produce(edm::Event&, const edm::EventSetup&) override;
0040
0041 private:
0042   int nEvt_;  // internal counter of events
0043
0044   std::string basicclusterCollection_;
0045   std::string superclusterCollection_;
0046
0047   edm::EDGetTokenT<EcalRecHitCollection> hitsToken_;
0048   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geoToken_;
0049   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> sevLvToken_;
0050
0051   HybridClusterAlgo* hybrid_p;  // clustering algorithm
0052   PositionCalc posCalculator_;  // position calculation algorithm
0053 };
0054
0055 #include "FWCore/Framework/interface/MakerMacros.h"
0056 DEFINE_FWK_MODULE(HybridClusterProducer);
0057
0058 HybridClusterProducer::HybridClusterProducer(const edm::ParameterSet& ps) {
0059   basicclusterCollection_ = ps.getParameter<std::string>("basicclusterCollection");
0060   superclusterCollection_ = ps.getParameter<std::string>("superclusterCollection");
0061   hitsToken_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("recHitsCollection"));
0062   geoToken_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0063   sevLvToken_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0064
0065   //Setup for core tools objects.
0066   edm::ParameterSet posCalcParameters = ps.getParameter<edm::ParameterSet>("posCalcParameters");
0067
0068   posCalculator_ = PositionCalc(posCalcParameters);
0069
0070   const std::vector<std::string> flagnames = ps.getParameter<std::vector<std::string> >("RecHitFlagToBeExcluded");
0071
0072   const std::vector<int> flagsexcl = StringToEnumValue<EcalRecHit::Flags>(flagnames);
0073
0074   const std::vector<std::string> severitynames =
0075       ps.getParameter<std::vector<std::string> >("RecHitSeverityToBeExcluded");
0076
0077   const std::vector<int> severitiesexcl = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynames);
0078
0079   hybrid_p = new HybridClusterAlgo(ps.getParameter<double>("HybridBarrelSeedThr"),
0080                                    ps.getParameter<int>("step"),
0081                                    ps.getParameter<double>("ethresh"),
0082                                    ps.getParameter<double>("eseed"),
0083                                    ps.getParameter<double>("xi"),
0084                                    ps.getParameter<bool>("useEtForXi"),
0085                                    ps.getParameter<double>("ewing"),
0086                                    flagsexcl,
0087                                    posCalculator_,
0088                                    ps.getParameter<bool>("dynamicEThresh"),
0089                                    ps.getParameter<double>("eThreshA"),
0090                                    ps.getParameter<double>("eThreshB"),
0091                                    severitiesexcl,
0092                                    ps.getParameter<bool>("excludeFlagged"));
0093   //bremRecoveryPset,
0094
0095   // get brem recovery parameters
0096   bool dynamicPhiRoad = ps.getParameter<bool>("dynamicPhiRoad");
0097   if (dynamicPhiRoad) {
0098     edm::ParameterSet bremRecoveryPset = ps.getParameter<edm::ParameterSet>("bremRecoveryPset");
0099     hybrid_p->setDynamicPhiRoad(bremRecoveryPset);
0100   }
0101
0102   produces<reco::BasicClusterCollection>(basicclusterCollection_);
0103   produces<reco::SuperClusterCollection>(superclusterCollection_);
0104   nEvt_ = 0;
0105 }
0106
0107 HybridClusterProducer::~HybridClusterProducer() { delete hybrid_p; }
0108
0109 void HybridClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0110   // get the hit collection from the event:
0111   edm::Handle<EcalRecHitCollection> rhcHandle;
0112
0113   evt.getByToken(hitsToken_, rhcHandle);
0114   if (!(rhcHandle.isValid())) {
0115     edm::LogError("MissingProduct") << "could not get a handle on the EcalRecHitCollection!";
0116     return;
0117   }
0118   const EcalRecHitCollection* hit_collection = rhcHandle.product();
0119
0120   // get the collection geometry:
0121   edm::ESHandle<CaloGeometry> geoHandle = es.getHandle(geoToken_);
0122   const CaloGeometry& geometry = *geoHandle;
0123   const CaloSubdetectorGeometry* geometry_p;
0124   std::unique_ptr<const CaloSubdetectorTopology> topology;
0125
0126   edm::ESHandle<EcalSeverityLevelAlgo> sevLv = es.getHandle(sevLvToken_);
0127
0128   geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0129   topology = std::make_unique<EcalBarrelTopology>(*geoHandle);
0130
0131   // make the Basic clusters!
0132   reco::BasicClusterCollection basicClusters;
0133   hybrid_p->makeClusters(
0134       hit_collection, geometry_p, basicClusters, sevLv.product(), false, std::vector<RectangularEtaPhiRegion>());
0135
0136   LogTrace("EcalClusters") << "Finished clustering - BasicClusterCollection returned to producer...";
0137
0138   // create a unique_ptr to a BasicClusterCollection, copy the clusters into it and put in the Event:
0139   auto basicclusters_p = std::make_unique<reco::BasicClusterCollection>();
0140   basicclusters_p->assign(basicClusters.begin(), basicClusters.end());
0141   edm::OrphanHandle<reco::BasicClusterCollection> bccHandle =
0142       evt.put(std::move(basicclusters_p), basicclusterCollection_);
0143
0144   //Basic clusters now in the event.
0145   LogTrace("EcalClusters") << "Basic Clusters now put into event.";
0146
0147   //Weird though it is, get the BasicClusters back out of the event.  We need the
0148   //edm::Ref to these guys to make our superclusters for Hybrid.
0149
0150   if (!(bccHandle.isValid())) {
0151     edm::LogError("Missing Product") << "could not get a handle on the BasicClusterCollection!";
0152     return;
0153   }
0154
0155   const reco::BasicClusterCollection& clusterCollection = *bccHandle;
0156
0157   LogTrace("EcalClusters") << "Got the BasicClusterCollection" << std::endl;
0158
0159   reco::CaloClusterPtrVector clusterPtrVector;
0160   for (unsigned int i = 0; i < clusterCollection.size(); i++) {
0161     clusterPtrVector.push_back(reco::CaloClusterPtr(bccHandle, i));
0162   }
0163
0164   reco::SuperClusterCollection superClusters = hybrid_p->makeSuperClusters(clusterPtrVector);
0165   LogTrace("EcalClusters") << "Found: " << superClusters.size() << " superclusters.";
0166
0167   auto superclusters_p = std::make_unique<reco::SuperClusterCollection>();
0168   superclusters_p->assign(superClusters.begin(), superClusters.end());
0169
0170   evt.put(std::move(superclusters_p), superclusterCollection_);
0171   LogTrace("EcalClusters") << "Hybrid Clusters (Basic/Super) added to the Event! :-)";
0172
0173   nEvt_++;
0174 }