File indexing completed on 2023-03-17 11:21:03
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
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
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
0100 declareDynArray(float, nhits, energies);
0101 unInitDynArray(int, nhits, qst);
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;
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
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:
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
0166 seedable[idx] = false;
0167 break;
0168 }
0169 }
0170 if (seedable[idx]) {
0171 for (auto neighbour : myNeighbours) {
0172
0173
0174
0175
0176 int seedlayer = (int)maybeseed.layer();
0177 switch (seedlayer) {
0178 case PFLayer::HCAL_BARREL1:
0179 case PFLayer::HCAL_ENDCAP:
0180 case PFLayer::HF_EM:
0181
0182 case PFLayer::HF_HAD:
0183
0184 auto const& nei = (*input)[neighbour];
0185 if (maybeseed.depth() != nei.depth())
0186 continue;
0187 if (std::abs(deltaPhi(maybeseed.positionREP().phi(), nei.positionREP().phi())) > dphicut &&
0188 std::abs(maybeseed.positionREP().eta() - nei.positionREP().eta()) > detacut)
0189 continue;
0190 break;
0191 }
0192
0193 usable[neighbour] = false;
0194
0195 }
0196 }
0197 }
0198
0199 LogDebug("LocalMaximumSeedFinder") << " found " << std::count(seedable.begin(), seedable.end(), true) << " seeds";
0200 }