CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
RBXAndHPDCleaner.cc
Go to the documentation of this file.
4 
5 #include <cmath>
6 #include <unordered_map>
7 
9 public:
12  RBXAndHPDCleaner(const RBXAndHPDCleaner&) = delete;
14 
15  void clean(const edm::Handle<reco::PFRecHitCollection>& input, std::vector<bool>& mask) override;
16 
17 private:
18  std::unordered_map<int, std::vector<unsigned> > _hpds, _rbxs;
19 };
20 
22 
23 namespace {
24  bool greaterByEnergy(const std::pair<unsigned, double>& a, const std::pair<unsigned, double>& b) {
25  return a.second > b.second;
26  }
27 } // namespace
28 
30  _hpds.clear();
31  _rbxs.clear();
32  //need to run over energy sorted rechits
33  std::vector<std::pair<unsigned, double> > ordered_hits;
34  ordered_hits.reserve(input->size());
35  for (unsigned i = 0; i < input->size(); ++i) {
36  std::pair<unsigned, double> val = std::make_pair(i, input->at(i).energy());
37  auto pos = std::upper_bound(ordered_hits.begin(), ordered_hits.end(), val, greaterByEnergy);
38  ordered_hits.insert(pos, val);
39  }
40 
41  for (const auto& idx_e : ordered_hits) {
42  if (!mask[idx_e.first])
43  continue;
44  const unsigned idx = idx_e.first;
45  const reco::PFRecHit& rechit = input->at(idx);
46  int layer = rechit.layer();
48  break;
49  }
50  HcalDetId theHcalDetId(rechit.detId());
51  int ieta = theHcalDetId.ieta();
52  int iphi = theHcalDetId.iphi();
53  int ihpd = 0, irbx = 0;
54  switch (layer) {
56  ihpd = (ieta < 0 ? -iphi : iphi);
57  irbx = (ieta < 0 ? -(iphi + 5) / 4 : (iphi + 5) / 4);
58  break;
60  ihpd = (ieta < 0 ? -(iphi + 1) / 2 - 100 : (iphi + 1) / 2 + 100);
61  irbx = (ieta < 0 ? -(iphi + 5) / 4 - 20 : (iphi + 5) / 4 + 20);
62  break;
63  default:
64  break;
65  }
66  switch (std::abs(irbx)) {
67  case 19:
68  irbx = (irbx < 0 ? -1 : 1);
69  break;
70  case 39:
71  irbx = (irbx < 0 ? -21 : 21);
72  break;
73  default:
74  break;
75  }
76  _hpds[ihpd].push_back(idx);
77  _rbxs[irbx].push_back(idx);
78  }
79  // loop on the rbx's we found and clean RBX's with tons of rechits
80  // and lots of energy
81  std::unordered_map<int, std::vector<unsigned> > theHPDs;
82  std::unordered_multimap<double, unsigned> theEnergies;
83  for (const auto& itrbx : _rbxs) {
84  if ((std::abs(itrbx.first) < 20 && itrbx.second.size() > 30) ||
85  (std::abs(itrbx.first) > 20 && itrbx.second.size() > 30)) {
86  const std::vector<unsigned>& rechits = itrbx.second;
87  theHPDs.clear();
88  theEnergies.clear();
89  int nSeeds0 = rechits.size();
90  for (unsigned jh = 0; jh < rechits.size(); ++jh) {
91  const reco::PFRecHit& rechit = (*input)[jh];
92  // check if rechit is a seed
93  unsigned nN = 0; // neighbours over threshold
94  bool isASeed = true;
95  auto const& neighbours4 = rechit.neighbours4();
96  for (auto k : neighbours4) {
97  auto const& neighbour = (*input)[k];
98  if (neighbour.energy() > rechit.energy()) {
99  --nSeeds0;
100  isASeed = false;
101  break;
102  } else {
103  if (neighbour.energy() > 0.4)
104  ++nN;
105  }
106  }
107  if (isASeed && !nN)
108  --nSeeds0;
109 
110  HcalDetId theHcalDetId(rechit.detId());
111  int iphi = theHcalDetId.iphi();
112  switch (rechit.layer()) {
114  theHPDs[iphi].push_back(rechits[jh]);
115  break;
117  theHPDs[(iphi - 1) / 2].push_back(rechits[jh]);
118  break;
119  default:
120  break;
121  }
122  const double rhenergy = rechit.energy();
123  theEnergies.emplace(rhenergy, rechits[jh]);
124  }
125  if (nSeeds0 > 6) {
126  unsigned nHPD15 = 0;
127  for (const auto& itHPD : theHPDs) {
128  int hpdN = itHPD.first;
129  const std::vector<unsigned>& hpdHits = itHPD.second;
130  if ((std::abs(hpdN) < 100 && hpdHits.size() > 14) || (std::abs(hpdN) > 100 && hpdHits.size() > 14))
131  ++nHPD15;
132  }
133  if (nHPD15 > 1) {
134  unsigned nn = 0;
135  double threshold = 1.0;
136  for (const auto& itEn : theEnergies) {
137  ++nn;
138  if (nn < 5) {
139  mask[itEn.second] = false;
140  } else if (nn == 5) {
141  threshold = itEn.first * 5;
142  mask[itEn.second] = false;
143  } else {
144  if (itEn.first < threshold)
145  mask[itEn.second] = false;
146  }
147  }
148  }
149  }
150  }
151  }
152 
153  // loop on the HPDs
154  std::unordered_map<int, std::vector<unsigned> >::iterator neighbour1;
155  std::unordered_map<int, std::vector<unsigned> >::iterator neighbour2;
156  std::unordered_map<int, std::vector<unsigned> >::iterator neighbour0;
157  std::unordered_map<int, std::vector<unsigned> >::iterator neighbour3;
158  unsigned size1 = 0, size2 = 0;
159  for (const auto& ithpd : _hpds) {
160  const std::vector<unsigned>& rechits = ithpd.second;
161  theEnergies.clear();
162  for (const unsigned rhidx : rechits) {
163  const reco::PFRecHit& rechit = input->at(rhidx);
164  theEnergies.emplace(rechit.energy(), rhidx);
165  }
166 
167  const int thehpd = ithpd.first;
168  switch (std::abs(thehpd)) {
169  case 1:
170  neighbour1 = (thehpd > 0 ? _hpds.find(72) : _hpds.find(-72));
171  break;
172  case 72:
173  neighbour2 = (thehpd > 0 ? _hpds.find(1) : _hpds.find(-1));
174  break;
175  case 101:
176  neighbour1 = (thehpd > 0 ? _hpds.find(136) : _hpds.find(-136));
177  break;
178  case 136:
179  neighbour2 = (thehpd > 0 ? _hpds.find(101) : _hpds.find(-101));
180  break;
181  default:
182  neighbour1 = (thehpd > 0 ? _hpds.find(thehpd - 1) : _hpds.find(thehpd + 1));
183  neighbour2 = (thehpd > 0 ? _hpds.find(thehpd + 1) : _hpds.find(thehpd - 1));
184  break;
185  }
186  if (neighbour1 != _hpds.end()) {
187  const int nb1 = neighbour1->first;
188  switch (std::abs(nb1)) {
189  case 1:
190  neighbour0 = (nb1 > 0 ? _hpds.find(72) : _hpds.find(-72));
191  break;
192  case 101:
193  neighbour0 = (nb1 > 0 ? _hpds.find(136) : _hpds.find(-136));
194  break;
195  default:
196  neighbour0 = (nb1 > 0 ? _hpds.find(nb1 - 1) : _hpds.find(nb1 + 1));
197  break;
198  }
199  } else {
200  neighbour0 = _hpds.end();
201  }
202 
203  if (neighbour2 != _hpds.end()) {
204  const int nb2 = neighbour2->first;
205  switch (std::abs(nb2)) {
206  case 72:
207  neighbour3 = (nb2 > 0 ? _hpds.find(1) : _hpds.find(-1));
208  break;
209  case 136:
210  neighbour3 = (nb2 > 0 ? _hpds.find(101) : _hpds.find(-101));
211  break;
212  default:
213  neighbour3 = (nb2 > 0 ? _hpds.find(nb2 + 1) : _hpds.find(nb2 - 1));
214  break;
215  }
216  } else {
217  neighbour3 = _hpds.end();
218  }
219 
220  size1 = neighbour1 != _hpds.end() ? neighbour1->second.size() : 0;
221  size2 = neighbour2 != _hpds.end() ? neighbour2->second.size() : 0;
222  if (size1 > 10) {
223  if ((abs(neighbour1->first) > 100 && neighbour1->second.size() > 15) ||
224  (abs(neighbour1->first) < 100 && neighbour1->second.size() > 12))
225  size1 = neighbour0 != _hpds.end() ? neighbour0->second.size() : 0;
226  }
227  if (size2 > 10) {
228  if ((abs(neighbour2->first) > 100 && neighbour2->second.size() > 15) ||
229  (abs(neighbour2->first) < 100 && neighbour2->second.size() > 12))
230  size2 = neighbour3 != _hpds.end() ? neighbour3->second.size() : 0;
231  }
232  if ((std::abs(ithpd.first) > 100 && ithpd.second.size() > 15) ||
233  (std::abs(ithpd.first) < 100 && ithpd.second.size() > 12)) {
234  if ((double)(size1 + size2) / (float)ithpd.second.size() < 1.0) {
235  unsigned nn = 0;
236  double threshold = 1.0;
237  for (const auto& itEn : theEnergies) {
238  if (nn < 5) {
239  mask[itEn.second] = false;
240  } else if (nn == 5) {
241  threshold = itEn.first * 2.5;
242  mask[itEn.second] = false;
243  } else {
244  if (itEn.first < threshold)
245  mask[itEn.second] = false;
246  }
247  }
248  }
249  }
250  }
251 }
pfDeepBoostedJetPreprocessParams_cfi.upper_bound
upper_bound
Definition: pfDeepBoostedJetPreprocessParams_cfi.py:16
mps_fire.i
i
Definition: mps_fire.py:428
input
static const std::string input
Definition: EdmProvDump.cc:48
reco::PFRecHit::energy
float energy() const
rechit energy
Definition: PFRecHit.h:99
pos
Definition: PixelAliasList.h:18
PFLayer::HCAL_ENDCAP
Definition: PFLayer.h:37
edm::Handle< reco::PFRecHitCollection >
RBXAndHPDCleaner::clean
void clean(const edm::Handle< reco::PFRecHitCollection > &input, std::vector< bool > &mask) override
Definition: RBXAndHPDCleaner.cc:29
RBXAndHPDCleaner::operator=
RBXAndHPDCleaner & operator=(const RBXAndHPDCleaner &)=delete
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
LEDCalibrationChannels.iphi
iphi
Definition: LEDCalibrationChannels.py:64
PFLayer::HCAL_BARREL1
Definition: PFLayer.h:35
HI_PhotonSkim_cff.rechits
rechits
Definition: HI_PhotonSkim_cff.py:76
dqmdumpme.k
k
Definition: dqmdumpme.py:60
DEFINE_EDM_PLUGIN
#define DEFINE_EDM_PLUGIN(factory, type, name)
Definition: PluginFactory.h:124
b
double b
Definition: hdecay.h:118
phase1PixelTopology::layer
constexpr std::array< uint8_t, layerIndexSize > layer
Definition: phase1PixelTopology.h:99
RBXAndHPDCleaner::_hpds
std::unordered_map< int, std::vector< unsigned > > _hpds
Definition: RBXAndHPDCleaner.cc:18
RecHitTopologicalCleanerBase.h
LEDCalibrationChannels.ieta
ieta
Definition: LEDCalibrationChannels.py:63
edm::ParameterSet
Definition: ParameterSet.h:47
a
double a
Definition: hdecay.h:119
edmplugin::PluginFactory
Definition: PluginFactory.h:34
HcalDetId.h
HcalDetId
Definition: HcalDetId.h:12
cc
reco::PFRecHit::neighbours4
Neighbours neighbours4() const
Definition: PFRecHit.h:81
groupFilesInBlocks.nn
nn
Definition: groupFilesInBlocks.py:150
reco::PFRecHit::detId
unsigned detId() const
rechit detId
Definition: PFRecHit.h:93
heppy_batch.val
val
Definition: heppy_batch.py:351
RBXAndHPDCleaner::_rbxs
std::unordered_map< int, std::vector< unsigned > > _rbxs
Definition: RBXAndHPDCleaner.cc:18
reco::PFRecHit
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
reco::PFRecHit::layer
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:96
RecHitTopologicalCleanerBase
Definition: RecHitTopologicalCleanerBase.h:13
ConsumesCollector.h
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
remoteMonitoring_LED_IterMethod_cfg.threshold
threshold
Definition: remoteMonitoring_LED_IterMethod_cfg.py:430
edm::ConsumesCollector
Definition: ConsumesCollector.h:45
RBXAndHPDCleaner::RBXAndHPDCleaner
RBXAndHPDCleaner(const edm::ParameterSet &conf, edm::ConsumesCollector &cc)
Definition: RBXAndHPDCleaner.cc:10
RBXAndHPDCleaner
Definition: RBXAndHPDCleaner.cc:8