PFClusterProducer/plugins/LocalMaximumSeedFinder.cc

0001 #include "CommonTools/Utils/interface/DynArray.h"
0002 #include "DataFormats/Math/interface/deltaPhi.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "RecoParticleFlow/PFClusterProducer/interface/SeedFinderBase.h"
0005 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0006 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0007
0008 #include <algorithm>
0009 #include <cfloat>
0010 #include <tuple>
0011 #include <unordered_map>
0012 #include <queue>
0013
0014 class LocalMaximumSeedFinder final : public SeedFinderBase {
0015 public:
0016   LocalMaximumSeedFinder(const edm::ParameterSet& conf);
0017   LocalMaximumSeedFinder(const LocalMaximumSeedFinder&) = delete;
0018   LocalMaximumSeedFinder& operator=(const LocalMaximumSeedFinder&) = delete;
0019
0020   void findSeeds(const edm::Handle<reco::PFRecHitCollection>& input,
0021                  const std::vector<bool>& mask,
0022                  std::vector<bool>& seedable,
0023                  const HcalPFCuts*) override;
0024
0025 private:
0026   const int _nNeighbours;
0027
0028   const std::unordered_map<std::string, int> _layerMap;
0029
0030   typedef std::tuple<std::vector<int>, std::vector<double>, std::vector<double> > I3tuple;
0031
0032   std::array<I3tuple, 35> _thresholds;
0033   static constexpr int layerOffset = 15;
0034
0035   static constexpr double detacut = 0.01;
0036   static constexpr double dphicut = 0.01;
0037 };
0038
0039 DEFINE_EDM_PLUGIN(SeedFinderFactory, LocalMaximumSeedFinder, "LocalMaximumSeedFinder");
0040
0041 namespace {
0042   const reco::PFRecHit::Neighbours _noNeighbours(nullptr, 0);
0043 }
0044
0045 LocalMaximumSeedFinder::LocalMaximumSeedFinder(const edm::ParameterSet& conf)
0046     : SeedFinderBase(conf),
0047       _nNeighbours(conf.getParameter<int>("nNeighbours")),
0048       _layerMap({{"PS2", (int)PFLayer::PS2},
0049                  {"PS1", (int)PFLayer::PS1},
0050                  {"ECAL_ENDCAP", (int)PFLayer::ECAL_ENDCAP},
0051                  {"ECAL_BARREL", (int)PFLayer::ECAL_BARREL},
0052                  {"NONE", (int)PFLayer::NONE},
0053                  {"HCAL_BARREL1", (int)PFLayer::HCAL_BARREL1},
0054                  {"HCAL_BARREL2_RING0", (int)PFLayer::HCAL_BARREL2},
0055                  // hack to deal with ring1 in HO
0056                  {"HCAL_BARREL2_RING1", 19},
0057                  {"HCAL_ENDCAP", (int)PFLayer::HCAL_ENDCAP},
0058                  {"HF_EM", (int)PFLayer::HF_EM},
0059                  {"HF_HAD", (int)PFLayer::HF_HAD}}) {
0060   const std::vector<edm::ParameterSet>& thresholds = conf.getParameterSetVector("thresholdsByDetector");
0061   for (const auto& pset : thresholds) {
0062     const std::string& det = pset.getParameter<std::string>("detector");
0063     std::vector<int> depths;
0064     std::vector<double> thresh_E;
0065     std::vector<double> thresh_pT;
0066     std::vector<double> thresh_pT2;
0067
0068     if (det == std::string("HCAL_BARREL1") || det == std::string("HCAL_ENDCAP")) {
0069       depths = pset.getParameter<std::vector<int> >("depths");
0070       thresh_E = pset.getParameter<std::vector<double> >("seedingThreshold");
0071       thresh_pT = pset.getParameter<std::vector<double> >("seedingThresholdPt");
0072       if (thresh_E.size() != depths.size() || thresh_pT.size() != depths.size()) {
0073         throw cms::Exception("InvalidGatheringThreshold") << "gatheringThresholds mismatch with the numbers of depths";
0074       }
0075     } else {
0076       depths.push_back(0);
0077       thresh_E.push_back(pset.getParameter<double>("seedingThreshold"));
0078       thresh_pT.push_back(pset.getParameter<double>("seedingThresholdPt"));
0079     }
0080
0081     for (unsigned int i = 0; i < thresh_pT.size(); ++i) {
0082       thresh_pT2.push_back(thresh_pT[i] * thresh_pT[i]);
0083     }
0084
0085     auto entry = _layerMap.find(det);
0086     if (entry == _layerMap.end()) {
0087       throw cms::Exception("InvalidDetectorLayer") << "Detector layer : " << det << " is not in the list of recognized"
0088                                                    << " detector layers!";
0089     }
0090
0091     _thresholds[entry->second + layerOffset] = std::make_tuple(depths, thresh_E, thresh_pT2);
0092   }
0093 }
0094
0095 // the starting state of seedable is all false!
0096 void LocalMaximumSeedFinder::findSeeds(const edm::Handle<reco::PFRecHitCollection>& input,
0097                                        const std::vector<bool>& mask,
0098                                        std::vector<bool>& seedable,
0099                                        const HcalPFCuts* hcalCuts) {
0100   auto nhits = input->size();
0101   initDynArray(bool, nhits, usable, true);
0102   //need to run over energy sorted rechits
0103   declareDynArray(float, nhits, energies);
0104   unInitDynArray(int, nhits, qst);  // queue storage
0105   auto cmp = [&](int i, int j) { return energies[i] < energies[j]; };
0106   std::priority_queue<int, DynArray<int>, decltype(cmp)> ordered_hits(cmp, std::move(qst));
0107
0108   for (unsigned i = 0; i < nhits; ++i) {
0109     if (!mask[i])
0110       continue;  // cannot seed masked objects
0111     auto const& maybeseed = (*input)[i];
0112     energies[i] = maybeseed.energy();
0113     int seedlayer = (int)maybeseed.layer();
0114     if (seedlayer == PFLayer::HCAL_BARREL2 && std::abs(maybeseed.positionREP().eta()) > 0.34) {
0115       seedlayer = 19;
0116     }
0117     auto const& thresholds = _thresholds[seedlayer + layerOffset];
0118
0119     double thresholdE = 0.;
0120     double thresholdPT2 = 0.;
0121
0122     for (unsigned int j = 0; j < (std::get<2>(thresholds)).size(); ++j) {
0123       int depth = std::get<0>(thresholds)[j];
0124       if ((seedlayer == PFLayer::HCAL_BARREL1 && maybeseed.depth() == depth) ||
0125           (seedlayer == PFLayer::HCAL_ENDCAP && maybeseed.depth() == depth) ||
0126           (seedlayer != PFLayer::HCAL_BARREL1 && seedlayer != PFLayer::HCAL_ENDCAP)) {
0127         thresholdE = std::get<1>(thresholds)[j];
0128         thresholdPT2 = std::get<2>(thresholds)[j];
0129       }
0130     }
0131
0132     if (hcalCuts != nullptr) {  // this means, cutsFromDB is set to True in PFClusterProducer.cc
0133       if ((seedlayer == PFLayer::HCAL_BARREL1) || (seedlayer == PFLayer::HCAL_ENDCAP)) {
0134         HcalDetId thisId = maybeseed.detId();
0135         const HcalPFCut* item = hcalCuts->getValues(thisId.rawId());
0136         thresholdE = item->seedThreshold();
0137       }
0138     }
0139
0140     if (maybeseed.energy() < thresholdE || maybeseed.pt2() < thresholdPT2)
0141       usable[i] = false;
0142     if (!usable[i])
0143       continue;
0144     ordered_hits.push(i);
0145   }
0146
0147   while (!ordered_hits.empty()) {
0148     auto idx = ordered_hits.top();
0149     ordered_hits.pop();
0150     if (!usable[idx])
0151       continue;
0152     //get the neighbours of this seed
0153     auto const& maybeseed = (*input)[idx];
0154     reco::PFRecHit::Neighbours myNeighbours;
0155     switch (_nNeighbours) {
0156       case -1:
0157         myNeighbours = maybeseed.neighbours();
0158         break;
0159       case 0:  // for HF clustering
0160         myNeighbours = _noNeighbours;
0161         break;
0162       case 4:
0163         myNeighbours = maybeseed.neighbours4();
0164         break;
0165       case 8:
0166         myNeighbours = maybeseed.neighbours8();
0167         break;
0168       default:
0169         throw cms::Exception("InvalidConfiguration") << "LocalMaximumSeedFinder only accepts nNeighbors = {-1,0,4,8}";
0170     }
0171     seedable[idx] = true;
0172     for (auto neighbour : myNeighbours) {
0173       if (!mask[neighbour])
0174         continue;
0175       if (energies[neighbour] > energies[idx]) {
0176         //        std::cout << "how this can be?" << std::endl;
0177         seedable[idx] = false;
0178         break;
0179       }
0180     }
0181     if (seedable[idx]) {
0182       for (auto neighbour : myNeighbours) {
0183         //
0184         // For HCAL,
0185         // even if channel a is a neighbor of channel b, channel b may not be a neighbor of channel a.
0186         // So, perform additional checks to ensure making a hit unusable for seeding is safe.
0187         int seedlayer = (int)maybeseed.layer();
0188         switch (seedlayer) {
0189           case PFLayer::HCAL_BARREL1:
0190           case PFLayer::HCAL_ENDCAP:
0191           case PFLayer::HF_EM:   // with the current HF setting, we won't see this case
0192                                  // but this can come in if we change _nNeighbours for HF.
0193           case PFLayer::HF_HAD:  // same as above
0194             // HO has only one depth and eta-phi segmentation is regular, so no need to make this check
0195             auto const& nei = (*input)[neighbour];
0196             if (maybeseed.depth() != nei.depth())
0197               continue;  // masking is done only if the neighbor is on the same depth layer as the seed
0198             if (std::abs(deltaPhi(maybeseed.positionREP().phi(), nei.positionREP().phi())) > dphicut &&
0199                 std::abs(maybeseed.positionREP().eta() - nei.positionREP().eta()) > detacut)
0200               continue;  // masking is done only if the neighbor is on the swiss-cross w.r.t. the seed
0201             break;
0202         }
0203
0204         usable[neighbour] = false;
0205
0206       }  // for-loop
0207     }
0208   }
0209
0210   LogDebug("LocalMaximumSeedFinder") << " found " << std::count(seedable.begin(), seedable.end(), true) << " seeds";
0211 }