CMS 3D CMS Logo
MatrixECALDetIds.cc
Go to the documentation of this file.
1 #include "DataFormats/EcalDetId/interface/EBDetId.h"
2 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
3 #include "DataFormats/EcalDetId/interface/EEDetId.h"
4 
5 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
6 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
7 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
8 
9 #include "FWCore/MessageLogger/interface/MessageLogger.h"
10 
11 #include <algorithm>
12 #include <iostream>
13 
14 namespace spr {
15 
16  void matrixECALIds(const DetId& det,
17  int ieta,
18  int iphi,
19  const CaloGeometry* geo,
20  const CaloTopology* caloTopology,
21  std::vector<DetId>& vdets,
22  bool debug,
23  bool ignoreTransition) {
24  const CaloSubdetectorTopology* barrelTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel));
25  const CaloSubdetectorTopology* endcapTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap));
26  const EcalBarrelGeometry* barrelGeom =
27  (dynamic_cast<const EcalBarrelGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel)));
28  const EcalEndcapGeometry* endcapGeom =
29  (dynamic_cast<const EcalEndcapGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap)));
30 
31  if (debug) {
32  edm::LogVerbatim("IsoTrack") << "matrixECALIds::Add " << ieta << " rows and " << iphi
33  << " columns of cells for 1 cell" << spr::debugEcalDet(0, det).str();
34  }
35  std::vector<DetId> dets(1, det);
36  std::vector<CaloDirection> dirs(1, NORTH);
37  vdets = spr::newECALIdNS(
38  dets, 0, ieta, iphi, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
39  dirs[0] = SOUTH;
40  std::vector<DetId> vdetS = spr::newECALIdNS(
41  dets, 0, ieta, iphi, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
42  for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
43  if (std::count(vdets.begin(), vdets.end(), vdetS[i1]) == 0)
44  vdets.push_back(vdetS[i1]);
45  }
46  unsigned int ndet = (2 * ieta + 1) * (2 * iphi + 1);
47  if (vdets.size() != ndet) {
48  std::vector<DetId> vdetExtra;
49  spr::extraIds(det, vdets, ieta, ieta, iphi, iphi, barrelGeom, endcapGeom, vdetExtra, debug);
50  if (!vdetExtra.empty())
51  vdets.insert(vdets.end(), vdetExtra.begin(), vdetExtra.end());
52  }
53 
54  if (debug) {
55  edm::LogVerbatim("IsoTrack") << "matrixECALIds::Total number of cells found is " << vdets.size();
56  spr::debugEcalDets(0, vdets);
57  }
58  }
59 
60  std::vector<DetId> matrixECALIds(const DetId& det,
61  int ieta,
62  int iphi,
63  const CaloGeometry* geo,
64  const CaloTopology* caloTopology,
65  bool debug,
66  bool ignoreTransition) {
67  std::vector<DetId> vdets;
68  spr::matrixECALIds(det, ieta, iphi, geo, caloTopology, vdets, debug, ignoreTransition);
69  return vdets;
70  }
71 
72  std::vector<DetId> matrixECALIds(const DetId& det,
73  double dR,
74  const GlobalVector& trackMom,
75  const CaloGeometry* geo,
76  const CaloTopology* caloTopology,
77  bool debug,
78  bool ignoreTransition) {
79  GlobalPoint core;
80  if (det.subdetId() == EcalEndcap) {
81  EEDetId EEid = EEDetId(det);
82  core = geo->getPosition(EEid);
83  } else {
84  EBDetId EBid = EBDetId(det);
85  core = geo->getPosition(EBid);
86  }
87  int ietaphi = (int)(dR / 2.0) + 1;
88  std::vector<DetId> vdets, vdetx;
89  spr::matrixECALIds(det, ietaphi, ietaphi, geo, caloTopology, vdets, debug, ignoreTransition);
90  for (unsigned int i = 0; i < vdets.size(); ++i) {
91  GlobalPoint rpoint;
92  if (vdets[i].subdetId() == EcalEndcap) {
93  EEDetId EEid = EEDetId(vdets[i]);
94  rpoint = geo->getPosition(EEid);
95  } else {
96  EBDetId EBid = EBDetId(vdets[i]);
97  rpoint = geo->getPosition(EBid);
98  }
99  if (spr::getDistInPlaneTrackDir(core, trackMom, rpoint) < dR) {
100  vdetx.push_back(vdets[i]);
101  }
102  }
103 
104  if (debug) {
105  edm::LogVerbatim("IsoTrack") << "matrixECALIds::Final List of cells for dR " << dR << " is with " << vdetx.size()
106  << " from original list of " << vdets.size();
107  spr::debugEcalDets(0, vdetx);
108  }
109  return vdetx;
110  }
111 
112  void matrixECALIds(const DetId& det,
113  int ietaE,
114  int ietaW,
115  int iphiN,
116  int iphiS,
117  const CaloGeometry* geo,
118  const CaloTopology* caloTopology,
119  std::vector<DetId>& vdets,
120  bool debug,
121  bool ignoreTransition) {
122  const CaloSubdetectorTopology* barrelTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel));
123  const CaloSubdetectorTopology* endcapTopo = (caloTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap));
124  const EcalBarrelGeometry* barrelGeom =
125  (dynamic_cast<const EcalBarrelGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel)));
126  const EcalEndcapGeometry* endcapGeom =
127  (dynamic_cast<const EcalEndcapGeometry*>(geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap)));
128 
129  if (debug) {
130  edm::LogVerbatim("IsoTrack") << "matrixECALIds::Add " << ietaE << "|" << ietaW << " rows and " << iphiN << "|"
131  << iphiS << " columns of cells for 1 cell" << spr::debugEcalDet(0, det).str();
132  }
133  std::vector<DetId> dets(1, det);
134  std::vector<CaloDirection> dirs(1, NORTH);
135  std::vector<int> jetaE(1, ietaE), jetaW(1, ietaW);
136  std::vector<int> jphiN(1, iphiN), jphiS(1, iphiS);
137  vdets = spr::newECALIdNS(dets,
138  0,
139  jetaE,
140  jetaW,
141  jphiN,
142  jphiS,
143  dirs,
144  barrelTopo,
145  endcapTopo,
146  barrelGeom,
147  endcapGeom,
148  debug,
149  ignoreTransition);
150  dirs[0] = SOUTH;
151  std::vector<DetId> vdetS = spr::newECALIdNS(dets,
152  0,
153  jetaE,
154  jetaW,
155  jphiN,
156  jphiS,
157  dirs,
158  barrelTopo,
159  endcapTopo,
160  barrelGeom,
161  endcapGeom,
162  debug,
163  ignoreTransition);
164  for (unsigned int i1 = 0; i1 < vdetS.size(); i1++) {
165  if (std::count(vdets.begin(), vdets.end(), vdetS[i1]) == 0)
166  vdets.push_back(vdetS[i1]);
167  }
168 
169  unsigned int ndet = (ietaE + ietaW + 1) * (iphiN + iphiS + 1);
170  if (vdets.size() != ndet) {
171  std::vector<DetId> vdetExtra;
172  spr::extraIds(det, vdets, ietaE, ietaW, iphiN, iphiS, barrelGeom, endcapGeom, vdetExtra, debug);
173  if (!vdetExtra.empty())
174  vdets.insert(vdets.end(), vdetExtra.begin(), vdetExtra.end());
175  }
176 
177  if (debug) {
178  edm::LogVerbatim("IsoTrack") << "matrixECALIds::Total number of cells found is " << vdets.size();
179  spr::debugEcalDets(0, vdets);
180  }
181  }
182 
183  std::vector<DetId> matrixECALIds(const DetId& det,
184  int ietaE,
185  int ietaW,
186  int iphiN,
187  int iphiS,
188  const CaloGeometry* geo,
189  const CaloTopology* caloTopology,
190  bool debug,
191  bool ignoreTransition) {
192  std::vector<DetId> vdets;
193  spr::matrixECALIds(det, ietaE, ietaW, iphiN, iphiS, geo, caloTopology, vdets, debug, ignoreTransition);
194  return vdets;
195  }
196 
197  std::vector<DetId> newECALIdNS(std::vector<DetId>& dets,
198  unsigned int last,
199  int ieta,
200  int iphi,
201  std::vector<CaloDirection>& dir,
202  const CaloSubdetectorTopology* barrelTopo,
203  const CaloSubdetectorTopology* endcapTopo,
204  const EcalBarrelGeometry* barrelGeom,
205  const EcalEndcapGeometry* endcapGeom,
206  bool debug,
207  bool ignoreTransition) {
208  if (debug) {
209  edm::LogVerbatim("IsoTrack") << "newECALIdNS::Add " << iphi << " columns of cells for " << (dets.size() - last)
210  << " cells (last " << last << ")";
211  spr::debugEcalDets(last, dets, dir);
212  }
213  std::vector<DetId> vdets;
214  std::vector<CaloDirection> dirs;
215  vdets.insert(vdets.end(), dets.begin(), dets.end());
216  dirs.insert(dirs.end(), dir.begin(), dir.end());
217 
218  std::vector<DetId> vdetE, vdetW;
219  if (last == 0) {
220  unsigned int ndet = vdets.size();
221  std::vector<CaloDirection> dirE(ndet, EAST), dirW(ndet, WEST);
222  vdetE = spr::newECALIdEW(
223  dets, last, ieta, dirE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
224  vdetW = spr::newECALIdEW(
225  dets, last, ieta, dirW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
226  for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
227  if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0) {
228  vdets.push_back(vdetW[i1]);
229  dirs.push_back(dir[0]);
230  }
231  }
232  for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
233  if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0) {
234  vdets.push_back(vdetE[i1]);
235  dirs.push_back(dir[0]);
236  }
237  }
238  if (debug) {
239  edm::LogVerbatim("IsoTrack") << "newECALIdNS::With Added cells along E/W results a set of "
240  << (vdets.size() - dets.size()) << " new cells";
241  spr::debugEcalDets(dets.size(), vdets, dirs);
242  }
243  }
244 
245  unsigned int last0 = vdets.size();
246  std::vector<DetId> vdetnew;
247  std::vector<CaloDirection> dirnew;
248  if (iphi > 0) {
249  std::vector<DetId> vdetn(1);
250  std::vector<CaloDirection> dirn(1);
251  std::vector<CaloDirection> dirnE(1, EAST), dirnW(1, WEST);
252  int flag = 0;
253  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
254  std::vector<DetId> cells;
255  spr::simpleMove(
256  dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
257  if (flag != 0) {
258  if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
259  vdetn[0] = cells[0];
260  vdetnew.push_back(vdetn[0]);
261  dirn[0] = dir[i1];
262  if (flag < 0) {
263  if (dirn[0] == NORTH)
264  dirn[0] = SOUTH;
265  else
266  dirn[0] = NORTH;
267  }
268  dirnew.push_back(dirn[0]);
269  vdetE = spr::newECALIdEW(
270  vdetn, 0, ieta, dirnE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
271  vdetW = spr::newECALIdEW(
272  vdetn, 0, ieta, dirnW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
273  for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
274  if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
275  std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0) {
276  vdets.push_back(vdetW[i2]);
277  dirs.push_back(dirn[0]);
278  }
279  }
280  for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
281  if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
282  std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0) {
283  vdets.push_back(vdetE[i2]);
284  dirs.push_back(dirn[0]);
285  }
286  }
287  }
288  }
289  }
290  iphi--;
291  last = vdets.size();
292  for (unsigned int i2 = 0; i2 < vdetnew.size(); i2++) {
293  if (std::count(vdets.begin(), vdets.end(), vdetnew[i2]) == 0) {
294  vdets.push_back(vdetnew[i2]);
295  dirs.push_back(dirnew[i2]);
296  }
297  }
298  if (debug) {
299  edm::LogVerbatim("IsoTrack") << "newECALIdNS::Addition results a set of " << (vdets.size() - last0)
300  << " new cells (last " << last0 << ", iphi " << iphi << ")";
301  spr::debugEcalDets(last0, vdets, dirs);
302  }
303  last0 = last;
304  }
305 
306  if (iphi > 0) {
307  last = last0;
308  return spr::newECALIdNS(
309  vdets, last, ieta, iphi, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
310  } else {
311  if (debug) {
312  edm::LogVerbatim("IsoTrack") << "newECALIdNS::Final list consists of " << vdets.size() << " cells";
313  spr::debugEcalDets(0, vdets);
314  }
315  return vdets;
316  }
317  }
318 
319  std::vector<DetId> newECALIdNS(std::vector<DetId>& dets,
320  unsigned int last,
321  std::vector<int>& ietaE,
322  std::vector<int>& ietaW,
323  std::vector<int>& iphiN,
324  std::vector<int>& iphiS,
325  std::vector<CaloDirection>& dir,
326  const CaloSubdetectorTopology* barrelTopo,
327  const CaloSubdetectorTopology* endcapTopo,
328  const EcalBarrelGeometry* barrelGeom,
329  const EcalEndcapGeometry* endcapGeom,
330  bool debug,
331  bool ignoreTransition) {
332  if (debug) {
333  edm::LogVerbatim("IsoTrack") << "newECALIdNS::Add columns of cells for " << (dets.size() - last)
334  << " cells (last) " << last;
335  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
336  edm::LogVerbatim("IsoTrack") << spr::debugEcalDet(i1, dets[i1]).str() << " along " << dir[i1] << " # "
337  << iphiN[i1] << "|" << iphiS[i1];
338  }
339  }
340  std::vector<DetId> vdets;
341  std::vector<CaloDirection> dirs;
342  std::vector<int> jetaE, jetaW, jphiN, jphiS;
343  vdets.insert(vdets.end(), dets.begin(), dets.end());
344  dirs.insert(dirs.end(), dir.begin(), dir.end());
345  jetaE.insert(jetaE.end(), ietaE.begin(), ietaE.end());
346  jetaW.insert(jetaW.end(), ietaW.begin(), ietaW.end());
347  jphiN.insert(jphiN.end(), iphiN.begin(), iphiN.end());
348  jphiS.insert(jphiS.end(), iphiS.begin(), iphiS.end());
349  std::vector<DetId> vdetE, vdetW;
350  if (last == 0) {
351  unsigned int ndet = vdets.size();
352  std::vector<CaloDirection> dirE(ndet, EAST), dirW(ndet, WEST);
353  vdetE = spr::newECALIdEW(
354  dets, last, ietaE, ietaW, dirE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
355  vdetW = spr::newECALIdEW(
356  dets, last, ietaE, ietaW, dirW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
357  for (unsigned int i1 = 0; i1 < vdetW.size(); i1++) {
358  if (std::count(vdets.begin(), vdets.end(), vdetW[i1]) == 0) {
359  vdets.push_back(vdetW[i1]);
360  dirs.push_back(dir[0]);
361  jetaE.push_back(0);
362  jetaW.push_back(0);
363  jphiN.push_back(iphiN[0]);
364  jphiS.push_back(iphiS[0]);
365  }
366  }
367  for (unsigned int i1 = 0; i1 < vdetE.size(); i1++) {
368  if (std::count(vdets.begin(), vdets.end(), vdetE[i1]) == 0) {
369  vdets.push_back(vdetE[i1]);
370  dirs.push_back(dir[0]);
371  jetaE.push_back(0);
372  jetaW.push_back(0);
373  jphiN.push_back(iphiN[0]);
374  jphiS.push_back(iphiS[0]);
375  }
376  }
377  if (debug) {
378  edm::LogVerbatim("IsoTrack") << "newECALIdNS::With Added cells along E/W results a set of "
379  << (vdets.size() - dets.size()) << " new cells";
380  spr::debugEcalDets(dets.size(), vdets, dirs);
381  }
382  }
383 
384  unsigned int last0 = vdets.size();
385  std::vector<DetId> vdetnew;
386  std::vector<CaloDirection> dirnew;
387  std::vector<int> kphiN, kphiS, ketaE, ketaW;
388  int kphi = 0;
389  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
390  int iphi = iphiS[i1];
391  if (dir[i1] == NORTH)
392  iphi = iphiN[i1];
393  if (iphi > 0) {
394  std::vector<DetId> vdetn(1);
395  std::vector<CaloDirection> dirn(1);
396  std::vector<CaloDirection> dirnE(1, EAST), dirnW(1, WEST);
397  int flag = 0;
398  std::vector<DetId> cells;
399  spr::simpleMove(
400  dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
401  iphi--;
402  if (iphi > kphi)
403  kphi = iphi;
404  if (dir[i1] == NORTH)
405  jphiN[i1] = iphi;
406  else
407  jphiS[i1] = iphi;
408  if (flag != 0) {
409  if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
410  int kfiN = iphiN[i1];
411  int kfiS = iphiS[i1];
412  vdetn[0] = cells[0];
413  vdetnew.push_back(vdetn[0]);
414  dirn[0] = dir[i1];
415  if (dir[i1] == NORTH)
416  kfiN = iphi;
417  else
418  kfiS = iphi;
419  if (flag < 0) {
420  int ktmp = kfiS;
421  kfiS = kfiN;
422  kfiN = ktmp;
423  if (dirn[0] == NORTH)
424  dirn[0] = SOUTH;
425  else
426  dirn[0] = NORTH;
427  }
428  dirnew.push_back(dirn[0]);
429  kphiN.push_back(kfiN);
430  ketaE.push_back(ietaE[i1]);
431  kphiS.push_back(kfiS);
432  ketaW.push_back(ietaW[i1]);
433  std::vector<int> ietE(1, ietaE[i1]), ietW(1, ietaW[i1]);
434  vdetE = spr::newECALIdEW(
435  vdetn, 0, ietE, ietW, dirnE, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
436  vdetW = spr::newECALIdEW(
437  vdetn, 0, ietE, ietW, dirnW, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
438  for (unsigned int i2 = 0; i2 < vdetW.size(); i2++) {
439  if (std::count(vdets.begin(), vdets.end(), vdetW[i2]) == 0 &&
440  std::count(vdetnew.begin(), vdetnew.end(), vdetW[i2]) == 0) {
441  vdets.push_back(vdetW[i2]);
442  dirs.push_back(dirn[0]);
443  jetaE.push_back(0);
444  jphiN.push_back(kfiN);
445  jetaW.push_back(0);
446  jphiS.push_back(kfiS);
447  }
448  }
449  for (unsigned int i2 = 0; i2 < vdetE.size(); i2++) {
450  if (std::count(vdets.begin(), vdets.end(), vdetE[i2]) == 0 &&
451  std::count(vdetnew.begin(), vdetnew.end(), vdetE[i2]) == 0) {
452  vdets.push_back(vdetE[i2]);
453  dirs.push_back(dirn[0]);
454  jetaE.push_back(0);
455  jphiN.push_back(kfiN);
456  jetaW.push_back(0);
457  jphiS.push_back(kfiS);
458  }
459  }
460  }
461  }
462  }
463  }
464  last = vdets.size();
465  for (unsigned int i2 = 0; i2 < vdetnew.size(); i2++) {
466  if (std::count(vdets.begin(), vdets.end(), vdetnew[i2]) == 0) {
467  vdets.push_back(vdetnew[i2]);
468  dirs.push_back(dirnew[i2]);
469  jetaE.push_back(ketaE[i2]);
470  jetaW.push_back(ketaW[i2]);
471  jphiN.push_back(kphiN[i2]);
472  jphiS.push_back(kphiS[i2]);
473  }
474  }
475  if (debug) {
476  edm::LogVerbatim("IsoTrack") << "newECALIdNS::Addition results a set of " << (vdets.size() - last0)
477  << " new cells (last " << last0 << ", iphi " << kphi << ")";
478  for (unsigned int i1 = last0; i1 < vdets.size(); i1++) {
479  edm::LogVerbatim("IsoTrack") << spr::debugEcalDet(i1, vdets[i1]).str() << " along " << dirs[i1] << " iphi "
480  << jphiN[i1] << "|" << jphiS[i1] << " ieta " << jetaE[i1] << "|" << jetaW[i1];
481  }
482  }
483  last0 = last;
484 
485  if (kphi > 0) {
486  last = last0;
487  return spr::newECALIdNS(vdets,
488  last,
489  jetaE,
490  jetaW,
491  jphiN,
492  jphiS,
493  dirs,
494  barrelTopo,
495  endcapTopo,
496  barrelGeom,
497  endcapGeom,
498  debug,
499  ignoreTransition);
500  } else {
501  if (debug) {
502  edm::LogVerbatim("IsoTrack") << "newECALIdNS::Final list consists of " << vdets.size() << " cells";
503  spr::debugEcalDets(0, vdets);
504  }
505  return vdets;
506  }
507  }
508 
509  std::vector<DetId> newECALIdEW(std::vector<DetId>& dets,
510  unsigned int last,
511  int ieta,
512  std::vector<CaloDirection>& dir,
513  const CaloSubdetectorTopology* barrelTopo,
514  const CaloSubdetectorTopology* endcapTopo,
515  const EcalBarrelGeometry* barrelGeom,
516  const EcalEndcapGeometry* endcapGeom,
517  bool debug,
518  bool ignoreTransition) {
519  if (debug) {
520  edm::LogVerbatim("IsoTrack") << "newECALIdEW::Add " << ieta << " rows of cells for " << last << ":" << dets.size()
521  << ":" << (dets.size() - last) << " cells" << std::endl;
522  spr::debugEcalDets(last, dets, dir);
523  }
524  std::vector<DetId> vdets;
525  vdets.clear();
526  std::vector<CaloDirection> dirs;
527  dirs.clear();
528  vdets.insert(vdets.end(), dets.begin(), dets.end());
529  dirs.insert(dirs.end(), dir.begin(), dir.end());
530 
531  if (ieta > 0) {
532  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
533  int flag = 0;
534  std::vector<DetId> cells;
535  spr::simpleMove(
536  dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
537  if (flag != 0) {
538  if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
539  CaloDirection dirn = dir[i1];
540  if (flag < 0) {
541  if (dirn == EAST)
542  dirn = WEST;
543  else
544  dirn = EAST;
545  }
546  vdets.push_back(cells[0]);
547  dirs.push_back(dirn);
548  }
549  }
550  }
551  ieta--;
552  }
553 
554  if (debug) {
555  edm::LogVerbatim("IsoTrack") << "newECALIdEW::Addition results a set of " << (vdets.size() - dets.size())
556  << " new cells";
557  spr::debugEcalDets(dets.size(), vdets, dirs);
558  }
559  if (ieta > 0) {
560  last = dets.size();
561  return spr::newECALIdEW(
562  vdets, last, ieta, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
563  } else {
564  if (debug) {
565  edm::LogVerbatim("IsoTrack") << "newECALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
566  spr::debugEcalDets(0, vdets);
567  }
568  return vdets;
569  }
570  }
571 
572  std::vector<DetId> newECALIdEW(std::vector<DetId>& dets,
573  unsigned int last,
574  std::vector<int>& ietaE,
575  std::vector<int>& ietaW,
576  std::vector<CaloDirection>& dir,
577  const CaloSubdetectorTopology* barrelTopo,
578  const CaloSubdetectorTopology* endcapTopo,
579  const EcalBarrelGeometry* barrelGeom,
580  const EcalEndcapGeometry* endcapGeom,
581  bool debug,
582  bool ignoreTransition) {
583  if (debug) {
584  edm::LogVerbatim("IsoTrack") << "newECALIdEW::Add " << ietaE[0] << "|" << ietaW[0] << " rows of cells for "
585  << (dets.size() - last) << " cells (last " << last << ")";
586  spr::debugEcalDets(last, dets, dir);
587  }
588  std::vector<DetId> vdets;
589  vdets.insert(vdets.end(), dets.begin(), dets.end());
590  std::vector<CaloDirection> dirs;
591  dirs.insert(dirs.end(), dir.begin(), dir.end());
592  std::vector<int> jetaE, jetaW;
593  jetaE.insert(jetaE.end(), ietaE.begin(), ietaE.end());
594  jetaW.insert(jetaW.end(), ietaW.begin(), ietaW.end());
595  int keta = 0;
596  for (unsigned int i1 = last; i1 < dets.size(); i1++) {
597  int ieta = ietaW[i1];
598  if (dir[i1] == EAST)
599  ieta = ietaE[i1];
600  if (ieta > 0) {
601  int flag = 0;
602  std::vector<DetId> cells;
603  spr::simpleMove(
604  dets[i1], dir[i1], barrelTopo, endcapTopo, barrelGeom, endcapGeom, cells, flag, debug, ignoreTransition);
605  ieta--;
606  if (ieta > keta)
607  keta = ieta;
608  if (dir[i1] == EAST)
609  jetaE[i1] = ieta;
610  else
611  jetaW[i1] = ieta;
612  if (flag != 0) {
613  if (std::count(vdets.begin(), vdets.end(), cells[0]) == 0) {
614  vdets.push_back(cells[0]);
615  CaloDirection dirn = dir[i1];
616  int ketaE = ietaE[i1];
617  int ketaW = ietaW[i1];
618  if (dirn == EAST)
619  ketaE = ieta;
620  else
621  ketaW = ieta;
622  if (flag < 0) {
623  int ktmp = ketaW;
624  ketaW = ketaE;
625  ketaE = ktmp;
626  if (dirn == EAST)
627  dirn = WEST;
628  else
629  dirn = EAST;
630  }
631  dirs.push_back(dirn);
632  jetaE.push_back(ketaE);
633  jetaW.push_back(ketaW);
634  }
635  }
636  }
637  }
638 
639  if (debug) {
640  edm::LogVerbatim("IsoTrack") << "newECALIdEW::Addition results a set of " << (vdets.size() - dets.size())
641  << " new cells (last " << dets.size() << ", ieta " << keta << ")";
642  spr::debugEcalDets(dets.size(), vdets);
643  }
644  if (keta > 0) {
645  last = dets.size();
646  return spr::newECALIdEW(
647  vdets, last, jetaE, jetaW, dirs, barrelTopo, endcapTopo, barrelGeom, endcapGeom, debug, ignoreTransition);
648  } else {
649  if (debug) {
650  edm::LogVerbatim("IsoTrack") << "newECALIdEW::Final list (EW) consists of " << vdets.size() << " cells";
651  spr::debugEcalDets(0, vdets);
652  }
653  return vdets;
654  }
655  }
656 
657  void simpleMove(DetId& det,
658  const CaloDirection& dir,
659  const CaloSubdetectorTopology* barrelTopo,
660  const CaloSubdetectorTopology* endcapTopo,
661  const EcalBarrelGeometry* barrelGeom,
662  const EcalEndcapGeometry* endcapGeom,
663  std::vector<DetId>& cells,
664  int& ok,
665  bool debug,
666  bool ignoreTransition) {
667  DetId cell;
668  ok = 0;
669  if (det.subdetId() == EcalBarrel) {
670  EBDetId detId = det;
671  std::vector<DetId> neighbours = barrelTopo->getNeighbours(detId, dir);
672  if (!neighbours.empty() && !neighbours[0].null()) {
673  cells.push_back(neighbours[0]);
674  cell = neighbours[0];
675  ok = 1;
676  } else {
677  const int ietaAbs(detId.ietaAbs()); // abs value of ieta
678  if (EBDetId::MAX_IETA == ietaAbs && (!ignoreTransition) && endcapGeom) {
679  // get ee nbrs for for end of barrel crystals
680  const EcalBarrelGeometry::OrderedListOfEEDetId& ol(*barrelGeom->getClosestEndcapCells(detId));
681  // take closest neighbour on the other side, that is in the endcap
682  cell = *(ol.begin());
683  neighbours = endcapTopo->getNeighbours(cell, dir);
684  if (!neighbours.empty() && !neighbours[0].null())
685  ok = 1;
686  else
687  ok = -1;
688  for (EcalBarrelGeometry::OrderedListOfEEDetId::const_iterator iptr = ol.begin(); iptr != ol.end(); ++iptr)
689  cells.push_back(*iptr);
690  }
691  }
692  } else if (det.subdetId() == EcalEndcap && endcapGeom) {
693  EEDetId detId = det;
694  std::vector<DetId> neighbours = endcapTopo->getNeighbours(detId, dir);
695  if (!neighbours.empty() && !neighbours[0].null()) {
696  cells.push_back(neighbours[0]);
697  cell = neighbours[0];
698  ok = 1;
699  } else {
700  // are we on the outer ring ?
701  const int iphi(detId.iPhiOuterRing());
702  // int isc = detId.isc();
703  if (iphi != 0 && (!ignoreTransition)) {
704  // get eb nbrs for for end of endcap crystals
705  const EcalEndcapGeometry::OrderedListOfEBDetId& ol(*endcapGeom->getClosestBarrelCells(detId));
706  // take closest neighbour on the other side, that is in the barrel.
707  cell = *(ol.begin());
708  neighbours = barrelTopo->getNeighbours(cell, dir);
709  if (!neighbours.empty() && !neighbours[0].null())
710  ok = 1;
711  else
712  ok = -1;
713  for (EcalEndcapGeometry::OrderedListOfEBDetId::const_iterator iptr = ol.begin(); iptr != ol.end(); ++iptr)
714  cells.push_back(*iptr);
715  }
716  }
717  }
718  if (debug) {
719  std::ostringstream st1;
720  for (unsigned int i1 = 0; i1 < cells.size(); ++i1)
721  st1 << " " << std::hex << cells[0]() << std::dec;
722  edm::LogVerbatim("IsoTrack") << "simpleMove:: Move DetId 0x" << std::hex << det() << std::dec << " along " << dir
723  << " to get 0x" << std::hex << cell() << std::dec << " with flag " << ok << " # "
724  << cells.size() << st1.str();
725  }
726  }
727 
728  void extraIds(const DetId& det,
729  std::vector<DetId>& dets,
730  int ietaE,
731  int ietaW,
732  int iphiN,
733  int iphiS,
734  const EcalBarrelGeometry* barrelGeom,
735  const EcalEndcapGeometry* endcapGeom,
736  std::vector<DetId>& cells,
737  bool debug) {
738  if (det.subdetId() == EcalBarrel) {
739  EBDetId id = det;
740  if (debug)
741  edm::LogVerbatim("IsoTrack") << "extraIds::Cell " << id << " rows " << ietaW << "|" << ietaE << " columns "
742  << iphiS << "|" << iphiN;
743  int etaC = id.ietaAbs();
744  int phiC = id.iphi();
745  int zsid = id.zside();
746  for (int eta = -ietaW; eta <= ietaE; ++eta) {
747  for (int phi = -iphiS; phi <= iphiN; ++phi) {
748  int iphi = phiC + phi;
749  if (iphi < 0)
750  iphi += 360;
751  else if (iphi > 360)
752  iphi -= 360;
753  int ieta = zsid * (etaC + eta);
754  if (EBDetId::validDetId(ieta, iphi)) {
755  id = EBDetId(ieta, iphi);
756  if (barrelGeom->present(id)) {
757  if (std::count(dets.begin(), dets.end(), (DetId)id) == 0) {
758  cells.push_back((DetId)id);
759  }
760  }
761  }
762  }
763  }
764  } else if (det.subdetId() == EcalEndcap && endcapGeom) {
765  EEDetId id = det;
766  if (debug)
767  edm::LogVerbatim("IsoTrack") << "extraIds::Cell " << id << " rows " << ietaW << "|" << ietaE << " columns "
768  << iphiS << "|" << iphiN;
769  int ixC = id.ix();
770  int iyC = id.iy();
771  int zsid = id.zside();
772  for (int kx = -ietaW; kx <= ietaE; ++kx) {
773  for (int ky = -iphiS; ky <= iphiN; ++ky) {
774  int ix = ixC + kx;
775  int iy = iyC + ky;
776  if (EEDetId::validDetId(ix, iy, zsid)) {
777  id = EEDetId(ix, iy, zsid);
778  if (endcapGeom->present(id)) {
779  if (std::count(dets.begin(), dets.end(), (DetId)id) == 0) {
780  cells.push_back((DetId)id);
781  }
782  }
783  }
784  }
785  }
786  }
787 
788  if (debug) {
789  edm::LogVerbatim("IsoTrack") << "extraIds:: finds " << cells.size() << " new cells";
790  spr::debugEcalDets(0, cells);
791  }
792  }
793 } // namespace spr
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
EBDetId::ietaAbs
int ietaAbs() const
get the absolute value of the crystal ieta
Definition: EBDetId.h:47
mps_fire.i
i
Definition: mps_fire.py:429
spr::matrixECALIds
void matrixECALIds(const DetId &det, int ieta, int iphi, const CaloGeometry *geo, const CaloTopology *caloTopology, std::vector< DetId > &vdets, bool debug=false, bool igNoreTransition=true)
Definition: MatrixECALDetIds.cc:16
core
Definition: __init__.py:1
testProducerWithPsetDescEmpty_cfi.i2
i2
Definition: testProducerWithPsetDescEmpty_cfi.py:46
EcalBarrelGeometry::getClosestEndcapCells
const OrderedListOfEEDetId * getClosestEndcapCells(EBDetId id) const
Definition: EcalBarrelGeometry.cc:338
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
EBDetId::validDetId
static bool validDetId(int i, int j)
check if a valid index combination
Definition: EBDetId.h:118
spr::debugEcalDets
void debugEcalDets(unsigned int, std::vector< DetId > &)
Definition: DebugInfo.cc:27
spr::newECALIdEW
std::vector< DetId > newECALIdEW(std::vector< DetId > &dets, unsigned int last, int ieta, std::vector< CaloDirection > &dir, const CaloSubdetectorTopology *barrelTopo, const CaloSubdetectorTopology *endcapTopo, const EcalBarrelGeometry *barrelGeom, const EcalEndcapGeometry *endcapGeom, bool debug=false, bool igNoreTransition=true)
Definition: MatrixECALDetIds.cc:509
spr::debugEcalDet
std::ostringstream debugEcalDet(unsigned int, const DetId &)
Definition: DebugInfo.cc:12
spr::newECALIdNS
std::vector< DetId > newECALIdNS(std::vector< DetId > &dets, unsigned int last, int ieta, int iphi, std::vector< CaloDirection > &dir, const CaloSubdetectorTopology *barrelTopo, const CaloSubdetectorTopology *endcapTopo, const EcalBarrelGeometry *barrelGeom, const EcalEndcapGeometry *endcapGeom, bool debug=false, bool igNoreTransition=true)
Definition: MatrixECALDetIds.cc:197
CaloGeometry::getPosition
GlobalPoint getPosition(const DetId &id) const
Get the position of a given detector id.
Definition: CaloGeometry.cc:50
spr::extraIds
void extraIds(const DetId &det, std::vector< DetId > &dets, int ietaE, int ietaW, int iphiN, int iphiS, const EcalBarrelGeometry *barrelGeom, const EcalEndcapGeometry *endcapGeom, std::vector< DetId > &cells, bool debug=false)
Definition: MatrixECALDetIds.cc:728
EcalEndcapGeometry::present
bool present(const DetId &id) const override
is this detid present in the geometry?
Definition: EcalEndcapGeometry.cc:450
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
CaloSubdetectorTopology::getNeighbours
virtual std::vector< DetId > getNeighbours(const DetId &id, const CaloDirection &dir) const
Definition: CaloSubdetectorTopology.h:77
cells
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ cells
Definition: CAHitNtupletGeneratorKernelsImpl.h:43
DetId
Definition: DetId.h:17
debug
#define debug
Definition: HDRShower.cc:19
CaloTopology::getSubdetectorTopology
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
Definition: CaloTopology.cc:17
EEDetId::validDetId
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:248
EBDetId::MAX_IETA
static const int MAX_IETA
Definition: EBDetId.h:136
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
testProducerWithPsetDescEmpty_cfi.i1
i1
Definition: testProducerWithPsetDescEmpty_cfi.py:45
EcalBarrelGeometry::present
bool present(const DetId &id) const override
is this detid present in the geometry?
Definition: EcalBarrelGeometry.cc:441
EZArrayFL::const_iterator
MgrType::const_iterator const_iterator
Definition: EZArrayFL.h:24
EcalEndcapGeometry::getClosestBarrelCells
const OrderedListOfEBDetId * getClosestBarrelCells(EEDetId id) const
Definition: EcalEndcapGeometry.cc:372
heppy_check.dirs
dirs
Definition: heppy_check.py:26
CaloDirection
CaloDirection
Codes the local directions in the cell lattice.
Definition: CaloDirection.h:9
spr::simpleMove
void simpleMove(DetId &det, const CaloDirection &dir, const CaloSubdetectorTopology *barrelTopo, const CaloSubdetectorTopology *endcapTopo, const EcalBarrelGeometry *barrelGeom, const EcalEndcapGeometry *endcapGeom, std::vector< DetId > &cells, int &flag, bool debug=false, bool igNoreTransition=true)
Definition: MatrixECALDetIds.cc:657
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:34
EEDetId::iPhiOuterRing
int iPhiOuterRing() const
Definition: EEDetId.cc:295
Generated for CMSSW Reference Manual by doxygen 1.8.14