16 std::vector<std::pair<unsigned, double> > ordered_hits;
17 ordered_hits.reserve(input->size());
18 for (
unsigned i = 0;
i < input->size(); ++
i) {
19 std::pair<unsigned, double>
val = std::make_pair(
i, input->at(
i).energy());
20 auto pos = std::upper_bound(ordered_hits.begin(), ordered_hits.end(),
val, greaterByEnergy);
21 ordered_hits.insert(pos, val);
24 for (
const auto& idx_e : ordered_hits) {
25 if (!mask[idx_e.first])
27 const unsigned idx = idx_e.first;
29 int layer = rechit.
layer();
34 int ieta = theHcalDetId.ieta();
35 int iphi = theHcalDetId.iphi();
36 int ihpd = 0, irbx = 0;
39 ihpd = (ieta < 0 ? -iphi :
iphi);
40 irbx = (ieta < 0 ? -(iphi + 5) / 4 : (iphi + 5) / 4);
43 ihpd = (ieta < 0 ? -(iphi + 1) / 2 - 100 : (iphi + 1) / 2 + 100);
44 irbx = (ieta < 0 ? -(iphi + 5) / 4 - 20 : (iphi + 5) / 4 + 20);
51 irbx = (irbx < 0 ? -1 : 1);
54 irbx = (irbx < 0 ? -21 : 21);
59 _hpds[ihpd].push_back(idx);
60 _rbxs[irbx].push_back(idx);
64 std::unordered_map<int, std::vector<unsigned> > theHPDs;
65 std::unordered_multimap<double, unsigned> theEnergies;
66 for (
const auto& itrbx :
_rbxs) {
67 if ((
std::abs(itrbx.first) < 20 && itrbx.second.size() > 30) ||
68 (
std::abs(itrbx.first) > 20 && itrbx.second.size() > 30)) {
69 const std::vector<unsigned>&
rechits = itrbx.second;
72 int nSeeds0 = rechits.size();
73 for (
unsigned jh = 0; jh < rechits.size(); ++jh) {
79 for (
auto k : neighbours4) {
80 auto const& neighbour = (*input)[
k];
81 if (neighbour.energy() > rechit.
energy()) {
86 if (neighbour.energy() > 0.4)
94 int iphi = theHcalDetId.
iphi();
95 switch (rechit.
layer()) {
97 theHPDs[
iphi].push_back(rechits[jh]);
100 theHPDs[(iphi - 1) / 2].push_back(rechits[jh]);
105 const double rhenergy = rechit.
energy();
106 theEnergies.emplace(rhenergy, rechits[jh]);
110 for (
const auto& itHPD : theHPDs) {
111 int hpdN = itHPD.first;
112 const std::vector<unsigned>& hpdHits = itHPD.second;
113 if ((
std::abs(hpdN) < 100 && hpdHits.size() > 14) || (
std::abs(hpdN) > 100 && hpdHits.size() > 14))
119 for (
const auto& itEn : theEnergies) {
122 mask[itEn.second] =
false;
123 }
else if (nn == 5) {
124 threshold = itEn.first * 5;
125 mask[itEn.second] =
false;
127 if (itEn.first < threshold)
128 mask[itEn.second] =
false;
137 std::unordered_map<int, std::vector<unsigned> >::iterator neighbour1;
138 std::unordered_map<int, std::vector<unsigned> >::iterator neighbour2;
139 std::unordered_map<int, std::vector<unsigned> >::iterator neighbour0;
140 std::unordered_map<int, std::vector<unsigned> >::iterator neighbour3;
141 unsigned size1 = 0, size2 = 0;
142 for (
const auto& ithpd :
_hpds) {
143 const std::vector<unsigned>& rechits = ithpd.second;
145 for (
const unsigned rhidx : rechits) {
147 theEnergies.emplace(rechit.
energy(), rhidx);
150 const int thehpd = ithpd.first;
153 neighbour1 = (thehpd > 0 ? _hpds.find(72) : _hpds.find(-72));
156 neighbour2 = (thehpd > 0 ? _hpds.find(1) : _hpds.find(-1));
159 neighbour1 = (thehpd > 0 ? _hpds.find(136) : _hpds.find(-136));
162 neighbour2 = (thehpd > 0 ? _hpds.find(101) : _hpds.find(-101));
165 neighbour1 = (thehpd > 0 ? _hpds.find(thehpd - 1) : _hpds.find(thehpd + 1));
166 neighbour2 = (thehpd > 0 ? _hpds.find(thehpd + 1) : _hpds.find(thehpd - 1));
169 if (neighbour1 != _hpds.end()) {
170 const int nb1 = neighbour1->first;
173 neighbour0 = (nb1 > 0 ? _hpds.find(72) : _hpds.find(-72));
176 neighbour0 = (nb1 > 0 ? _hpds.find(136) : _hpds.find(-136));
179 neighbour0 = (nb1 > 0 ? _hpds.find(nb1 - 1) : _hpds.find(nb1 + 1));
183 neighbour0 = _hpds.end();
186 if (neighbour2 != _hpds.end()) {
187 const int nb2 = neighbour2->first;
190 neighbour3 = (nb2 > 0 ? _hpds.find(1) : _hpds.find(-1));
193 neighbour3 = (nb2 > 0 ? _hpds.find(101) : _hpds.find(-101));
196 neighbour3 = (nb2 > 0 ? _hpds.find(nb2 + 1) : _hpds.find(nb2 - 1));
200 neighbour3 = _hpds.end();
203 size1 = neighbour1 != _hpds.end() ? neighbour1->second.size() : 0;
204 size2 = neighbour2 != _hpds.end() ? neighbour2->second.size() : 0;
206 if ((
abs(neighbour1->first) > 100 && neighbour1->second.size() > 15) ||
207 (
abs(neighbour1->first) < 100 && neighbour1->second.size() > 12))
208 size1 = neighbour0 != _hpds.end() ? neighbour0->second.size() : 0;
211 if ((
abs(neighbour2->first) > 100 && neighbour2->second.size() > 15) ||
212 (
abs(neighbour2->first) < 100 && neighbour2->second.size() > 12))
213 size2 = neighbour3 != _hpds.end() ? neighbour3->second.size() : 0;
215 if ((
std::abs(ithpd.first) > 100 && ithpd.second.size() > 15) ||
216 (
std::abs(ithpd.first) < 100 && ithpd.second.size() > 12)) {
217 if ((
double)(size1 + size2) / (
float)ithpd.second.size() < 1.0) {
219 double threshold = 1.0;
220 for (
const auto& itEn : theEnergies) {
222 mask[itEn.second] =
false;
223 }
else if (nn == 5) {
224 threshold = itEn.first * 2.5;
225 mask[itEn.second] =
false;
227 if (itEn.first < threshold)
228 mask[itEn.second] =
false;
std::unordered_map< int, std::vector< unsigned > > _rbxs
unsigned detId() const
rechit detId
PFLayer::Layer layer() const
rechit layer
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Abs< T >::type abs(const T &t)
float energy() const
rechit energy
int iphi() const
get the cell iphi
std::unordered_map< int, std::vector< unsigned > > _hpds
Neighbours neighbours4() const