CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalNumberingFromDDD.cc
Go to the documentation of this file.
1 // File: HcalNumberingFromDDD.cc
3 // Description: Usage of DDD to get to numbering scheme for hadron calorimeter
5 
6 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
7 
8 #include "DetectorDescription/Base/interface/DDException.h"
9 #include "DetectorDescription/Core/interface/DDFilter.h"
10 #include "DetectorDescription/Core/interface/DDFilteredView.h"
11 #include "DetectorDescription/Core/interface/DDSolid.h"
12 #include "DetectorDescription/Core/interface/DDValue.h"
13 
14 #include "FWCore/MessageLogger/interface/MessageLogger.h"
15 
16 #include "CLHEP/Units/GlobalPhysicalConstants.h"
17 #include "CLHEP/Units/GlobalSystemOfUnits.h"
18 #include <iostream>
19 
20 //#define DebugLog
21 
22 HcalNumberingFromDDD::HcalNumberingFromDDD(std::string & name,
23  const DDCompactView & cpv) {
24  edm::LogInfo("HCalGeom") << "Creating HcalNumberingFromDDD";
25  initialize(name, cpv);
26 }
27 
28 HcalNumberingFromDDD::~HcalNumberingFromDDD() {
29  edm::LogInfo("HCalGeom") << "Deleting HcalNumberingFromDDD";
30 }
31 
32 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
33  CLHEP::Hep3Vector point,
34  int depth,
35  int lay) const {
36 
37 
38  double hx = point.x();
39  double hy = point.y();
40  double hz = point.z();
41  double hR = sqrt(hx*hx+hy*hy+hz*hz);
42  double htheta = (hR == 0. ? 0. : acos(std::max(std::min(hz/hR,1.0),-1.0)));
43  double hsintheta = sin(htheta);
44  double hphi = (hR*hsintheta == 0. ? 0. :atan2(hy,hx));
45  double heta = (fabs(hsintheta) == 1.? 0. : -log(fabs(tan(htheta/2.))) );
46 
47  int hsubdet=0;
48  double etaR;
49 
50  //First eta index
51  if (det == 5) { // Forward HCal
52  hsubdet = static_cast<int>(HcalForward);
53  hR = sqrt(hx*hx+hy*hy);
54  etaR = (heta >= 0. ? hR : -hR);
55  } else { // Barrel or Endcap
56  etaR = heta;
57  if (det == 3) {
58  hsubdet = static_cast<int>(HcalBarrel);
59  } else {
60  hsubdet = static_cast<int>(HcalEndcap);
61  }
62  }
63 
64 #ifdef DebugLog
65  LogDebug("HCalGeom") << "HcalNumberingFromDDD: point = " << point << " det "
66  << hsubdet << " eta/R " << etaR << " phi " << hphi;
67 #endif
68  HcalNumberingFromDDD::HcalID tmp = unitID(hsubdet,etaR,hphi,depth,lay);
69  return tmp;
70 
71 }
72 
73 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(double eta,double fi,
74  int depth,
75  int lay) const {
76 
77  int ieta = 0;
78  double heta = fabs(eta);
79  for (int i = 0; i < nEta; i++)
80  if (heta > etaTable[i]) ieta = i + 1;
81  int hsubdet=0;
82  double etaR;
83  if (ieta <= etaMin[1]) {
84  if ((ieta <= etaMin[1] && depth==3) || ieta > etaMax[0]) {
85  hsubdet = static_cast<int>(HcalEndcap);
86  } else {
87  hsubdet = static_cast<int>(HcalBarrel);
88  }
89  etaR = eta;
90  } else {
91  hsubdet = static_cast<int>(HcalForward);
92  double theta = 2.*atan(exp(-heta));
93  double hR = zVcal*tan(theta);
94  etaR = (eta >= 0. ? hR : -hR);
95  }
96 
97  HcalNumberingFromDDD::HcalID tmp = unitID(hsubdet,etaR,fi,depth,lay);
98  return tmp;
99 }
100 
101 
102 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
103  double etaR,
104  double phi,
105  int depth,
106  int lay) const {
107 
108  int ieta=0;
109  double fioff, fibin;
110  double hetaR = fabs(etaR);
111 
112  //First eta index
113  if (det == static_cast<int>(HcalForward)) { // Forward HCal
114  fioff = phioff[2];
115  ieta = etaMax[2];
116  for (int i = nR-1; i > 0; i--)
117  if (hetaR < rTable[i]) ieta = etaMin[2] + nR - i - 1;
118  fibin = phibin[nEta+ieta-etaMin[2]-1];
119  if (ieta > etaMax[2]-2 ) { // HF double-phi
120  fioff += 0.5*fibin;
121  }
122  } else { // Barrel or Endcap
123  ieta = 1;
124  for (int i = 0; i < nEta-1; i++)
125  if (hetaR > etaTable[i]) ieta = i + 1;
126  if (det == static_cast<int>(HcalBarrel)) {
127  fioff = phioff[0];
128  if (ieta > etaMax[0]) ieta = etaMax[0];
129  if (lay == 18 && nOff.size() > 13) {
130  if (hetaR > etaHO[1] && ieta == nOff[13]) ieta++;
131  }
132  } else {
133  fioff = phioff[1];
134  if (ieta <= etaMin[1]) ieta = etaMin[1];
135  }
136  fibin = phibin[ieta-1];
137  }
138 
139  int nphi = int((CLHEP::twopi+0.1*fibin)/fibin);
140  int zside = etaR>0 ? 1: 0;
141  double hphi = phi+fioff;
142  if (hphi < 0) hphi += CLHEP::twopi;
143  int iphi = int(hphi/fibin) + 1;
144  if (iphi > nphi) iphi = 1;
145 
146 #ifdef DebugLog
147  LogDebug("HCalGeom") << "HcalNumberingFromDDD: etaR = " << etaR << " : "
148  << zside << "/" << ieta << " phi " << hphi << " : "
149  << iphi;
150 #endif
151  HcalNumberingFromDDD::HcalID tmp = unitID(det,zside,depth,ieta,iphi,lay);
152  return tmp;
153 
154 }
155 
156 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det, int zside,
157  int depth, int etaR,
158  int phi,
159  int lay) const {
160 
161  //Modify the depth index
162  if (det == static_cast<int>(HcalForward)) { // Forward HCal
163  } else {
164  if (lay >= 0) {
165  double fibin = phibin[etaR-1];
166  int depth0 = depth1[etaR-1];
167  int kphi = phi + int((phioff[3]+0.1)/fibin);
168  kphi = (kphi-1)%4 + 1;
169  if (etaR == nOff[0] && (kphi == 2 || kphi == 3)) depth0--;
170  if (lay <= depth2[etaR-1]) {
171  if (lay <= depth0) depth = 1;
172  else depth = 2;
173  } else if (lay <= depth3[etaR-1]) {
174  depth = 3;
175  } else depth = 4;
176  } else if (det == static_cast<int>(HcalBarrel)) {
177  if (depth==3) depth = 2;
178  }
179  if (det != static_cast<int>(HcalBarrel)) {
180  if (etaR <= etaMin[1]) depth = 3;
181  }
182  }
183  if (etaR == nOff[1] && depth > 2 && det == static_cast<int>(HcalEndcap))
184  etaR = nOff[1]-1;
185  if (det == static_cast<int>(HcalBarrel) && depth == 4) {
186  det = static_cast<int>(HcalOuter);
187  }
188 
189  int units = unitPhi(det, etaR);
190  int iphi_skip = phi;
191  if (units==2) iphi_skip = (phi-1)*2+1;
192  else if (units==4) iphi_skip = (phi-1)*4-1;
193  if (iphi_skip < 0) iphi_skip += 72;
194 
195 #ifdef DebugLog
196  LogDebug("HCalGeom") << "HcalNumberingFromDDD: phi units=" << units
197  << " iphi_skip=" << iphi_skip;
198 #endif
199  HcalNumberingFromDDD::HcalID tmp(det,zside,depth,etaR,phi,iphi_skip,lay);
200 
201 #ifdef DebugLog
202  LogDebug("HCalGeom") << "HcalNumberingFromDDD: det = " << det << " "
203  << tmp.subdet << " zside = " << tmp.zside << " depth = "
204  << tmp.depth << " eta/R = " << tmp.etaR << " phi = "
205  << tmp.phi << " layer = " << tmp.lay;
206 #endif
207  return tmp;
208 }
209 
210 HcalCellType::HcalCell HcalNumberingFromDDD::cell(int det, int zside,
211  int depth, int etaR,
212  int iphi, bool corr) const {
213 
214  int idet = det;
215  double etaMn = etaMin[0];
216  double etaMx = etaMax[0];
217  if (idet==static_cast<int>(HcalEndcap)) {
218  etaMn = etaMin[1]; etaMx = etaMax[1];
219  } else if (idet==static_cast<int>(HcalForward)) {
220  etaMn = etaMin[2]; etaMx = etaMax[2];
221  }
222  if (corr) {
223  if (etaR >= nOff[2] && depth == 3 && idet == static_cast<int>(HcalBarrel))
224  idet = static_cast<int>(HcalEndcap);
225  }
226  double eta = 0, deta = 0, phi = 0, dphi = 0, rz = 0, drz = 0;
227  bool ok = false, flagrz = true;
228  if ((idet==static_cast<int>(HcalBarrel)||idet==static_cast<int>(HcalEndcap)||
229  idet==static_cast<int>(HcalOuter)||idet==static_cast<int>(HcalForward))
230  && etaR >=etaMn && etaR <= etaMx)
231  ok = true;
232  if (idet == static_cast<int>(HcalEndcap)) {
233  if (depth < 3 && etaR <= etaMin[1]) ok = false;
234  else if (depth > 2 && etaR == nOff[1]) ok = false;
235  }
236  if (ok) {
237  int maxlay = (int)(rHB.size());
238  if (idet == static_cast<int>(HcalEndcap)) maxlay = (int)(zHE.size());
239  eta = getEta(idet, etaR, zside, depth);
240  deta = deltaEta(idet, etaR, depth);
241  double fibin, fioff;
242  if (idet == static_cast<int>(HcalBarrel)||
243  idet == static_cast<int>(HcalOuter)) {
244  fioff = phioff[0];
245  fibin = phibin[etaR-1];
246  } else if (idet == static_cast<int>(HcalEndcap)) {
247  fioff = phioff[1];
248  fibin = phibin[etaR-1];
249  } else {
250  fioff = phioff[2];
251  fibin = phibin[nEta+etaR-etaMin[2]-1];
252  if (etaR > etaMax[2]-2 ) fioff += 0.5*fibin;
253  }
254  phi = fioff + (iphi - 0.5)*fibin;
255  dphi = 0.5*fibin;
256  if (idet == static_cast<int>(HcalForward)) {
257  int ir = nR + etaMin[2] - etaR - 1;
258  if (ir > 0 && ir < nR) {
259  rz = 0.5*(rTable[ir]+rTable[ir-1]);
260  drz = 0.5*(rTable[ir]-rTable[ir-1]);
261  } else {
262  ok = false;
263 #ifdef DebugLog
264  LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong eta " << etaR
265  << " (" << ir << "/" << nR << ") Detector "
266  << idet;
267 #endif
268  }
269  if (depth != 1 && depth != 2) {
270  ok = false;
271 #ifdef DebugLog
272  LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong depth " << depth
273  << " in Detector " << idet;
274 #endif
275  }
276  } else if (etaR <= nEta) {
277  int depth0 = depth1[etaR-1];
278  int kphi = iphi + int((phioff[3]+0.1)/fibin);
279  kphi = (kphi-1)%4 + 1;
280  if (etaR == nOff[0] && (kphi == 2 || kphi == 3)) depth0--;
281  int laymin, laymax;
282  if (depth == 1) {
283  laymin = 1;
284  if (idet==static_cast<int>(HcalEndcap)) laymin = 2;
285  laymax = depth0;
286  if (nOff.size() > 12) {
287  if (etaR == nOff[6]) {
288  laymin = nOff[7];
289  laymax = nOff[8];
290  } else if (etaR == nOff[9]) {
291  laymin = nOff[10];
292  }
293  }
294  } else if (depth == 2) {
295  laymin = depth0+1;
296  laymax = depth2[etaR-1];
297  if (etaR==etaMax[0] && idet==static_cast<int>(HcalBarrel) &&
298  nOff.size()>3) laymax = nOff[3];
299  if (nOff.size() > 12) {
300  if (etaR == nOff[9]) laymax = nOff[11];
301  if (etaR == nOff[6]) laymax = nOff[12];
302  }
303  } else if (depth == 3) {
304  laymin = depth2[etaR-1]+1;
305  laymax = depth3[etaR-1];
306  if (etaR<=etaMin[1] && idet==static_cast<int>(HcalEndcap)) {
307  if (nOff.size() > 4) laymin = nOff[4];
308  if (nOff.size() > 5) laymax = nOff[5];
309  }
310  } else {
311  laymin = depth3[etaR-1]+1;
312  laymax = maxlay;
313  }
314  if (idet == static_cast<int>(HcalOuter) && nOff.size() > 13) {
315  if (etaR > nOff[13] && laymin <= laymax) laymin = laymax;
316  }
317  double d1=0, d2=0;
318  if (laymin <= maxlay && laymax <= maxlay && laymin <= laymax) {
319  if (idet == static_cast<int>(HcalEndcap)) {
320  flagrz = false;
321  if (depth == 1 || laymin <= 1) d1 = zHE[laymin-1] - dzHE[laymin-1];
322  else d1 = zHE[laymin-2] + dzHE[laymin-2];
323  d2 = zHE[laymax-1] + dzHE[laymax-1];
324  } else {
325  if (idet == static_cast<int>(HcalOuter) ||
326  depth == 1 || laymin <=1) d1 = rHB[laymin-1] - drHB[laymin-1];
327  else d1 = rHB[laymin-2] + drHB[laymin-1];
328  d2 = rHB[laymax-1] + drHB[laymax-1];
329  }
330  rz = 0.5*(d2+d1);
331  drz = 0.5*(d2-d1);
332  } else {
333  ok = false;
334 #ifdef DebugLog
335  LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong depth " << depth
336  << " (Layer minimum " << laymin << " maximum "
337  << laymax << " maxLay " << maxlay << ")";
338 #endif
339  }
340  } else {
341  ok = false;
342 #ifdef DebugLog
343  LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong eta " << etaR
344  << "/" << nEta << " Detector " << idet;
345 #endif
346  }
347  } else {
348  ok = false;
349 #ifdef DebugLog
350  LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong eta " << etaR
351  << " det " << idet;
352 #endif
353  }
354  HcalCellType::HcalCell tmp(ok,eta,deta,phi,dphi,rz,drz,flagrz);
355 
356 #ifdef DebugLog
357  LogDebug("HCalGeom") << "HcalNumberingFromDDD: det/side/depth/etaR/phi "
358  << det << "/" << zside << "/" << depth << "/" << etaR
359  << "/" << iphi << " Cell Flag " << tmp.ok << " "
360  << tmp.eta << " " << tmp.deta << " phi " << tmp.phi
361  << " " << tmp.dphi << " r(z) " << tmp.rz << " "
362  << tmp.drz << " " << tmp.flagrz;
363 #endif
364  return tmp;
365 }
366 
367 std::vector<double> HcalNumberingFromDDD::getEtaTable() const {
368 
369  std::vector<double> tmp = etaTable;
370  return tmp;
371 }
372 
373 unsigned int HcalNumberingFromDDD::numberOfCells(HcalSubdetector subdet) const{
374 
375  unsigned int num = 0;
376  std::vector<HcalCellType> cellTypes = HcalCellTypes(subdet);
377  for (unsigned int i=0; i<cellTypes.size(); i++) {
378  num += (unsigned int)(cellTypes[i].nPhiBins());
379  if (cellTypes[i].nHalves() > 1)
380  num += (unsigned int)(cellTypes[i].nPhiBins());
381  num -= (unsigned int)(cellTypes[i].nPhiMissingBins());
382  }
383 #ifdef DebugLog
384  LogDebug ("HCalGeom") << "HcalNumberingFromDDD:numberOfCells "
385  << cellTypes.size() << " " << num
386  << " for subdetector " << subdet;
387 #endif
388  return num;
389 }
390 
391 std::vector<HcalCellType> HcalNumberingFromDDD::HcalCellTypes() const{
392 
393  std::vector<HcalCellType> cellTypes =HcalCellTypes(HcalBarrel);
394 #ifdef DebugLog
395  LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << cellTypes.size()
396  << " cells of type HCal Barrel";
397  for (unsigned int i=0; i<cellTypes.size(); i++)
398  LogDebug ("HCalGeom") << "Cell " << i << " " << cellTypes[i];
399 #endif
400 
401  std::vector<HcalCellType> hoCells =HcalCellTypes(HcalOuter);
402 #ifdef DebugLog
403  LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << hoCells.size()
404  << " cells of type HCal Outer";
405  for (unsigned int i=0; i<hoCells.size(); i++)
406  LogDebug ("HCalGeom") << "Cell " << i << " " << hoCells[i];
407 #endif
408  cellTypes.insert(cellTypes.end(), hoCells.begin(), hoCells.end());
409 
410  std::vector<HcalCellType> heCells =HcalCellTypes(HcalEndcap);
411 #ifdef DebugLog
412  LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << heCells.size()
413  << " cells of type HCal Endcap";
414  for (unsigned int i=0; i<heCells.size(); i++)
415  LogDebug ("HCalGeom") << "Cell " << i << " " << heCells[i];
416 #endif
417  cellTypes.insert(cellTypes.end(), heCells.begin(), heCells.end());
418 
419  std::vector<HcalCellType> hfCells =HcalCellTypes(HcalForward);
420 #ifdef DebugLog
421  LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << hfCells.size()
422  << " cells of type HCal Forward";
423  for (unsigned int i=0; i<hfCells.size(); i++)
424  LogDebug ("HCalGeom") << "Cell " << i << " " << hfCells[i];
425 #endif
426  cellTypes.insert(cellTypes.end(), hfCells.begin(), hfCells.end());
427 
428  return cellTypes;
429 }
430 
431 std::vector<HcalCellType> HcalNumberingFromDDD::HcalCellTypes(HcalSubdetector subdet) const {
432 
433  std::vector<HcalCellType> cellTypes;
434  if (subdet == HcalForward) {
435  if (dzVcal < 0) return cellTypes;
436  }
437 
438  int dmin, dmax, indx, nz, nmod;
439  double hsize = 0;
440  switch(subdet) {
441  case HcalEndcap:
442  dmin = 1; dmax = 3; indx = 1; nz = nzHE; nmod = nmodHE;
443  break;
444  case HcalForward:
445  dmin = 1; dmax = 2; indx = 2; nz = 2; nmod = 18;
446  break;
447  case HcalOuter:
448  dmin = 4; dmax = 4; indx = 0; nz = nzHB; nmod = nmodHB;
449  break;
450  default:
451  dmin = 1; dmax = 3; indx = 0; nz = nzHB; nmod = nmodHB;
452  break;
453  }
454 
455  int phi = 1, zside = 1;
456  bool cor = false;
457 
458  // Get the Cells
459  int subdet0 = static_cast<int>(subdet);
460  for (int depth=dmin; depth<=dmax; depth++) {
461  int shift = getShift(subdet, depth);
462  double gain = getGain (subdet, depth);
463  if (subdet == HcalForward) {
464  if (depth == 1) hsize = dzVcal;
465  else hsize = dzVcal-0.5*dlShort;
466  }
467  for (int eta=etaMin[indx]; eta<= etaMax[indx]; eta++) {
468  HcalCellType::HcalCell temp1 = cell(subdet0,zside,depth,eta,phi,cor);
469  if (temp1.ok) {
470  int units = unitPhi (subdet0, eta);
471  HcalCellType temp2(subdet, eta, phi, depth, temp1,
472  shift, gain, nz, nmod, hsize, units);
473  if (subdet == HcalOuter && nOff.size() > 17) {
474  if (eta == nOff[15]) {
475  std::vector<int> missPlus, missMinus;
476  int kk = 18;
477  for (int miss=0; miss<nOff[16]; miss++) {
478  missPlus.push_back(nOff[kk]);
479  kk++;
480  }
481  for (int miss=0; miss<nOff[17]; miss++) {
482  missMinus.push_back(nOff[kk]);
483  kk++;
484  }
485  temp2.setMissingPhi(missPlus, missMinus);
486  }
487  }
488  cellTypes.push_back(temp2);
489  }
490  }
491  }
492  return cellTypes;
493 }
494 
495 void HcalNumberingFromDDD::printTile() {
496 
497  std::cout << "Tile Information for HB:\n" << "========================\n\n";
498  for (int eta=etaMin[0]; eta<= etaMax[0]; eta++) {
499  int dmax = 1;
500  if (depth1[eta-1] < 17) dmax = 2;
501  for (int depth=1; depth<=dmax; depth++)
502  tileHB(eta, depth);
503  }
504 
505  std::cout << "\nTile Information for HE:\n" <<"========================\n\n";
506  for (int eta=etaMin[1]; eta<= etaMax[1]; eta++) {
507  int dmin=1, dmax=3;
508  if (eta == etaMin[1]) {
509  dmin = 3;
510  } else if (depth1[eta-1] > 18) {
511  dmax = 1;
512  } else if (depth2[eta-1] > 18) {
513  dmax = 2;
514  }
515  for (int depth=dmin; depth<=dmax; depth++)
516  tileHE(eta, depth);
517  }
518 }
519 
520 double HcalNumberingFromDDD::getEta(int det, int etaR, int zside,
521  int depth) const {
522 
523  double tmp = 0;
524  if (det == static_cast<int>(HcalForward)) {
525  int ir = nR + etaMin[2] - etaR - 1;
526  if (ir > 0 && ir < nR) {
527  double z = zVcal;
528  if (depth != 1) z += dlShort;
529  tmp = 0.5*(getEta(rTable[ir-1],z)+getEta(rTable[ir],z));
530  }
531  } else {
532  if (etaR > 0 && etaR < nEta) {
533  if (etaR == nOff[1]-1 && depth > 2) {
534  tmp = 0.5*(etaTable[etaR+1]+etaTable[etaR-1]);
535  } else if (det == static_cast<int>(HcalOuter) && nOff.size() > 13) {
536  if (etaR == nOff[13]) {
537  tmp = 0.5*(etaHO[0]+etaTable[etaR-1]);
538  } else if (etaR == nOff[13]+1) {
539  tmp = 0.5*(etaTable[etaR]+etaHO[1]);
540  } else if (etaR == nOff[14]) {
541  tmp = 0.5*(etaHO[2]+etaTable[etaR-1]);
542  } else if (etaR == nOff[14]+1) {
543  tmp = 0.5*(etaTable[etaR]+etaHO[3]);
544  } else {
545  tmp = 0.5*(etaTable[etaR]+etaTable[etaR-1]);
546  }
547  } else {
548  tmp = 0.5*(etaTable[etaR]+etaTable[etaR-1]);
549  }
550  }
551  }
552  if (zside == 0) tmp = -tmp;
553 #ifdef DebugLog
554  LogDebug("HCalGeom") << "HcalNumberingFromDDD::getEta " << etaR << " "
555  << zside << " " << depth << " ==> " << tmp;
556 #endif
557  return tmp;
558 }
559 
560 double HcalNumberingFromDDD::getEta(double r, double z) const {
561 
562  double tmp = 0;
563  if (z != 0) tmp = -log(tan(0.5*atan(r/z)));
564 #ifdef DebugLog
565  LogDebug("HCalGeom") << "HcalNumberingFromDDD::getEta " << r << " " << z
566  << " ==> " << tmp;
567 #endif
568  return tmp;
569 }
570 
571 double HcalNumberingFromDDD::deltaEta(int det, int etaR, int depth) const {
572 
573  double tmp = 0;
574  if (det == static_cast<int>(HcalForward)) {
575  int ir = nR + etaMin[2] - etaR - 1;
576  if (ir > 0 && ir < nR) {
577  double z = zVcal;
578  if (depth != 1) z += dlShort;
579  tmp = 0.5*(getEta(rTable[ir-1],z)-getEta(rTable[ir],z));
580  }
581  } else {
582  if (etaR > 0 && etaR < nEta) {
583  if (etaR == nOff[1]-1 && depth > 2) {
584  tmp = 0.5*(etaTable[etaR+1]-etaTable[etaR-1]);
585  } else if (det == static_cast<int>(HcalOuter) && nOff.size() > 13) {
586  if (etaR == nOff[13]) {
587  tmp = 0.5*(etaHO[0]-etaTable[etaR-1]);
588  } else if (etaR == nOff[13]+1) {
589  tmp = 0.5*(etaTable[etaR]-etaHO[1]);
590  } else if (etaR == nOff[14]) {
591  tmp = 0.5*(etaHO[2]-etaTable[etaR-1]);
592  } else if (etaR == nOff[14]+1) {
593  tmp = 0.5*(etaTable[etaR]-etaHO[3]);
594  } else {
595  tmp = 0.5*(etaTable[etaR]-etaTable[etaR-1]);
596  }
597  } else {
598  tmp = 0.5*(etaTable[etaR]-etaTable[etaR-1]);
599  }
600  }
601  }
602 #ifdef DebugLog
603  LogDebug("HCalGeom") << "HcalNumberingFromDDD::deltaEta " << etaR << " "
604  << depth << " ==> " << tmp;
605 #endif
606  return tmp;
607 }
608 
609 void HcalNumberingFromDDD::initialize(std::string & name,
610  const DDCompactView & cpv) {
611 
612  std::string attribute = "ReadOutName";
613  edm::LogInfo("HCalGeom") << "HcalNumberingFromDDD: Initailise for " << name
614  << " as " << attribute;
615 
616  DDSpecificsFilter filter;
617  DDValue ddv(attribute,name,0);
618  filter.setCriteria(ddv,DDSpecificsFilter::equals);
619  DDFilteredView fv(cpv);
620  fv.addFilter(filter);
621  bool ok = fv.firstChild();
622 
623  if (ok) {
624  //Load the SpecPars
625  loadSpecPars(fv);
626 
627  //Load the Geometry parameters
628  loadGeometry(fv);
629  } else {
630  edm::LogError("HCalGeom") << "HcalNumberingFromDDD: cannot get filtered "
631  << " view for " << attribute << " matching "
632  << name;
633  throw DDException("HcalNumberingFromDDD: cannot match "+attribute+" to "+name);
634  }
635 
636 #ifdef DebugLog
637  std::vector<HcalCellType> cellTypes = HcalCellTypes();
638  LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << cellTypes.size()
639  << " cells of type HCal (All)";
640  for (unsigned int i=0; i<cellTypes.size(); i++)
641  LogDebug ("HCalGeom") << "Cell " << i << " " << cellTypes[i];
642 #endif
643 }
644 
645 void HcalNumberingFromDDD::loadSpecPars(DDFilteredView fv) {
646 
647  DDsvalues_type sv(fv.mergedSpecifics());
648 
649  // Phi Offset
650  int i, nphi=4;
651  std::vector<double> tmp1 = getDDDArray("phioff",sv,nphi);
652  phioff.resize(tmp1.size());
653  for (i=0; i<nphi; i++) {
654  phioff[i] = tmp1[i];
655 #ifdef DebugLog
656  LogDebug("HCalGeom") << "HcalNumberingFromDDD: phioff[" << i << "] = "
657  << phioff[i]/CLHEP::deg;
658 #endif
659  }
660 
661  //Eta table
662  nEta = -1;
663  std::vector<double> tmp2 = getDDDArray("etaTable",sv,nEta);
664  etaTable.resize(tmp2.size());
665  for (i=0; i<nEta; i++) {
666  etaTable[i] = tmp2[i];
667 #ifdef DebugLog
668  LogDebug("HCalGeom") << "HcalNumberingFromDDD: etaTable[" << i << "] = "
669  << etaTable[i];
670 #endif
671  }
672 
673  //R table
674  nR = -1;
675  std::vector<double> tmp3 = getDDDArray("rTable",sv,nR);
676  rTable.resize(tmp3.size());
677  for (i=0; i<nR; i++) {
678  rTable[i] = tmp3[i];
679 #ifdef DebugLog
680  LogDebug("HCalGeom") << "HcalNumberingFromDDD: rTable[" << i << "] = "
681  << rTable[i]/CLHEP::cm;
682 #endif
683  }
684 
685  //Phi bins
686  nPhi = nEta + nR - 2;
687  std::vector<double> tmp4 = getDDDArray("phibin",sv,nPhi);
688  phibin.resize(tmp4.size());
689  for (i=0; i<nPhi; i++) {
690  phibin[i] = tmp4[i];
691 #ifdef DebugLog
692  LogDebug("HCalGeom") << "HcalNumberingFromDDD: phibin[" << i << "] = "
693  << phibin[i]/CLHEP::deg;
694 #endif
695  }
696 
697  //Layer boundaries for depths 1, 2, 3, 4
698  nDepth = nEta - 1;
699  std::vector<double> d1 = getDDDArray("depth1",sv,nDepth);
700  nDepth = nEta - 1;
701  std::vector<double> d2 = getDDDArray("depth2",sv,nDepth);
702  nDepth = nEta - 1;
703  std::vector<double> d3 = getDDDArray("depth3",sv,nDepth);
704 #ifdef DebugLog
705  LogDebug("HCalGeom") << "HcalNumberingFromDDD: " << nDepth << " Depths";
706 #endif
707  depth1.resize(nDepth);
708  depth2.resize(nDepth);
709  depth3.resize(nDepth);
710  for (i=0; i<nDepth; i++) {
711  depth1[i] = static_cast<int>(d1[i]);
712  depth2[i] = static_cast<int>(d2[i]);
713  depth3[i] = static_cast<int>(d3[i]);
714 #ifdef DebugLog
715  LogDebug("HCalGeom") << "HcalNumberingFromDDD: depth1[" << i << "] = "
716  << depth1[i] << " depth2[" << i << "] = "<< depth2[i]
717  << " depth3[" << i << "] = " << depth3[i];
718 #endif
719  }
720 
721  // Minimum and maximum eta boundaries
722  int ndx = 3;
723  std::vector<double> tmp5 = getDDDArray("etaMin",sv,ndx);
724  std::vector<double> tmp6 = getDDDArray("etaMax",sv,ndx);
725  etaMin.resize(ndx);
726  etaMax.resize(ndx);
727  for (i=0; i<ndx; i++) {
728  etaMin[i] = static_cast<int>(tmp5[i]);
729  etaMax[i] = static_cast<int>(tmp6[i]);
730  }
731  etaMin[0] = 1;
732  etaMax[1] = nEta-1;
733  etaMax[2] = etaMin[2]+nR-2;
734 #ifdef DebugLog
735  for (i=0; i<ndx; i++)
736  LogDebug("HCalGeom") << "HcalNumberingFromDDD: etaMin[" << i << "] = "
737  << etaMin[i] << " etaMax[" << i << "] = "<< etaMax[i];
738 #endif
739 
740  // Geometry parameters for HF
741  int ngpar = 7;
742  std::vector<double> gpar = getDDDArray("gparHF",sv,ngpar);
743  dlShort = gpar[0];
744  zVcal = gpar[4];
745 #ifdef DebugLog
746  LogDebug("HCalGeom") << "HcalNumberingFromDDD: dlShort " << dlShort
747  << " zVcal " << zVcal;
748 #endif
749 
750  // nOff
751  int noff = 3;
752  std::vector<double> nvec = getDDDArray("noff",sv,noff);
753  nOff.resize(noff);
754  for (i=0; i<noff; i++) {
755  nOff[i] = static_cast<int>(nvec[i]);
756 #ifdef DebugLog
757  LogDebug("HCalGeom") << "HcalNumberingFromDDD: nOff[" << i << "] = "
758  << nOff[i];
759 #endif
760  }
761 
762  //Gains and Shifts for HB depths
763  ndx = 4;
764  gainHB = getDDDArray("HBGains",sv,ndx);
765  std::vector<double> tmp7 = getDDDArray("HBShift",sv,ndx);
766  shiftHB.resize(ndx);
767 #ifdef DebugLog
768  LogDebug("HCalGeom") << "HcalNumberingFromDDD:: Gain factor and Shift for "
769  << "HB depth layers:";
770 #endif
771  for (i=0; i<ndx; i++) {
772  shiftHB[i] = static_cast<int>(tmp7[i]);
773 #ifdef DebugLog
774  LogDebug("HCalGeom") <<"HcalNumberingFromDDD:: gainHB[" << i << "] = "
775  << gainHB[i] << " shiftHB[" << i << "] = "
776  << shiftHB[i];
777 #endif
778  }
779 
780  //Gains and Shifts for HB depths
781  ndx = 4;
782  gainHE = getDDDArray("HEGains",sv,ndx);
783  std::vector<double> tmp8 = getDDDArray("HEShift",sv,ndx);
784  shiftHE.resize(ndx);
785 #ifdef DebugLog
786  LogDebug("HCalGeom") << "HcalNumberingFromDDD:: Gain factor and Shift for "
787  << "HE depth layers:";
788 #endif
789  for (i=0; i<ndx; i++) {
790  shiftHE[i] = static_cast<int>(tmp8[i]);
791 #ifdef DebugLog
792  LogDebug("HCalGeom") <<"HcalNumberingFromDDD:: gainHE[" << i << "] = "
793  << gainHE[i] << " shiftHE[" << i << "] = "
794  << shiftHE[i];
795 #endif
796  }
797 
798  //Gains and Shifts for HF depths
799  ndx = 4;
800  gainHF = getDDDArray("HFGains",sv,ndx);
801  std::vector<double> tmp9 = getDDDArray("HFShift",sv,ndx);
802  shiftHF.resize(ndx);
803 #ifdef DebugLog
804  LogDebug("HCalGeom") << "HcalNumberingFromDDD:: Gain factor and Shift for "
805  << "HF depth layers:";
806 #endif
807  for (i=0; i<ndx; i++) {
808  shiftHF[i] = static_cast<int>(tmp9[i]);
809 #ifdef DebugLog
810  LogDebug("HCalGeom") <<"HcalNumberingFromDDD:: gainHF[" << i << "] = "
811  << gainHF[i] << " shiftHF[" << i << "] = "
812  << shiftHF[i];
813 #endif
814  }
815 }
816 
817 void HcalNumberingFromDDD::loadGeometry(DDFilteredView fv) {
818 
819  bool dodet=true, hf=false;
820  std::vector<double> rb(20,0.0), ze(20,0.0), thkb(20,-1.0), thke(20,-1.0);
821  std::vector<double> zho;
822  std::vector<int> ib(20,0), ie(20,0);
823  std::vector<int> izb, phib, ize, phie, izf, phif;
824  std::vector<double> rxb;
825  rhoxb.clear(); zxb.clear(); dyxb.clear(); dzxb.clear();
826  layb.clear(); laye.clear();
827  zxe.clear(); rhoxe.clear(); dyxe.clear(); dx1e.clear(); dx2e.clear();
828  double zf = 0;
829  dzVcal = -1.;
830 
831  while (dodet) {
832  DDTranslation t = fv.translation();
833  std::vector<int> copy = fv.copyNumbers();
834  const DDSolid & sol = fv.logicalPart().solid();
835  int idet = 0, lay = -1;
836  int nsiz = (int)(copy.size());
837  if (nsiz>0) lay = copy[nsiz-1]/10;
838  if (nsiz>1) idet = copy[nsiz-2]/1000;
839  double dx=0, dy=0, dz=0, dx1=0, dx2=0;
840  if (sol.shape() == 1) {
841  const DDBox & box = static_cast<DDBox>(fv.logicalPart().solid());
842  dx = box.halfX();
843  dy = box.halfY();
844  dz = box.halfZ();
845  } else if (sol.shape() == 3) {
846  const DDTrap & trp = static_cast<DDTrap>(fv.logicalPart().solid());
847  dx1= trp.x1();
848  dx2= trp.x2();
849  dx = 0.25*(trp.x1()+trp.x2()+trp.x3()+trp.x4());
850  dy = 0.5*(trp.y1()+trp.y2());
851  dz = trp.halfZ();
852  } else if (sol.shape() == 2) {
853  const DDTubs & tub = static_cast<DDTubs>(fv.logicalPart().solid());
854  dx = tub.rIn();
855  dy = tub.rOut();
856  dz = tub.zhalf();
857  }
858  if (idet == 3) {
859  // HB
860 #ifdef DebugLog
861  LogDebug("HCalGeom") << "HB " << sol.name() << " Shape " << sol.shape()
862  << " Layer " << lay << " R " << t.Rho();
863 #endif
864  if (lay >=0 && lay < 20) {
865  ib[lay]++;
866  rb[lay] += t.Rho();
867  if (thkb[lay] <= 0) {
868  if (lay < 17) thkb[lay] = dx;
869  else thkb[lay] = std::min(dx,dy);
870  }
871  if (lay < 17) {
872  bool found = false;
873  for (unsigned int k=0; k<rxb.size(); k++) {
874  if (std::abs(rxb[k]-t.Rho()) < 0.01) {
875  found = true;
876  break;
877  }
878  }
879  if (!found) {
880  rxb.push_back(t.Rho());
881  rhoxb.push_back(t.Rho()*std::cos(t.phi()));
882  zxb.push_back(std::abs(t.z()));
883  dyxb.push_back(2.*dy);
884  dzxb.push_back(2.*dz);
885  layb.push_back(lay);
886  }
887  }
888  }
889  if (lay == 2) {
890  int iz = copy[nsiz-5];
891  int fi = copy[nsiz-4];
892  unsigned int it1 = find(iz, izb);
893  if (it1 == izb.size()) izb.push_back(iz);
894  unsigned int it2 = find(fi, phib);
895  if (it2 == phib.size()) phib.push_back(fi);
896  }
897  if (lay == 18) {
898  int ifi=-1, ich=-1;
899  if (nsiz>2) ifi = copy[nsiz-3];
900  if (nsiz>3) ich = copy[nsiz-4];
901  double z1 = std::abs((t.z()) + dz);
902  double z2 = std::abs((t.z()) - dz);
903  if (std::abs(z1-z2) < 0.01) z1 = 0;
904  if (ifi == 1 && ich == 4) {
905  if (z1 > z2) {
906  double tmp = z1;
907  z1 = z2;
908  z2 = tmp;
909  }
910  bool sok = true;
911  for (unsigned int kk=0; kk<zho.size(); kk++) {
912  if (std::abs(z2-zho[kk]) < 0.01) {
913  sok = false;
914  break;
915  } else if (z2 < zho[kk]) {
916  zho.resize(zho.size()+2);
917  for (unsigned int kz=zho.size()-1; kz>kk+1; kz=kz-2) {
918  zho[kz] = zho[kz-2];
919  zho[kz-1] = zho[kz-3];
920  }
921  zho[kk+1] = z2;
922  zho[kk] = z1;
923  sok = false;
924  break;
925  }
926  }
927  if (sok) {
928  zho.push_back(z1);
929  zho.push_back(z2);
930  }
931  LogDebug("HCalGeom") << "Detector " << idet << " Lay " << lay << " fi " << ifi << " " << ich << " z " << z1 << " " << z2;
932  }
933  }
934  } else if (idet == 4) {
935  // HE
936 #ifdef DebugLog
937  LogDebug("HCalGeom") << "HE " << sol.name() << " Shape " << sol.shape()
938  << " Layer " << lay << " Z " << t.z();
939 #endif
940  if (lay >=0 && lay < 20) {
941  ie[lay]++;
942  ze[lay] += std::abs(t.z());
943  if (thke[lay] <= 0) thke[lay] = dz;
944  bool found = false;
945  for (unsigned int k=0; k<zxe.size(); k++) {
946  if (std::abs(zxe[k]-std::abs(t.z())) < 0.01) {
947  found = true;
948  break;
949  }
950  }
951  if (!found) {
952  zxe.push_back(std::abs(t.z()));
953  rhoxe.push_back(t.Rho()*std::cos(t.phi()));
954  dyxe.push_back(dy*std::cos(t.phi()));
955  dx1 -= 0.5*(t.rho()-dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
956  dx2 -= 0.5*(t.rho()+dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
957  dx1e.push_back(-dx1);
958  dx2e.push_back(-dx2);
959  laye.push_back(lay);
960  }
961  }
962  if (copy[nsiz-1] == 21) {
963  int iz = copy[nsiz-7];
964  int fi = copy[nsiz-5];
965  unsigned int it1 = find(iz, ize);
966  if (it1 == ize.size()) ize.push_back(iz);
967  unsigned int it2 = find(fi, phie);
968  if (it2 == phie.size()) phie.push_back(fi);
969  }
970  } else if (idet == 5) {
971  // HF
972  if (!hf) {
973  const std::vector<double> & paras = sol.parameters();
974 #ifdef DebugLog
975  LogDebug("HCalGeom") << "HF " << sol.name() << " Shape " << sol.shape()
976  << " Z " << t.z() << " with " << paras.size()
977  << " Parameters";
978  for (unsigned j=0; j<paras.size(); j++)
979  LogDebug("HCalGeom") << "HF Parameter[" << j << "] = " << paras[j];
980 #endif
981  zf = fabs(t.z());
982  if (sol.shape() == ddpolycone_rrz) {
983  int nz = (int)(paras.size())-3;
984  zf += paras[3];
985  dzVcal = 0.5*(paras[nz]-paras[3]);
986  hf = true;
987  } else if (sol.shape() == ddtubs || sol.shape() == ddcons) {
988  dzVcal = paras[0];
989  zf -= paras[0];
990  hf = true;
991  }
992  }
993 #ifdef DebugLog
994  } else {
995  LogDebug("HCalGeom") << "Unknown Detector " << idet << " for "
996  << sol.name() << " Shape " << sol.shape() << " R "
997  << t.Rho() << " Z " << t.z();
998 #endif
999  }
1000  dodet = fv.next();
1001  }
1002 
1003  int ibmx = 0, iemx = 0;
1004  for (int i = 0; i < 20; i++) {
1005  if (ib[i]>0) {
1006  rb[i] /= (double)(ib[i]);
1007  ibmx = i+1;
1008  }
1009  if (ie[i]>0) {
1010  ze[i] /= (double)(ie[i]);
1011  iemx = i+1;
1012  }
1013 #ifdef DebugLog
1014  LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
1015  << rb[i] << " Endcap " << ie[i] << " " << ze[i];
1016 #endif
1017  }
1018  for (int i = 4; i >= 0; i--) {
1019  if (ib[i] == 0) {rb[i] = rb[i+1]; thkb[i] = thkb[i+1];}
1020  if (ie[i] == 0) {ze[i] = ze[i+1]; thke[i] = thke[i+1];}
1021 #ifdef DebugLog
1022  if (ib[i] == 0 || ie[i] == 0)
1023  LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
1024  << rb[i] << " Endcap " << ie[i] << " " << ze[i];
1025 #endif
1026  }
1027 
1028 #ifdef DebugLog
1029  for (unsigned int k=0; k<layb.size(); ++k)
1030  std::cout << "HB: " << layb[k] << " R " << rxb[k] << " " << rhoxb[k] << " Z " << zxb[k] << " DY " << dyxb[k] << " DZ " << dzxb[k] << "\n";
1031  for (unsigned int k=0; k<laye.size(); ++k)
1032  std::cout << "HE: " << laye[k] << " R " << rhoxe[k] << " Z " << zxe[k] << " X1|X2 " << dx1e[k] << "|" << dx2e[k] << " DY " << dyxe[k] << "\n";
1033 
1034  printTile();
1035  LogDebug("HCalGeom") << "HcalNumberingFromDDD: Maximum Layer for HB "
1036  << ibmx << " for HE " << iemx << " Z for HF " << zf
1037  << " extent " << dzVcal;
1038 #endif
1039 
1040  if (ibmx > 0) {
1041  rHB.resize(ibmx);
1042  drHB.resize(ibmx);
1043  for (int i=0; i<ibmx; i++) {
1044  rHB[i] = rb[i];
1045  drHB[i] = thkb[i];
1046 #ifdef DebugLog
1047  LogDebug("HCalGeom") << "HcalNumberingFromDDD: rHB[" << i << "] = "
1048  << rHB[i] << " drHB[" << i << "] = " << drHB[i];
1049 #endif
1050  }
1051  }
1052  if (iemx > 0) {
1053  zHE.resize(iemx);
1054  dzHE.resize(iemx);
1055  for (int i=0; i<iemx; i++) {
1056  zHE[i] = ze[i];
1057  dzHE[i] = thke[i];
1058 #ifdef DebugLog
1059  LogDebug("HCalGeom") << "HcalNumberingFromDDD: zHE[" << i << "] = "
1060  << zHE[i] << " dzHE[" << i << "] = " << dzHE[i];
1061 #endif
1062  }
1063  }
1064 
1065  nzHB = (int)(izb.size());
1066  nmodHB = (int)(phib.size());
1067 #ifdef DebugLog
1068  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHB
1069  << " barrel half-sectors";
1070  for (int i=0; i<nzHB; i++)
1071  LogDebug("HCalGeom") << "Section " << i << " Copy number " << izb[i];
1072  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHB
1073  << " barrel modules";
1074  for (int i=0; i<nmodHB; i++)
1075  LogDebug("HCalGeom") << "Module " << i << " Copy number " << phib[i];
1076 #endif
1077 
1078  nzHE = (int)(ize.size());
1079  nmodHE = (int)(phie.size());
1080 #ifdef DebugLog
1081  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHE
1082  << " endcap half-sectors";
1083  for (int i=0; i<nzHE; i++)
1084  LogDebug("HCalGeom") << "Section " << i << " Copy number " << ize[i];
1085  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHE
1086  << " endcap modules";
1087  for (int i=0; i<nmodHE; i++)
1088  LogDebug("HCalGeom") << "Module " << i << " Copy number " << phie[i];
1089 #endif
1090 
1091  LogDebug("HCalGeom") << "HO has Z of size " << zho.size();
1092  for (unsigned int kk=0; kk<zho.size(); kk++)
1093  LogDebug("HCalGeom") << "ZHO[" << kk << "] = " << zho[kk];
1094  if (ibmx > 17 && zho.size() > 2) {
1095  etaHO[0] = getEta(0.5*(rHB[17]+rHB[18]), zho[1]);
1096  etaHO[1] = getEta(rHB[18]+drHB[18], zho[2]);
1097  etaHO[2] = getEta(rHB[18]-drHB[18], zho[3]);
1098  etaHO[3] = getEta(rHB[18]+drHB[18], zho[4]);
1099  } else {
1100  etaHO[0] = etaTable[4];
1101  etaHO[1] = etaTable[4];
1102  etaHO[2] = etaTable[10];
1103  etaHO[3] = etaTable[10];
1104  }
1105  LogDebug("HCalGeom") << "HO Eta boundaries " << etaHO[0] << " " << etaHO[1]
1106  << " " << etaHO[2] << " " << etaHO[3];
1107 }
1108 
1109 std::vector<double> HcalNumberingFromDDD::getDDDArray(const std::string & str,
1110  const DDsvalues_type & sv,
1111  int & nmin) const {
1112 #ifdef DebugLog
1113  LogDebug("HCalGeom") << "HcalNumberingFromDDD:getDDDArray called for "
1114  << str << " with nMin " << nmin;
1115 #endif
1116  DDValue value(str);
1117  if (DDfetch(&sv,value)) {
1118 #ifdef DebugLog
1119  LogDebug("HCalGeom") << "HcalNumberingFromDDD: " << value;
1120 #endif
1121  const std::vector<double> & fvec = value.doubles();
1122  int nval = fvec.size();
1123  if (nmin > 0) {
1124  if (nval < nmin) {
1125  edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str
1126  << " bins " << nval << " < " << nmin
1127  << " ==> illegal";
1128  throw DDException("HcalNumberingFromDDD: cannot get array "+str);
1129  }
1130  } else {
1131  if (nval < 2) {
1132  edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str
1133  << " bins " << nval << " < 2 ==> illegal"
1134  << " (nmin=" << nmin << ")";
1135  throw DDException("HcalNumberingFromDDD: cannot get array "+str);
1136  }
1137  }
1138  nmin = nval;
1139  return fvec;
1140  } else {
1141  edm::LogError("HCalGeom") << "HcalNumberingFromDDD: cannot get array "
1142  << str;
1143  throw DDException("HcalNumberingFromDDD: cannot get array "+str);
1144  }
1145 }
1146 
1147 int HcalNumberingFromDDD::getShift(HcalSubdetector subdet, int depth) const {
1148 
1149  int shift;
1150  switch(subdet) {
1151  case HcalEndcap:
1152  shift = shiftHE[depth-1];
1153  break;
1154  case HcalForward:
1155  shift = shiftHF[depth-1];
1156  break;
1157  default:
1158  shift = shiftHB[depth-1];
1159  break;
1160  }
1161  return shift;
1162 }
1163 
1164 double HcalNumberingFromDDD::getGain(HcalSubdetector subdet, int depth) const {
1165 
1166  double gain;
1167  switch(subdet) {
1168  case HcalEndcap:
1169  gain = gainHE[depth-1];
1170  break;
1171  case HcalForward:
1172  gain = gainHF[depth-1];
1173  break;
1174  default:
1175  gain = gainHB[depth-1];
1176  break;
1177  }
1178  return gain;
1179 }
1180 
1181 unsigned int HcalNumberingFromDDD::find(int element,
1182  std::vector<int> array) const {
1183 
1184  unsigned int id = array.size();
1185  for (unsigned int i = 0; i < array.size(); i++) {
1186  if (element == array[i]) {
1187  id = i;
1188  break;
1189  }
1190  }
1191  return id;
1192 }
1193 
1194 int HcalNumberingFromDDD::unitPhi(int det, int etaR) const {
1195 
1196  const double fiveDegInRad = 2*M_PI/72;
1197  int units=0;
1198  if (det == static_cast<int>(HcalForward))
1199  units=int(phibin[nEta+etaR-etaMin[2]-1]/fiveDegInRad+0.5);
1200  else
1201  units=int(phibin[etaR-1]/fiveDegInRad+0.5);
1202 
1203  return units;
1204 }
1205 
1206 void HcalNumberingFromDDD::tileHB(int eta, int depth) {
1207 
1208  double etaL = etaTable[eta-1];
1209  double thetaL = 2.*atan(exp(-etaL));
1210  double etaH = etaTable[eta];
1211  double thetaH = 2.*atan(exp(-etaH));
1212  int layL=0, layH=0;
1213  if (depth == 1) {
1214  layH = depth1[eta-1];
1215  } else {
1216  layL = depth1[eta-1];
1217  layH = depth2[eta-1];
1218  }
1219  std::cout << "\ntileHB:: eta|depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << "\n";
1220  for (int lay=layL; lay<layH; ++lay) {
1221  std::vector<double> area(2,0);
1222  int kk=0;
1223  for (unsigned int k=0; k<layb.size(); ++k) {
1224  if (lay == layb[k]) {
1225  double zmin = rhoxb[k]*std::cos(thetaL)/std::sin(thetaL);
1226  double zmax = rhoxb[k]*std::cos(thetaH)/std::sin(thetaH);
1227  double dz = (std::min(zmax,dzxb[k]) - zmin);
1228  if (dz > 0) {
1229  area[kk] = dz*dyxb[k];
1230  kk++;
1231  }
1232  }
1233  }
1234  if (area[0] > 0) std::cout << std::setw(2) << lay << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
1235  }
1236 }
1237 
1238 void HcalNumberingFromDDD::tileHE(int eta, int depth) {
1239 
1240  double etaL = etaTable[eta-1];
1241  double thetaL = 2.*atan(exp(-etaL));
1242  double etaH = etaTable[eta];
1243  double thetaH = 2.*atan(exp(-etaH));
1244  int layL=0, layH=0;
1245  if (eta == 16) {
1246  layH = depth3[eta-1];
1247  } else if (depth == 1) {
1248  layH = depth1[eta-1];
1249  } else if (depth == 2) {
1250  layL = depth1[eta-1];
1251  layH = depth2[eta-1];
1252  } else {
1253  layL = depth2[eta-1];
1254  layH = depth3[eta-1];
1255  }
1256  double phib = phibin[eta-1];
1257  int nphi = 2;
1258  if (phib > 6*CLHEP::deg) nphi = 1;
1259  std::cout << "\ntileHE:: Eta/depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << " phi " << nphi << "\n";
1260  for (int lay=layL; lay<layH; ++lay) {
1261  std::vector<double> area(4,0);
1262  int kk=0;
1263  for (unsigned int k=0; k<laye.size(); ++k) {
1264  if (lay == laye[k]) {
1265  double rmin = zxe[k]*std::tan(thetaH);
1266  double rmax = zxe[k]*std::tan(thetaL);
1267  if ((lay != 0 || eta == 18) &&
1268  (lay != 1 || (eta == 18 && rhoxe[k]-dyxe[k] > 1000) || (eta != 18 && rhoxe[k]-dyxe[k] < 1000)) &&
1269  rmin+30 < rhoxe[k]+dyxe[k] && rmax > rhoxe[k]-dyxe[k]) {
1270  rmin = std::max(rmin,rhoxe[k]-dyxe[k]);
1271  rmax = std::min(rmax,rhoxe[k]+dyxe[k]);
1272  double dx1 = rmin*std::tan(phib);
1273  double dx2 = rmax*std::tan(phib);
1274  double ar1=0, ar2=0;
1275  if (nphi == 1) {
1276  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-4.*dx1e[k]);
1277  } else {
1278  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-2.*dx1e[k]);
1279  ar2 = 0.5*(rmax-rmin)*((rmax+rmin)*tan(10.*CLHEP::deg)-4*dx1e[k])-ar1;
1280  }
1281  area[kk] = ar1;
1282  area[kk+2] = ar2;
1283  kk++;
1284  }
1285  }
1286  }
1287  if (area[0] > 0 && area[1] > 0) {
1288  int lay0 = lay-1;
1289  if (eta == 18) lay0++;
1290  if (nphi == 1) {
1291  std::cout << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
1292  } else {
1293  std::cout << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << ":" << std::setw(8) << area[2] << " " << std::setw(8) << area[3] << "\n";
1294  }
1295  }
1296  }
1297 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:510
DDTrap::halfZ
double halfZ() const
half of the z-Axis
Definition: DDSolid.cc:170
HcalNumberingFromDDD::shiftHE
std::vector< int > shiftHE
Definition: HcalNumberingFromDDD.h:82
DDSolid::parameters
const std::vector< double > & parameters() const
Don&#39;t use (only meant to be used by DDbox(), DDtub(), ...)
Definition: DDSolid.cc:153
HcalNumberingFromDDD::shiftHB
std::vector< int > shiftHB
Definition: HcalNumberingFromDDD.h:80
HcalNumberingFromDDD::dx2e
std::vector< double > dx2e
Definition: HcalNumberingFromDDD.h:96
DDBox::halfZ
double halfZ() const
Definition: DDSolid.cc:261
DDTrap::y2
double y2() const
Half-length along y of the face at +pDz.
Definition: DDSolid.cc:184
i
int i
Definition: DBlmapReader.cc:9
DDFilteredView::logicalPart
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
Definition: DDFilteredView.cc:18
DDValue::doubles
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:151
DDFilteredView::addFilter
void addFilter(const DDFilter &, log_op op=AND)
Definition: DDFilteredView.cc:24
DDBase::name
const N & name() const
Definition: DDBase.h:88
HcalNumberingFromDDD::depth2
std::vector< int > depth2
Definition: HcalNumberingFromDDD.h:76
DDTrap::x2
double x2() const
Half-length along x of the side at y=+pDy1 of the face at -pDz.
Definition: DDSolid.cc:180
DDSolid::shape
DDSolidShape shape() const
The type of the solid.
Definition: DDSolid.cc:147
DDFilteredView::copyNumbers
nav_type copyNumbers() const
return the stack of copy numbers
Definition: DDFilteredView.cc:282
HcalNumberingFromDDD::dyxe
std::vector< double > dyxe
Definition: HcalNumberingFromDDD.h:96
DDTubs::rOut
double rOut() const
Definition: DDSolid.cc:512
HcalNumberingFromDDD::tileHB
void tileHB(int eta, int depth)
Definition: HcalNumberingFromDDD.cc:1206
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
DDException
An exception for DDD errors.
Definition: DDException.h:23
HcalCellType::setMissingPhi
void setMissingPhi(std::vector< int >, std::vector< int >)
Definition: HcalCellType.cc:57
HcalNumberingFromDDD::zxe
std::vector< double > zxe
Definition: HcalNumberingFromDDD.h:96
HcalNumberingFromDDD::etaTable
std::vector< double > etaTable
Definition: HcalNumberingFromDDD.h:67
HcalNumberingFromDDD::layb
std::vector< int > layb
Definition: HcalNumberingFromDDD.h:95
abs
#define abs(x)
Definition: mlp_lapack.h:159
HcalNumberingFromDDD::getShift
int getShift(HcalSubdetector subdet, int depth) const
Definition: HcalNumberingFromDDD.cc:1147
HcalNumberingFromDDD::zxb
std::vector< double > zxb
Definition: HcalNumberingFromDDD.h:94
theta
Geom::Theta< T > theta() const
Definition: newBasic3DVector.h:115
min
#define min(a, b)
Definition: mlp_lapack.h:161
HcalNumberingFromDDD::unitPhi
int unitPhi(int det, int etaR) const
Definition: HcalNumberingFromDDD.cc:1194
HcalNumberingFromDDD::rHB
std::vector< double > rHB
Definition: HcalNumberingFromDDD.h:89
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
HcalNumberingFromDDD::HcalCellTypes
std::vector< HcalCellType > HcalCellTypes() const
Definition: HcalNumberingFromDDD.cc:391
HcalNumberingFromDDD::drHB
std::vector< double > drHB
Definition: HcalNumberingFromDDD.h:89
HcalNumberingFromDDD::getGain
double getGain(HcalSubdetector subdet, int depth) const
Definition: HcalNumberingFromDDD.cc:1164
debug_cff.d1
tuple d1
Definition: debug_cff.py:7
DDCompactView
type of data representation of DDCompactView
Definition: DDCompactView.h:81
HcalNumberingFromDDD::phibin
std::vector< double > phibin
Definition: HcalNumberingFromDDD.h:73
DDfetch
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:102
eta
T eta() const
Definition: newBasic3DVector.h:122
HcalNumberingFromDDD::rhoxe
std::vector< double > rhoxe
Definition: HcalNumberingFromDDD.h:96
DDSolid
A DDSolid represents the shape of a part.
Definition: DDSolid.h:42
HcalNumberingFromDDD::numberOfCells
unsigned int numberOfCells(HcalSubdetector) const
Definition: HcalNumberingFromDDD.cc:373
AlCaRecoCosmics_cfg.name
tuple name
Definition: AlCaRecoCosmics_cfg.py:24
DDTranslation
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > DDTranslation
Definition: DDTranslation.h:7
HcalNumberingFromDDD::HcalNumberingFromDDD
HcalNumberingFromDDD(std::string &name, const DDCompactView &cpv)
Definition: HcalNumberingFromDDD.cc:22
HcalNumberingFromDDD::dzxb
std::vector< double > dzxb
Definition: HcalNumberingFromDDD.h:94
HcalNumberingFromDDD::gainHF
std::vector< double > gainHF
Definition: HcalNumberingFromDDD.h:83
HcalNumberingFromDDD::shiftHF
std::vector< int > shiftHF
Definition: HcalNumberingFromDDD.h:84
DDBox::halfX
double halfX() const
Definition: DDSolid.cc:257
HcalNumberingFromDDD::dyxb
std::vector< double > dyxb
Definition: HcalNumberingFromDDD.h:94
z
Definition: DDAxes.h:10
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
DDFilteredView::next
bool next()
set current node to the next node in the filtered tree
Definition: DDFilteredView.cc:104
DDBox::halfY
double halfY() const
Definition: DDSolid.cc:259
HcalNumberingFromDDD::phioff
std::vector< double > phioff
Definition: HcalNumberingFromDDD.h:66
HcalNumberingFromDDD::dx1e
std::vector< double > dx1e
Definition: HcalNumberingFromDDD.h:96
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
HcalNumberingFromDDD::initialize
void initialize(std::string &name, const DDCompactView &cpv)
Definition: HcalNumberingFromDDD.cc:609
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
HcalNumberingFromDDD::laye
std::vector< int > laye
Definition: HcalNumberingFromDDD.h:95
DDsvalues_type
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
std::maps an index to a DDValue. The index corresponds to the index assigned to the name of the std::...
Definition: DDsvalues.h:19
dmin
#define dmin(a, b)
Definition: mlp_lapack.h:163
DDTrap
Interface to a Trapezoid.
Definition: DDSolid.h:115
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
HcalSubdetector
HcalSubdetector
Definition: HcalAssistant.h:32
j
int j
Definition: DBlmapReader.cc:9
HcalNumberingFromDDD::nOff
std::vector< int > nOff
Definition: HcalNumberingFromDDD.h:88
ExpressReco_HICollisions_FallBack.id
tuple id
Definition: ExpressReco_HICollisions_FallBack.py:9307
HcalNumberingFromDDD::depth3
std::vector< int > depth3
Definition: HcalNumberingFromDDD.h:77
corr
JetCorrectorParameters corr
Definition: classes.h:9
DDTubs::rIn
double rIn() const
Definition: DDSolid.cc:510
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
HcalNumberingFromDDD::rhoxb
std::vector< double > rhoxb
Definition: HcalNumberingFromDDD.h:94
HcalNumberingFromDDD::find
unsigned find(int element, std::vector< int > array) const
Definition: HcalNumberingFromDDD.cc:1181
DDTrap::y1
double y1() const
Half-length along y of the face at -pDz.
Definition: DDSolid.cc:176
HcalNumberingFromDDD::etaMin
std::vector< int > etaMin
Definition: HcalNumberingFromDDD.h:71
HcalNumberingFromDDD::etaMax
std::vector< int > etaMax
Definition: HcalNumberingFromDDD.h:72
DDTrap::x4
double x4() const
Half-length along x of the side at y=+pDy2 of the face at +pDz.
Definition: DDSolid.cc:188
M_PI
#define M_PI
Definition: BFit3D.cc:3
HcalNumberingFromDDD::loadGeometry
void loadGeometry(DDFilteredView)
Definition: HcalNumberingFromDDD.cc:817
funct::log
Log< T >::type log(const T &t)
Definition: Log.h:22
HcalNumberingFromDDD::deltaEta
double deltaEta(int det, int eta, int depth) const
Definition: HcalNumberingFromDDD.cc:571
align_tpl.filter
tuple filter
USE THIS FOR SKIMMED TRACKS process.p = cms.Path(process.hltLevel1GTSeed*process.skimming*process.offlineBeamSpot*process.TrackRefitter2) OTHERWISE USE THIS.
Definition: align_tpl.py:86
HcalNumberingFromDDD::getDDDArray
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &) const
Definition: HcalNumberingFromDDD.cc:1109
DDBox
Definition: DDSolid.h:185
HcalNumberingFromDDD::getEta
double getEta(int det, int etaR, int zside, int depth=1) const
Definition: HcalNumberingFromDDD.cc:520
HcalNumberingFromDDD::gainHB
std::vector< double > gainHB
Definition: HcalNumberingFromDDD.h:79
evf::utils::num
long long int num
Definition: procUtils.cc:71
HcalNumberingFromDDD::gainHE
std::vector< double > gainHE
Definition: HcalNumberingFromDDD.h:81
HcalNumberingFromDDD::depth1
std::vector< int > depth1
Definition: HcalNumberingFromDDD.h:75
DDFilteredView::mergedSpecifics
DDsvalues_type mergedSpecifics() const
Definition: DDFilteredView.cc:66
HcalNumberingFromDDD::zHE
std::vector< double > zHE
Definition: HcalNumberingFromDDD.h:90
HcalNumberingFromDDD::tileHE
void tileHE(int eta, int depth)
Definition: HcalNumberingFromDDD.cc:1238
DDTubs::zhalf
double zhalf() const
Definition: DDSolid.cc:508
HcalNumberingFromDDD::dzHE
std::vector< double > dzHE
Definition: HcalNumberingFromDDD.h:90
DDTrap::x3
double x3() const
Half-length along x of the side at y=-pDy2 of the face at +pDz.
Definition: DDSolid.cc:186
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
DDLogicalPart::solid
const DDSolid & solid() const
Returns a reference object of the solid being the shape of this LogicalPart.
Definition: DDLogicalPart.cc:148
DDFilteredView::firstChild
bool firstChild()
set the current node to the first child ...
Definition: DDFilteredView.cc:131
HcalNumberingFromDDD::loadSpecPars
void loadSpecPars(DDFilteredView)
Definition: HcalNumberingFromDDD.cc:645
edm::shift
static unsigned int const shift
Definition: LuminosityBlockID.cc:9
DDSpecificsFilter::setCriteria
void setCriteria(const DDValue &nameVal, comp_op, log_op l=AND, bool asString=true, bool merged=true)
Definition: DDFilter.cc:285
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
HcalNumberingFromDDD::getEtaTable
std::vector< double > getEtaTable() const
Definition: HcalNumberingFromDDD.cc:367
dmax
#define dmax(a, b)
Definition: mlp_lapack.h:164
HcalNumberingFromDDD::cell
HcalCellType::HcalCell cell(int det, int zside, int depth, int etaR, int iphi, bool corr=true) const
Definition: HcalNumberingFromDDD.cc:210
DDFilteredView::translation
const DDTranslation & translation() const
The absolute translation of the current node.
Definition: DDFilteredView.cc:32
HcalNumberingFromDDD::unitID
HcalID unitID(int det, CLHEP::Hep3Vector pos, int depth, int lay=-1) const
Definition: HcalNumberingFromDDD.cc:32
DDTrap::x1
double x1() const
Half-length along x of the side at y=-pDy1 of the face at -pDz.
Definition: DDSolid.cc:178
point
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
HcalNumberingFromDDD::rTable
std::vector< double > rTable
Definition: HcalNumberingFromDDD.h:69
phi
Definition: DDAxes.h:10
DDSpecificsFilter
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.
Definition: DDFilter.h:37
Generated for CMSSW Reference Manual by doxygen 1.8.5