PFClusterProducer/plugins/Basic2DGenericTopoClusterizer.cc

0001 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
0004
0005 class Basic2DGenericTopoClusterizer : public InitialClusteringStepBase {
0006   typedef Basic2DGenericTopoClusterizer B2DGT;
0007
0008 public:
0009   Basic2DGenericTopoClusterizer(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0010       : InitialClusteringStepBase(conf, cc), _useCornerCells(conf.getParameter<bool>("useCornerCells")) {}
0011   ~Basic2DGenericTopoClusterizer() override = default;
0012   Basic2DGenericTopoClusterizer(const B2DGT&) = delete;
0013   B2DGT& operator=(const B2DGT&) = delete;
0014
0015   void buildClusters(const edm::Handle<reco::PFRecHitCollection>&,
0016                      const std::vector<bool>&,
0017                      const std::vector<bool>&,
0018                      reco::PFClusterCollection&) override;
0019
0020 private:
0021   const bool _useCornerCells;
0022   void buildTopoCluster(const edm::Handle<reco::PFRecHitCollection>&,
0023                         const std::vector<bool>&,  // masked rechits
0024                         unsigned int,              //present rechit
0025                         std::vector<bool>&,        // hit usage state
0026                         reco::PFCluster&);         // the topocluster
0027 };
0028
0029 DEFINE_EDM_PLUGIN(InitialClusteringStepFactory, Basic2DGenericTopoClusterizer, "Basic2DGenericTopoClusterizer");
0030
0031 #ifdef PFLOW_DEBUG
0032 #define LOGVERB(x) edm::LogVerbatim(x)
0033 #define LOGWARN(x) edm::LogWarning(x)
0034 #define LOGERR(x) edm::LogError(x)
0035 #define LOGDRESSED(x) edm::LogInfo(x)
0036 #else
0037 #define LOGVERB(x) LogTrace(x)
0038 #define LOGWARN(x) edm::LogWarning(x)
0039 #define LOGERR(x) edm::LogError(x)
0040 #define LOGDRESSED(x) LogDebug(x)
0041 #endif
0042
0043 void Basic2DGenericTopoClusterizer::buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
0044                                                   const std::vector<bool>& rechitMask,
0045                                                   const std::vector<bool>& seedable,
0046                                                   reco::PFClusterCollection& output) {
0047   auto const& hits = *input;
0048   std::vector<bool> used(hits.size(), false);
0049   std::vector<unsigned int> seeds;
0050
0051   // get the seeds and sort them descending in energy
0052   seeds.reserve(hits.size());
0053   for (unsigned int i = 0; i < hits.size(); ++i) {
0054     if (!rechitMask[i] || !seedable[i] || used[i])
0055       continue;
0056     seeds.emplace_back(i);
0057   }
0058   // maxHeap would be better
0059   std::sort(
0060       seeds.begin(), seeds.end(), [&](unsigned int i, unsigned int j) { return hits[i].energy() > hits[j].energy(); });
0061
0062   reco::PFCluster temp;
0063   for (auto seed : seeds) {
0064     if (!rechitMask[seed] || !seedable[seed] || used[seed])
0065       continue;
0066     temp.reset();
0067     buildTopoCluster(input, rechitMask, seed, used, temp);
0068     if (!temp.recHitFractions().empty())
0069       output.push_back(temp);
0070   }
0071 }
0072
0073 void Basic2DGenericTopoClusterizer::buildTopoCluster(const edm::Handle<reco::PFRecHitCollection>& input,
0074                                                      const std::vector<bool>& rechitMask,
0075                                                      unsigned int kcell,
0076                                                      std::vector<bool>& used,
0077                                                      reco::PFCluster& topocluster) {
0078   auto const& cell = (*input)[kcell];
0079   int cell_layer = (int)cell.layer();
0080   if (cell_layer == PFLayer::HCAL_BARREL2 && std::abs(cell.positionREP().eta()) > 0.34) {
0081     cell_layer *= 100;
0082   }
0083
0084   auto const& thresholds = _thresholds.find(cell_layer)->second;
0085   double thresholdE = 0.;
0086   double thresholdPT2 = 0.;
0087
0088   for (unsigned int j = 0; j < (std::get<1>(thresholds)).size(); ++j) {
0089     int depth = std::get<0>(thresholds)[j];
0090
0091     if ((cell_layer == PFLayer::HCAL_BARREL1 && cell.depth() == depth) ||
0092         (cell_layer == PFLayer::HCAL_ENDCAP && cell.depth() == depth) ||
0093         (cell_layer != PFLayer::HCAL_BARREL1 && cell_layer != PFLayer::HCAL_ENDCAP)) {
0094       thresholdE = std::get<1>(thresholds)[j];
0095       thresholdPT2 = std::get<2>(thresholds)[j];
0096     }
0097   }
0098
0099   if (cell.energy() < thresholdE || cell.pt2() < thresholdPT2) {
0100     LOGDRESSED("GenericTopoCluster::buildTopoCluster()")
0101         << "RecHit " << cell.detId() << " with enegy " << cell.energy() << " GeV was rejected!." << std::endl;
0102     return;
0103   }
0104
0105   auto k = kcell;
0106   used[k] = true;
0107   auto ref = makeRefhit(input, k);
0108   topocluster.addRecHitFraction(reco::PFRecHitFraction(ref, 1.0));
0109
0110   auto const& neighbours = (_useCornerCells ? cell.neighbours8() : cell.neighbours4());
0111
0112   for (auto nb : neighbours) {
0113     if (used[nb] || !rechitMask[nb]) {
0114       LOGDRESSED("GenericTopoCluster::buildTopoCluster()")
0115           << "  RecHit " << cell.detId() << "\'s"
0116           << " neighbor RecHit " << input->at(nb).detId() << " with enegy " << input->at(nb).energy()
0117           << " GeV was rejected!"
0118           << " Reasons : " << used[nb] << " (used) " << !rechitMask[nb] << " (masked)." << std::endl;
0119       continue;
0120     }
0121     buildTopoCluster(input, rechitMask, nb, used, topocluster);
0122   }
0123 }