CMS 3D CMS Logo
MatrixHCALDetIds.cc
Go to the documentation of this file.
1 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
2 
3 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
4 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
5 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
6 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
7 #include "FWCore/MessageLogger/interface/MessageLogger.h"
8 
9 #include <algorithm>
10 
11 namespace spr {
12 
13  std::vector<DetId> matrixHCALIds(
14  std::vector<DetId>& dets, const HcalTopology* topology, int ieta, int iphi, bool includeHO, bool debug) {
15  if (debug) {
16  edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Add " << ieta << " rows and " << iphi << " columns of cells for "
17  << dets.size() << " cells";
18  spr::debugHcalDets(0, dets);
19  }
20 
21  std::vector<DetId> vdetN = spr::newHCALIdNS(dets, 0, topology, true, ieta, iphi, debug);
22  std::vector<DetId> vdetS = spr::newHCALIdNS(dets, 0, topology, false, ieta, iphi, debug);
23  for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
24  if (std::count(vdetN.begin(), vdetN.end(), vdetS[i1]) == 0)
25  vdetN.push_back(vdetS[i1]);
26  }
27 
28  vdetS = spr::matrixHCALIdsDepth(vdetN, topology, includeHO, debug);
29 
30  if (debug) {
31  edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdetS.size();
32  spr::debugHcalDets(0, vdetS);
33  }
34  return vdetS;
35  }
36 
37  std::vector<DetId> matrixHCALIds(const DetId& det,
38  const CaloGeometry* geo,
39  const HcalTopology* topology,
40  double dR,
41  const GlobalVector& trackMom,
42  bool includeHO,
43  bool debug) {
44  HcalDetId hcdet = HcalDetId(det);
45  GlobalPoint core = (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(hcdet)))->getPosition(hcdet);
46  std::vector<DetId> dets, vdetx;
47  dets.push_back(det);
48  int ietaphi = (int)(dR / 15.0) + 1;
49  std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ietaphi, ietaphi, includeHO, debug);
50  for (unsigned int i = 0; i < vdets.size(); ++i) {
51  HcalDetId hcdet = HcalDetId(vdets[i]);
52  GlobalPoint rpoint = (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(hcdet)))->getPosition(hcdet);
53  if (spr::getDistInPlaneTrackDir(core, trackMom, rpoint) < dR) {
54  vdetx.push_back(vdets[i]);
55  }
56  }
57 
58  if (debug) {
59  edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Final List of cells for dR " << dR << " is with " << vdetx.size()
60  << " from original list of " << vdets.size() << " cells";
61  spr::debugHcalDets(0, vdetx);
62  }
63  return vdetx;
64  }
65 
66  std::vector<DetId> matrixHCALIds(std::vector<DetId>& dets,
67  const HcalTopology* topology,
68  int ietaE,
69  int ietaW,
70  int iphiN,
71  int iphiS,
72  bool includeHO,
73  bool debug) {
74  if (debug) {
75  edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Add " << ietaE << "|" << ietaW << " rows and " << iphiN << "|"
76  << iphiS << " columns of cells for " << dets.size() << " cells";
77  spr::debugHcalDets(0, dets);
78  }
79 
80  std::vector<DetId> vdetN = spr::newHCALIdNS(dets, 0, topology, true, ietaE, ietaW, iphiN, iphiS, debug);
81  std::vector<DetId> vdetS = spr::newHCALIdNS(dets, 0, topology, false, ietaE, ietaW, iphiN, iphiS, debug);
82  for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
83  if (std::count(vdetN.begin(), vdetN.end(), vdetS[i1]) == 0)
84  vdetN.push_back(vdetS[i1]);
85  }
86 
87  vdetS = spr::matrixHCALIdsDepth(vdetN, topology, includeHO, debug);
88 
89  if (debug) {
90  edm::LogVerbatim("IsoTrack") << "matrixHCALIds::Total number of cells found is " << vdetS.size();
91  spr::debugHcalDets(0, vdetS);
92  }
93  return vdetS;
94  }
95 
96  std::vector<DetId> newHCALIdNS(std::vector<DetId>& dets,
97  unsigned int last,
98  const HcalTopology* topology,
99  bool shiftNorth,
100  int ieta,
101  int iphi,
102  bool debug) {
103  if (debug) {
104  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Add " << iphi << " columns of cells along " << shiftNorth << " for "
105  << (dets.size() - last) << " cells";
106  spr::debugHcalDets(last, dets);
107  }
108 
109  std::vector<DetId> vdets;
110  vdets.insert(vdets.end(), dets.begin(), dets.end());
111  std::vector<DetId> vdetE, vdetW;
112  if (last == 0) {
113  vdetE = spr::newHCALIdEW(dets, last, topology, true, ieta, debug);
114  vdetW = spr::newHCALIdEW(dets, last, topology, false, ieta, debug);
115  for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
116  if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0)
117  vdets.push_back(vdetW[i1]);
118  }
119  for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
120  if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0)
121  vdets.push_back(vdetE[i1]);
122  }
123 
124  if (debug) {
125  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::With Added cells along E/W results a set of "
126  << (vdets.size() - dets.size()) << " new cells";
127  spr::debugHcalDets(dets.size(), vdets);
128  }
129  }
130  unsigned int last0 = vdets.size();
131  if (iphi > 0) {
132  std::vector<DetId> vdetnew;
133  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
134  std::vector<DetId> vdet;
135  if (shiftNorth)
136  vdet = topology->north(dets[i1]);
137  else
138  vdet = topology->south(dets[i1]);
139  for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
140  if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
141  vdetnew.push_back(vdet[i2]);
142  }
143  }
144  iphi--;
145  vdetE = spr::newHCALIdEW(vdetnew, 0, topology, true, ieta, debug);
146  vdetW = spr::newHCALIdEW(vdetnew, 0, topology, false, ieta, debug);
147  for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
148  if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
149  std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0)
150  vdets.push_back(vdetW[i2]);
151  }
152  for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
153  if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
154  std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0)
155  vdets.push_back(vdetE[i2]);
156  }
157  last = vdets.size();
158  vdets.insert(vdets.end(), vdetnew.begin(), vdetnew.end());
159 
160  if (debug) {
161  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Addition results a set of " << (vdets.size() - last0)
162  << " new cells";
163  spr::debugHcalDets(last0, vdets);
164  }
165  last0 = last;
166  }
167 
168  if (iphi > 0) {
169  last = last0;
170  return spr::newHCALIdNS(vdets, last, topology, shiftNorth, ieta, iphi, debug);
171  } else {
172  if (debug) {
173  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Final list consists of " << vdets.size() << " cells";
174  spr::debugHcalDets(0, vdets);
175  }
176  return vdets;
177  }
178  }
179 
180  std::vector<DetId> newHCALIdNS(std::vector<DetId>& dets,
181  unsigned int last,
182  const HcalTopology* topology,
183  bool shiftNorth,
184  int ietaE,
185  int ietaW,
186  int iphiN,
187  int iphiS,
188  bool debug) {
189  if (debug) {
190  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Add " << iphiN << "|" << iphiS << " columns of cells along "
191  << shiftNorth << " for " << (dets.size() - last) << " cells";
192  spr::debugHcalDets(last, dets);
193  }
194 
195  std::vector<DetId> vdets;
196  vdets.insert(vdets.end(), dets.begin(), dets.end());
197  std::vector<DetId> vdetE, vdetW;
198  if (last == 0) {
199  vdetE = spr::newHCALIdEW(dets, last, topology, true, ietaE, ietaW, debug);
200  vdetW = spr::newHCALIdEW(dets, last, topology, false, ietaE, ietaW, debug);
201  for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
202  if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0)
203  vdets.push_back(vdetW[i1]);
204  }
205  for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
206  if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0)
207  vdets.push_back(vdetE[i1]);
208  }
209 
210  if (debug) {
211  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::With Added cells along E/W results a set of "
212  << (vdets.size() - dets.size()) << " new cells";
213  spr::debugHcalDets(dets.size(), vdets);
214  }
215  }
216  unsigned int last0 = vdets.size();
217  int iphi = iphiS;
218  if (shiftNorth)
219  iphi = iphiN;
220  if (iphi > 0) {
221  std::vector<DetId> vdetnew;
222  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
223  std::vector<DetId> vdet;
224  if (shiftNorth)
225  vdet = topology->north(dets[i1]);
226  else
227  vdet = topology->south(dets[i1]);
228  for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
229  if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
230  vdetnew.push_back(vdet[i2]);
231  }
232  }
233  iphi--;
234  vdetE = spr::newHCALIdEW(vdetnew, 0, topology, true, ietaE, ietaW, debug);
235  vdetW = spr::newHCALIdEW(vdetnew, 0, topology, false, ietaE, ietaW, debug);
236  for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
237  if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
238  std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0)
239  vdets.push_back(vdetW[i2]);
240  }
241  for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
242  if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
243  std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0)
244  vdets.push_back(vdetE[i2]);
245  }
246  last = vdets.size();
247  vdets.insert(vdets.end(), vdetnew.begin(), vdetnew.end());
248 
249  if (debug) {
250  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Addition results a set of " << (vdets.size() - last0)
251  << " new cells";
252  spr::debugHcalDets(last0, vdets);
253  }
254  last0 = last;
255  }
256  if (shiftNorth)
257  iphiN = iphi;
258  else
259  iphiS = iphi;
260 
261  if (iphi > 0) {
262  last = last0;
263  return spr::newHCALIdNS(vdets, last, topology, shiftNorth, ietaE, ietaW, iphiN, iphiS, debug);
264  } else {
265  if (debug) {
266  edm::LogVerbatim("IsoTrack") << "newHCALIdNS::Final list consists of " << vdets.size() << " cells";
267  spr::debugHcalDets(0, vdets);
268  }
269  return vdets;
270  }
271  }
272 
273  std::vector<DetId> newHCALIdEW(
274  std::vector<DetId>& dets, unsigned int last, const HcalTopology* topology, bool shiftEast, int ieta, bool debug) {
275  if (debug) {
276  edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Add " << ieta << " rows of cells along " << shiftEast << " for "
277  << (dets.size() - last) << " cells";
278  spr::debugHcalDets(last, dets);
279  }
280 
281  std::vector<DetId> vdets;
282  vdets.insert(vdets.end(), dets.begin(), dets.end());
283  if (ieta > 0) {
284  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
285  std::vector<DetId> vdet;
286  if (shiftEast)
287  vdet = topology->east(dets[i1]);
288  else
289  vdet = topology->west(dets[i1]);
290  for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
291  if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
292  vdets.push_back(vdet[i2]);
293  }
294  }
295  ieta--;
296  }
297 
298  if (debug) {
299  edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Addition results a set of " << (vdets.size() - dets.size())
300  << " new cells";
301  spr::debugHcalDets(dets.size(), vdets);
302  }
303  if (ieta > 0) {
304  last = dets.size();
305  return spr::newHCALIdEW(vdets, last, topology, shiftEast, ieta, debug);
306  } else {
307  if (debug) {
308  edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
309  spr::debugHcalDets(0, vdets);
310  }
311  return vdets;
312  }
313  }
314 
315  std::vector<DetId> newHCALIdEW(std::vector<DetId>& dets,
316  unsigned int last,
317  const HcalTopology* topology,
318  bool shiftEast,
319  int ietaE,
320  int ietaW,
321  bool debug) {
322  if (debug) {
323  edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Add " << ietaE << "|" << ietaW << " rows of cells along "
324  << shiftEast << " for " << (dets.size() - last) << " cells";
325  spr::debugHcalDets(last, dets);
326  }
327 
328  int ieta = ietaW;
329  if (shiftEast)
330  ieta = ietaE;
331  std::vector<DetId> vdets;
332  vdets.insert(vdets.end(), dets.begin(), dets.end());
333  if (ieta > 0) {
334  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
335  std::vector<DetId> vdet;
336  if (shiftEast)
337  vdet = topology->east(dets[i1]);
338  else
339  vdet = topology->west(dets[i1]);
340  for (unsigned int i2 = 0; i2 < vdet.size(); i2++) {
341  if (std::count(vdets.begin(), vdets.end(), vdet[i2]) == 0)
342  vdets.push_back(vdet[i2]);
343  }
344  }
345  ieta--;
346  }
347  if (shiftEast)
348  ietaE = ieta;
349  else
350  ietaW = ieta;
351 
352  if (debug) {
353  edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Addition results a set of " << (vdets.size() - dets.size())
354  << " new cells";
355  spr::debugHcalDets(dets.size(), vdets);
356  }
357 
358  if (ieta > 0) {
359  last = dets.size();
360  return spr::newHCALIdEW(vdets, last, topology, shiftEast, ietaE, ietaW, debug);
361  } else {
362  if (debug) {
363  edm::LogVerbatim("IsoTrack") << "newHCALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
364  spr::debugHcalDets(0, vdets);
365  }
366  return vdets;
367  }
368  }
369 
370  std::vector<DetId> matrixHCALIdsDepth(std::vector<DetId>& dets,
371  const HcalTopology* topology,
372  bool includeHO,
373  bool debug) {
374  if (debug) {
375  edm::LogVerbatim("IsoTrack") << "matrixHCALIdsDepth::Add cells with higher depths with HO"
376  << "Flag set to " << includeHO << " to existing " << dets.size() << " cells";
377  spr::debugHcalDets(0, dets);
378  }
379  std::vector<DetId> vdets(dets);
380  for (unsigned int i1 = 0; i1 < dets.size(); i1++) {
381  HcalDetId vdet = dets[i1];
382  for (int idepth = 0; idepth < 3; idepth++) {
383  std::vector<DetId> vUpDetId = topology->up(vdet);
384  if (!vUpDetId.empty()) {
385  if (includeHO || vUpDetId[0].subdetId() != (int)(HcalOuter)) {
386  int n = std::count(vdets.begin(), vdets.end(), vUpDetId[0]);
387  if (n == 0) {
388  if (debug)
389  edm::LogVerbatim("IsoTrack")
390  << "matrixHCALIdsDepth:: Depth " << idepth << " " << vdet << " " << (HcalDetId)vUpDetId[0];
391  vdets.push_back(vUpDetId[0]);
392  }
393  }
394  vdet = vUpDetId[0];
395  }
396  }
397  }
398 
399  if (debug) {
400  edm::LogVerbatim("IsoTrack") << "matrixHCALIdsDepth::Final list contains " << vdets.size() << " cells";
401  spr::debugHcalDets(0, vdets);
402  }
403  return vdets;
404  }
405 
406 } // namespace spr
spr::newHCALIdNS
std::vector< DetId > newHCALIdNS(std::vector< DetId > &dets, unsigned int last, const HcalTopology *topology, bool shiftNorth, int ieta, int iphi, bool debug=false)
Definition: MatrixHCALDetIds.cc:96
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
spr::getDistInPlaneTrackDir
double getDistInPlaneTrackDir(const GlobalPoint &caloPoint, const GlobalVector &caloVector, const GlobalPoint &rechitPoint, bool debug=false)
Definition: FindDistCone.cc:9
mps_fire.i
i
Definition: mps_fire.py:428
HcalTopology::south
std::vector< DetId > south(const DetId &id) const override
Definition: HcalTopology.cc:436
core
Definition: __init__.py:1
testProducerWithPsetDescEmpty_cfi.i2
i2
Definition: testProducerWithPsetDescEmpty_cfi.py:46
spr::debugHcalDets
void debugHcalDets(unsigned int, std::vector< DetId > &)
Definition: DebugInfo.cc:39
HcalTopology::up
std::vector< DetId > up(const DetId &id) const override
Definition: HcalTopology.cc:447
HcalTopology::east
std::vector< DetId > east(const DetId &id) const override
Definition: HcalTopology.cc:403
DetId
Definition: DetId.h:17
debug
#define debug
Definition: HDRShower.cc:19
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
testProducerWithPsetDescEmpty_cfi.i1
i1
Definition: testProducerWithPsetDescEmpty_cfi.py:45
HcalTopology::north
std::vector< DetId > north(const DetId &id) const override
Definition: HcalTopology.cc:425
spr::matrixHCALIdsDepth
std::vector< DetId > matrixHCALIdsDepth(std::vector< DetId > &dets, const HcalTopology *topology, bool includeHO=false, bool debug=false)
Definition: MatrixHCALDetIds.cc:370
spr::matrixHCALIds
std::vector< DetId > matrixHCALIds(std::vector< DetId > &dets, const HcalTopology *topology, int ieta, int iphi, bool includeHO=false, bool debug=false)
Definition: MatrixHCALDetIds.cc:13
HcalTopology::west
std::vector< DetId > west(const DetId &id) const override
Definition: HcalTopology.cc:414
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:34
spr::newHCALIdEW
std::vector< DetId > newHCALIdEW(std::vector< DetId > &dets, unsigned int last, const HcalTopology *topology, bool shiftEast, int ieta, bool debug=false)
Definition: MatrixHCALDetIds.cc:273
Generated for CMSSW Reference Manual by doxygen 1.8.14