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