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