EgammaClusterProducers/src/Multi5x5SuperClusterProducer.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0003 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0004 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0005 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "RecoEcal/EgammaClusterAlgos/interface/Multi5x5BremRecoveryClusterAlgo.h"
0013
0014 #include <iostream>
0015 #include <memory>
0016 #include <sstream>
0017 #include <vector>
0018
0019 class Multi5x5SuperClusterProducer : public edm::stream::EDProducer<> {
0020 public:
0021   Multi5x5SuperClusterProducer(const edm::ParameterSet& ps);
0022
0023   void produce(edm::Event&, const edm::EventSetup&) override;
0024   void endStream() override;
0025
0026 private:
0027   edm::EDGetTokenT<reco::BasicClusterCollection> eeClustersToken_;
0028   edm::EDGetTokenT<reco::BasicClusterCollection> ebClustersToken_;
0029   edm::EDPutTokenT<reco::SuperClusterCollection> endcapPutToken_;
0030   edm::EDPutTokenT<reco::SuperClusterCollection> barrelPutToken_;
0031
0032   const float barrelEtaSearchRoad_;
0033   const float barrelPhiSearchRoad_;
0034   const float endcapEtaSearchRoad_;
0035   const float endcapPhiSearchRoad_;
0036   const float seedTransverseEnergyThreshold_;
0037
0038   const bool doBarrel_;
0039   const bool doEndcaps_;
0040
0041   std::unique_ptr<Multi5x5BremRecoveryClusterAlgo> bremAlgo_p;
0042
0043   double totalE;
0044   int noSuperClusters;
0045
0046   reco::CaloClusterPtrVector getClusterPtrVector(
0047       edm::Event& evt, const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken) const;
0048
0049   void produceSuperclustersForECALPart(edm::Event& evt,
0050                                        const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken,
0051                                        const edm::EDPutTokenT<reco::SuperClusterCollection>& putToken);
0052
0053   void outputValidationInfo(reco::SuperClusterCollection& superclusterCollection);
0054 };
0055
0056 #include "FWCore/Framework/interface/MakerMacros.h"
0057 DEFINE_FWK_MODULE(Multi5x5SuperClusterProducer);
0058
0059 Multi5x5SuperClusterProducer::Multi5x5SuperClusterProducer(const edm::ParameterSet& ps)
0060     : barrelEtaSearchRoad_{static_cast<float>(ps.getParameter<double>("barrelEtaSearchRoad"))},
0061       barrelPhiSearchRoad_{static_cast<float>(ps.getParameter<double>("barrelPhiSearchRoad"))},
0062       endcapEtaSearchRoad_{static_cast<float>(ps.getParameter<double>("endcapEtaSearchRoad"))},
0063       endcapPhiSearchRoad_{static_cast<float>(ps.getParameter<double>("endcapPhiSearchRoad"))},
0064       seedTransverseEnergyThreshold_{static_cast<float>(ps.getParameter<double>("seedTransverseEnergyThreshold"))},
0065       doBarrel_{ps.getParameter<bool>("doBarrel")},
0066       doEndcaps_{ps.getParameter<bool>("doEndcaps")},
0067       totalE{0.},
0068       noSuperClusters{0} {
0069   if (doEndcaps_) {
0070     eeClustersToken_ = consumes<reco::BasicClusterCollection>(ps.getParameter<edm::InputTag>("endcapClusterTag"));
0071   }
0072   if (doBarrel_) {
0073     ebClustersToken_ = consumes<reco::BasicClusterCollection>(ps.getParameter<edm::InputTag>("barrelClusterTag"));
0074   }
0075
0076   const edm::ParameterSet bremRecoveryPset = ps.getParameter<edm::ParameterSet>("bremRecoveryPset");
0077   bool dynamicPhiRoad = ps.getParameter<bool>("dynamicPhiRoad");
0078
0079   bremAlgo_p = std::make_unique<Multi5x5BremRecoveryClusterAlgo>(bremRecoveryPset,
0080                                                                  barrelEtaSearchRoad_,
0081                                                                  barrelPhiSearchRoad_,
0082                                                                  endcapEtaSearchRoad_,
0083                                                                  endcapPhiSearchRoad_,
0084                                                                  dynamicPhiRoad,
0085                                                                  seedTransverseEnergyThreshold_);
0086
0087   if (doEndcaps_) {
0088     endcapPutToken_ =
0089         produces<reco::SuperClusterCollection>(ps.getParameter<std::string>("endcapSuperclusterCollection"));
0090   }
0091   if (doBarrel_) {
0092     barrelPutToken_ =
0093         produces<reco::SuperClusterCollection>(ps.getParameter<std::string>("barrelSuperclusterCollection"));
0094   }
0095 }
0096
0097 void Multi5x5SuperClusterProducer::endStream() {
0098   double averEnergy = 0.;
0099   std::ostringstream str;
0100   str << "Multi5x5SuperClusterProducer::endJob()\n"
0101       << "  total # reconstructed super clusters: " << noSuperClusters << "\n"
0102       << "  total energy of all clusters: " << totalE << "\n";
0103   if (noSuperClusters > 0) {
0104     averEnergy = totalE / noSuperClusters;
0105     str << "  average SuperCluster energy = " << averEnergy << "\n";
0106   }
0107   edm::LogInfo("Multi5x5SuperClusterProducerInfo") << str.str() << "\n";
0108 }
0109
0110 void Multi5x5SuperClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0111   if (doEndcaps_)
0112     produceSuperclustersForECALPart(evt, eeClustersToken_, endcapPutToken_);
0113
0114   if (doBarrel_)
0115     produceSuperclustersForECALPart(evt, ebClustersToken_, barrelPutToken_);
0116 }
0117
0118 void Multi5x5SuperClusterProducer::produceSuperclustersForECALPart(
0119     edm::Event& evt,
0120     const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken,
0121     const edm::EDPutTokenT<reco::SuperClusterCollection>& putToken) {
0122   // get the cluster collection out and turn it to a BasicClusterRefVector:
0123   reco::CaloClusterPtrVector clusterPtrVector_p = getClusterPtrVector(evt, clustersToken);
0124
0125   // run the brem recovery and get the SC collection
0126   reco::SuperClusterCollection superclusters_ap(bremAlgo_p->makeSuperClusters(clusterPtrVector_p));
0127
0128   // count the total energy and the number of superclusters
0129   for (auto const& sc : superclusters_ap) {
0130     totalE += sc.energy();
0131     noSuperClusters++;
0132   }
0133
0134   // put the SC collection in the event
0135   evt.emplace(putToken, std::move(superclusters_ap));
0136 }
0137
0138 reco::CaloClusterPtrVector Multi5x5SuperClusterProducer::getClusterPtrVector(
0139     edm::Event& evt, const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken) const {
0140   reco::CaloClusterPtrVector clusterPtrVector_p;
0141   edm::Handle<reco::BasicClusterCollection> bccHandle;
0142   evt.getByToken(clustersToken, bccHandle);
0143
0144   const reco::BasicClusterCollection* clusterCollection_p = bccHandle.product();
0145   clusterPtrVector_p.reserve(clusterCollection_p->size());
0146   for (unsigned int i = 0; i < clusterCollection_p->size(); i++) {
0147     clusterPtrVector_p.push_back(reco::CaloClusterPtr(bccHandle, i));
0148   }
0149   return clusterPtrVector_p;
0150 }