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