HiEgammaAlgos/plugins/HiSuperClusterProducer.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/PluginManager/interface/ModuleDef.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 #include "RecoHI/HiEgammaAlgos/interface/HiBremRecoveryClusterAlgo.h"
0014
0015 #include <iostream>
0016 #include <memory>
0017 #include <sstream>
0018 #include <vector>
0019
0020 class HiSuperClusterProducer : public edm::stream::EDProducer<> {
0021 public:
0022   HiSuperClusterProducer(const edm::ParameterSet& ps);
0023
0024   void produce(edm::Event&, const edm::EventSetup&) override;
0025   virtual void endJob();
0026
0027 private:
0028   const edm::EDPutTokenT<reco::SuperClusterCollection> endcapSuperclusterPut_;
0029   const edm::EDPutTokenT<reco::SuperClusterCollection> barrelSuperclusterPut_;
0030
0031   const edm::EDGetTokenT<reco::BasicClusterCollection> eeClustersToken_;
0032   const edm::EDGetTokenT<reco::BasicClusterCollection> ebClustersToken_;
0033
0034   const float barrelEtaSearchRoad_;
0035   const float barrelPhiSearchRoad_;
0036   const float endcapEtaSearchRoad_;
0037   const float endcapPhiSearchRoad_;
0038   const float seedTransverseEnergyThreshold_;
0039   const float barrelBCEnergyThreshold_;
0040   const float endcapBCEnergyThreshold_;
0041
0042   const bool doBarrel_;
0043   const bool doEndcaps_;
0044
0045   std::unique_ptr<HiBremRecoveryClusterAlgo> bremAlgo_p;
0046
0047   double totalE = 0.0;
0048   int noSuperClusters = 0;
0049
0050   void getClusterPtrVector(edm::Event& evt,
0051                            const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken,
0052                            reco::CaloClusterPtrVector*);
0053
0054   void produceSuperclustersForECALPart(edm::Event& evt,
0055                                        const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken,
0056                                        edm::EDPutTokenT<reco::SuperClusterCollection> const& putToken);
0057
0058   void outputValidationInfo(reco::SuperClusterCollection& superclusterCollection);
0059 };
0060
0061 //define this as a plug-in
0062 #include "FWCore/Framework/interface/MakerMacros.h"
0063 DEFINE_FWK_MODULE(HiSuperClusterProducer);
0064
0065 HiSuperClusterProducer::HiSuperClusterProducer(const edm::ParameterSet& ps)
0066     : endcapSuperclusterPut_{produces<reco::SuperClusterCollection>(
0067           ps.getParameter<std::string>("endcapSuperclusterCollection"))},
0068       barrelSuperclusterPut_{
0069           produces<reco::SuperClusterCollection>(ps.getParameter<std::string>("barrelSuperclusterCollection"))},
0070       eeClustersToken_{consumes<reco::BasicClusterCollection>(
0071           edm::InputTag(ps.getParameter<std::string>("endcapClusterProducer"),
0072                         ps.getParameter<std::string>("endcapClusterCollection")))},
0073       ebClustersToken_{consumes<reco::BasicClusterCollection>(
0074           edm::InputTag(ps.getParameter<std::string>("barrelClusterProducer"),
0075                         ps.getParameter<std::string>("barrelClusterCollection")))},
0076       barrelEtaSearchRoad_{float(ps.getParameter<double>("barrelEtaSearchRoad"))},
0077       barrelPhiSearchRoad_{float(ps.getParameter<double>("barrelPhiSearchRoad"))},
0078       endcapEtaSearchRoad_{float(ps.getParameter<double>("endcapEtaSearchRoad"))},
0079       endcapPhiSearchRoad_{float(ps.getParameter<double>("endcapPhiSearchRoad"))},
0080       seedTransverseEnergyThreshold_{float(ps.getParameter<double>("seedTransverseEnergyThreshold"))},
0081       barrelBCEnergyThreshold_{float(ps.getParameter<double>("barrelBCEnergyThreshold"))},
0082       endcapBCEnergyThreshold_{float(ps.getParameter<double>("endcapBCEnergyThreshold"))},
0083       doBarrel_{ps.getParameter<bool>("doBarrel")},
0084       doEndcaps_{ps.getParameter<bool>("doEndcaps")}
0085
0086 {
0087   // The verbosity level
0088   HiBremRecoveryClusterAlgo::VerbosityLevel verbosity;
0089   std::string verbosityString = ps.getParameter<std::string>("VerbosityLevel");
0090   if (verbosityString == "DEBUG")
0091     verbosity = HiBremRecoveryClusterAlgo::pDEBUG;
0092   else if (verbosityString == "WARNING")
0093     verbosity = HiBremRecoveryClusterAlgo::pWARNING;
0094   else if (verbosityString == "INFO")
0095     verbosity = HiBremRecoveryClusterAlgo::pINFO;
0096   else
0097     verbosity = HiBremRecoveryClusterAlgo::pERROR;
0098
0099   if (verbosityString == "INFO") {
0100     std::cout << "Barrel BC Energy threshold = " << barrelBCEnergyThreshold_ << std::endl;
0101     std::cout << "Endcap BC Energy threshold = " << endcapBCEnergyThreshold_ << std::endl;
0102   }
0103
0104   bremAlgo_p = std::make_unique<HiBremRecoveryClusterAlgo>(barrelEtaSearchRoad_,
0105                                                            barrelPhiSearchRoad_,
0106                                                            endcapEtaSearchRoad_,
0107                                                            endcapPhiSearchRoad_,
0108                                                            seedTransverseEnergyThreshold_,
0109                                                            barrelBCEnergyThreshold_,
0110                                                            endcapBCEnergyThreshold_,
0111                                                            verbosity);
0112 }
0113
0114 void HiSuperClusterProducer::endJob() {
0115   double averEnergy = 0.;
0116   std::ostringstream str;
0117   str << "HiSuperClusterProducer::endJob()\n"
0118       << "  total # reconstructed super clusters: " << noSuperClusters << "\n"
0119       << "  total energy of all clusters: " << totalE << "\n";
0120   if (noSuperClusters > 0) {
0121     averEnergy = totalE / noSuperClusters;
0122     str << "  average SuperCluster energy = " << averEnergy << "\n";
0123   }
0124   edm::LogInfo("HiSuperClusterProducerInfo") << str.str() << "\n";
0125 }
0126
0127 void HiSuperClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0128   if (doEndcaps_)
0129     produceSuperclustersForECALPart(evt, eeClustersToken_, endcapSuperclusterPut_);
0130
0131   if (doBarrel_)
0132     produceSuperclustersForECALPart(evt, ebClustersToken_, barrelSuperclusterPut_);
0133 }
0134
0135 void HiSuperClusterProducer::produceSuperclustersForECALPart(
0136     edm::Event& evt,
0137     const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken,
0138     edm::EDPutTokenT<reco::SuperClusterCollection> const& putToken) {
0139   // get the cluster collection out and turn it to a BasicClusterRefVector:
0140   reco::CaloClusterPtrVector clusterPtrVector_p{};
0141   getClusterPtrVector(evt, clustersToken, &clusterPtrVector_p);
0142
0143   // run the brem recovery and get the SC collection
0144   reco::SuperClusterCollection superclusters_ap = bremAlgo_p->makeSuperClusters(clusterPtrVector_p);
0145
0146   // count the total energy and the number of superclusters
0147   for (auto const& cluster : superclusters_ap) {
0148     totalE += cluster.energy();
0149     noSuperClusters++;
0150   }
0151
0152   // put the SC collection in the event
0153   evt.emplace(putToken, std::move(superclusters_ap));
0154 }
0155
0156 void HiSuperClusterProducer::getClusterPtrVector(edm::Event& evt,
0157                                                  const edm::EDGetTokenT<reco::BasicClusterCollection>& clustersToken,
0158                                                  reco::CaloClusterPtrVector* clusterPtrVector_p) {
0159   edm::Handle<reco::BasicClusterCollection> bccHandle = evt.getHandle(clustersToken);
0160
0161   if (!(bccHandle.isValid())) {
0162     edm::LogError("HiSuperClusterProducerError") << "could not get a handle on the BasicCluster Collection!";
0163     clusterPtrVector_p = nullptr;
0164   }
0165
0166   const reco::BasicClusterCollection* clusterCollection_p = bccHandle.product();
0167   for (unsigned int i = 0; i < clusterCollection_p->size(); i++) {
0168     clusterPtrVector_p->push_back(reco::CaloClusterPtr(bccHandle, i));
0169   }
0170 }