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