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/Core/interface/DDFilter.h"
9 #include "DetectorDescription/Core/interface/DDFilteredView.h"
10 #include "DetectorDescription/Core/interface/DDSolid.h"
11 #include "DetectorDescription/Core/interface/DDValue.h"
12 
13 #include "FWCore/MessageLogger/interface/MessageLogger.h"
14 
15 #include "CLHEP/Units/GlobalPhysicalConstants.h"
16 #include "CLHEP/Units/GlobalSystemOfUnits.h"
17 #include <iostream>
18 
19 //#define DebugLog
20 
21 HcalNumberingFromDDD::HcalNumberingFromDDD(std::string & name,
22  const DDCompactView & cpv) {
23  edm::LogInfo("HCalGeom") << "Creating HcalNumberingFromDDD";
24  initialize(name, cpv);
25 }
26 
27 HcalNumberingFromDDD::~HcalNumberingFromDDD() {
28  edm::LogInfo("HCalGeom") << "Deleting HcalNumberingFromDDD";
29 }
30 
31 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
32  CLHEP::Hep3Vector point,
33  int depth,
34  int lay) const {
35 
36 
37  double hx = point.x();
38  double hy = point.y();
39  double hz = point.z();
40  double hR = sqrt(hx*hx+hy*hy+hz*hz);
41  double htheta = (hR == 0. ? 0. : acos(std::max(std::min(hz/hR,1.0),-1.0)));
42  double hsintheta = sin(htheta);
43  double hphi = (hR*hsintheta == 0. ? 0. :atan2(hy,hx));
44  double heta = (fabs(hsintheta) == 1.? 0. : -log(fabs(tan(htheta/2.))) );
45 
46  int hsubdet=0;
47  double etaR;
48 
49  //First eta index
50  if (det == 5) { // Forward HCal
51  hsubdet = static_cast<int>(HcalForward);
52  hR = sqrt(hx*hx+hy*hy);
53  etaR = (heta >= 0. ? hR : -hR);
54  } else { // Barrel or Endcap
55  etaR = heta;
56  if (det == 3) {
57  hsubdet = static_cast<int>(HcalBarrel);
58  if (zho.size() > 4) etaR = getEtaHO(heta,hx,hy,hz);
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 cms::Exception("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<int> ib(20,0), ie(20,0);
822  std::vector<int> izb, phib, ize, phie, izf, phif;
823  std::vector<double> rxb;
824  rhoxb.clear(); zxb.clear(); dyxb.clear(); dzxb.clear();
825  layb.clear(); laye.clear();
826  zxe.clear(); rhoxe.clear(); dyxe.clear(); dx1e.clear(); dx2e.clear();
827  double zf = 0;
828  dzVcal = -1.;
829 
830  while (dodet) {
831  DDTranslation t = fv.translation();
832  std::vector<int> copy = fv.copyNumbers();
833  const DDSolid & sol = fv.logicalPart().solid();
834  int idet = 0, lay = -1;
835  int nsiz = (int)(copy.size());
836  if (nsiz>0) lay = copy[nsiz-1]/10;
837  if (nsiz>1) idet = copy[nsiz-2]/1000;
838  double dx=0, dy=0, dz=0, dx1=0, dx2=0;
839  if (sol.shape() == 1) {
840  const DDBox & box = static_cast<DDBox>(fv.logicalPart().solid());
841  dx = box.halfX();
842  dy = box.halfY();
843  dz = box.halfZ();
844  } else if (sol.shape() == 3) {
845  const DDTrap & trp = static_cast<DDTrap>(fv.logicalPart().solid());
846  dx1= trp.x1();
847  dx2= trp.x2();
848  dx = 0.25*(trp.x1()+trp.x2()+trp.x3()+trp.x4());
849  dy = 0.5*(trp.y1()+trp.y2());
850  dz = trp.halfZ();
851  } else if (sol.shape() == 2) {
852  const DDTubs & tub = static_cast<DDTubs>(fv.logicalPart().solid());
853  dx = tub.rIn();
854  dy = tub.rOut();
855  dz = tub.zhalf();
856  }
857  if (idet == 3) {
858  // HB
859 #ifdef DebugLog
860  LogDebug("HCalGeom") << "HB " << sol.name() << " Shape " << sol.shape()
861  << " Layer " << lay << " R " << t.Rho();
862 #endif
863  if (lay >=0 && lay < 20) {
864  ib[lay]++;
865  rb[lay] += t.Rho();
866  if (thkb[lay] <= 0) {
867  if (lay < 17) thkb[lay] = dx;
868  else thkb[lay] = std::min(dx,dy);
869  }
870  if (lay < 17) {
871  bool found = false;
872  for (unsigned int k=0; k<rxb.size(); k++) {
873  if (std::abs(rxb[k]-t.Rho()) < 0.01) {
874  found = true;
875  break;
876  }
877  }
878  if (!found) {
879  rxb.push_back(t.Rho());
880  rhoxb.push_back(t.Rho()*std::cos(t.phi()));
881  zxb.push_back(std::abs(t.z()));
882  dyxb.push_back(2.*dy);
883  dzxb.push_back(2.*dz);
884  layb.push_back(lay);
885  }
886  }
887  }
888  if (lay == 2) {
889  int iz = copy[nsiz-5];
890  int fi = copy[nsiz-4];
891  unsigned int it1 = find(iz, izb);
892  if (it1 == izb.size()) izb.push_back(iz);
893  unsigned int it2 = find(fi, phib);
894  if (it2 == phib.size()) phib.push_back(fi);
895  }
896  if (lay == 18) {
897  int ifi=-1, ich=-1;
898  if (nsiz>2) ifi = copy[nsiz-3];
899  if (nsiz>3) ich = copy[nsiz-4];
900  double z1 = std::abs((t.z()) + dz);
901  double z2 = std::abs((t.z()) - dz);
902  if (std::abs(z1-z2) < 0.01) z1 = 0;
903  if (ifi == 1 && ich == 4) {
904  if (z1 > z2) {
905  double tmp = z1;
906  z1 = z2;
907  z2 = tmp;
908  }
909  bool sok = true;
910  for (unsigned int kk=0; kk<zho.size(); kk++) {
911  if (std::abs(z2-zho[kk]) < 0.01) {
912  sok = false;
913  break;
914  } else if (z2 < zho[kk]) {
915  zho.resize(zho.size()+2);
916  for (unsigned int kz=zho.size()-1; kz>kk+1; kz=kz-2) {
917  zho[kz] = zho[kz-2];
918  zho[kz-1] = zho[kz-3];
919  }
920  zho[kk+1] = z2;
921  zho[kk] = z1;
922  sok = false;
923  break;
924  }
925  }
926  if (sok) {
927  zho.push_back(z1);
928  zho.push_back(z2);
929  }
930 #ifdef DebugLog
931  LogDebug("HCalGeom") << "Detector " << idet << " Lay " << lay << " fi " << ifi << " " << ich << " z " << z1 << " " << z2;
932 #endif
933  }
934  }
935  } else if (idet == 4) {
936  // HE
937 #ifdef DebugLog
938  LogDebug("HCalGeom") << "HE " << sol.name() << " Shape " << sol.shape()
939  << " Layer " << lay << " Z " << t.z();
940 #endif
941  if (lay >=0 && lay < 20) {
942  ie[lay]++;
943  ze[lay] += std::abs(t.z());
944  if (thke[lay] <= 0) thke[lay] = dz;
945  bool found = false;
946  for (unsigned int k=0; k<zxe.size(); k++) {
947  if (std::abs(zxe[k]-std::abs(t.z())) < 0.01) {
948  found = true;
949  break;
950  }
951  }
952  if (!found) {
953  zxe.push_back(std::abs(t.z()));
954  rhoxe.push_back(t.Rho()*std::cos(t.phi()));
955  dyxe.push_back(dy*std::cos(t.phi()));
956  dx1 -= 0.5*(t.rho()-dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
957  dx2 -= 0.5*(t.rho()+dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
958  dx1e.push_back(-dx1);
959  dx2e.push_back(-dx2);
960  laye.push_back(lay);
961  }
962  }
963  if (copy[nsiz-1] == 21) {
964  int iz = copy[nsiz-7];
965  int fi = copy[nsiz-5];
966  unsigned int it1 = find(iz, ize);
967  if (it1 == ize.size()) ize.push_back(iz);
968  unsigned int it2 = find(fi, phie);
969  if (it2 == phie.size()) phie.push_back(fi);
970  }
971  } else if (idet == 5) {
972  // HF
973  if (!hf) {
974  const std::vector<double> & paras = sol.parameters();
975 #ifdef DebugLog
976  LogDebug("HCalGeom") << "HF " << sol.name() << " Shape " << sol.shape()
977  << " Z " << t.z() << " with " << paras.size()
978  << " Parameters";
979  for (unsigned j=0; j<paras.size(); j++)
980  LogDebug("HCalGeom") << "HF Parameter[" << j << "] = " << paras[j];
981 #endif
982  zf = fabs(t.z());
983  if (sol.shape() == ddpolycone_rrz) {
984  int nz = (int)(paras.size())-3;
985  zf += paras[3];
986  dzVcal = 0.5*(paras[nz]-paras[3]);
987  hf = true;
988  } else if (sol.shape() == ddtubs || sol.shape() == ddcons) {
989  dzVcal = paras[0];
990  zf -= paras[0];
991  hf = true;
992  }
993  }
994 #ifdef DebugLog
995  } else {
996  LogDebug("HCalGeom") << "Unknown Detector " << idet << " for "
997  << sol.name() << " Shape " << sol.shape() << " R "
998  << t.Rho() << " Z " << t.z();
999 #endif
1000  }
1001  dodet = fv.next();
1002  }
1003 
1004  int ibmx = 0, iemx = 0;
1005  for (int i = 0; i < 20; i++) {
1006  if (ib[i]>0) {
1007  rb[i] /= (double)(ib[i]);
1008  ibmx = i+1;
1009  }
1010  if (ie[i]>0) {
1011  ze[i] /= (double)(ie[i]);
1012  iemx = i+1;
1013  }
1014 #ifdef DebugLog
1015  LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
1016  << rb[i] << " Endcap " << ie[i] << " " << ze[i];
1017 #endif
1018  }
1019  for (int i = 4; i >= 0; i--) {
1020  if (ib[i] == 0) {rb[i] = rb[i+1]; thkb[i] = thkb[i+1];}
1021  if (ie[i] == 0) {ze[i] = ze[i+1]; thke[i] = thke[i+1];}
1022 #ifdef DebugLog
1023  if (ib[i] == 0 || ie[i] == 0)
1024  LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
1025  << rb[i] << " Endcap " << ie[i] << " " << ze[i];
1026 #endif
1027  }
1028 
1029 #ifdef DebugLog
1030  for (unsigned int k=0; k<layb.size(); ++k)
1031  std::cout << "HB: " << layb[k] << " R " << rxb[k] << " " << rhoxb[k] << " Z " << zxb[k] << " DY " << dyxb[k] << " DZ " << dzxb[k] << "\n";
1032  for (unsigned int k=0; k<laye.size(); ++k)
1033  std::cout << "HE: " << laye[k] << " R " << rhoxe[k] << " Z " << zxe[k] << " X1|X2 " << dx1e[k] << "|" << dx2e[k] << " DY " << dyxe[k] << "\n";
1034 
1035  printTile();
1036  LogDebug("HCalGeom") << "HcalNumberingFromDDD: Maximum Layer for HB "
1037  << ibmx << " for HE " << iemx << " Z for HF " << zf
1038  << " extent " << dzVcal;
1039 #endif
1040 
1041  if (ibmx > 0) {
1042  rHB.resize(ibmx);
1043  drHB.resize(ibmx);
1044  for (int i=0; i<ibmx; i++) {
1045  rHB[i] = rb[i];
1046  drHB[i] = thkb[i];
1047 #ifdef DebugLog
1048  LogDebug("HCalGeom") << "HcalNumberingFromDDD: rHB[" << i << "] = "
1049  << rHB[i] << " drHB[" << i << "] = " << drHB[i];
1050 #endif
1051  }
1052  }
1053  if (iemx > 0) {
1054  zHE.resize(iemx);
1055  dzHE.resize(iemx);
1056  for (int i=0; i<iemx; i++) {
1057  zHE[i] = ze[i];
1058  dzHE[i] = thke[i];
1059 #ifdef DebugLog
1060  LogDebug("HCalGeom") << "HcalNumberingFromDDD: zHE[" << i << "] = "
1061  << zHE[i] << " dzHE[" << i << "] = " << dzHE[i];
1062 #endif
1063  }
1064  }
1065 
1066  nzHB = (int)(izb.size());
1067  nmodHB = (int)(phib.size());
1068 #ifdef DebugLog
1069  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHB
1070  << " barrel half-sectors";
1071  for (int i=0; i<nzHB; i++)
1072  LogDebug("HCalGeom") << "Section " << i << " Copy number " << izb[i];
1073  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHB
1074  << " barrel modules";
1075  for (int i=0; i<nmodHB; i++)
1076  LogDebug("HCalGeom") << "Module " << i << " Copy number " << phib[i];
1077 #endif
1078 
1079  nzHE = (int)(ize.size());
1080  nmodHE = (int)(phie.size());
1081 #ifdef DebugLog
1082  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHE
1083  << " endcap half-sectors";
1084  for (int i=0; i<nzHE; i++)
1085  LogDebug("HCalGeom") << "Section " << i << " Copy number " << ize[i];
1086  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHE
1087  << " endcap modules";
1088  for (int i=0; i<nmodHE; i++)
1089  LogDebug("HCalGeom") << "Module " << i << " Copy number " << phie[i];
1090 #endif
1091 
1092 #ifdef DebugLog
1093  LogDebug("HCalGeom") << "HO has Z of size " << zho.size();
1094  for (unsigned int kk=0; kk<zho.size(); kk++)
1095  LogDebug("HCalGeom") << "ZHO[" << kk << "] = " << zho[kk];
1096 #endif
1097  if (ibmx > 17 && zho.size() > 4) {
1098  rminHO = rHB[17]-100.0;
1099  etaHO[0] = getEta(0.5*(rHB[17]+rHB[18]), zho[1]);
1100  etaHO[1] = getEta(rHB[18]+drHB[18], zho[2]);
1101  etaHO[2] = getEta(rHB[18]-drHB[18], zho[3]);
1102  etaHO[3] = getEta(rHB[18]+drHB[18], zho[4]);
1103  } else {
1104  rminHO =-1.0;
1105  etaHO[0] = etaTable[4];
1106  etaHO[1] = etaTable[4];
1107  etaHO[2] = etaTable[10];
1108  etaHO[3] = etaTable[10];
1109  }
1110 #ifdef DebugLog
1111  LogDebug("HCalGeom") << "HO Eta boundaries " << etaHO[0] << " " << etaHO[1]
1112  << " " << etaHO[2] << " " << etaHO[3];
1113  std::cout << "HO Parameters " << rminHO << " " << zho.size();
1114  for (int i=0; i<4; ++i) std::cout << " eta[" << i << "] = " << etaHO[i];
1115  for (unsigned int i=0; i<zho.size(); ++i) std::cout << " zho[" << i << "] = " << zho[i];
1116  std::cout << std::endl;
1117 #endif
1118 }
1119 
1120 std::vector<double> HcalNumberingFromDDD::getDDDArray(const std::string & str,
1121  const DDsvalues_type & sv,
1122  int & nmin) const {
1123 #ifdef DebugLog
1124  LogDebug("HCalGeom") << "HcalNumberingFromDDD:getDDDArray called for "
1125  << str << " with nMin " << nmin;
1126 #endif
1127  DDValue value(str);
1128  if (DDfetch(&sv,value)) {
1129 #ifdef DebugLog
1130  LogDebug("HCalGeom") << "HcalNumberingFromDDD: " << value;
1131 #endif
1132  const std::vector<double> & fvec = value.doubles();
1133  int nval = fvec.size();
1134  if (nmin > 0) {
1135  if (nval < nmin) {
1136  edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str
1137  << " bins " << nval << " < " << nmin
1138  << " ==> illegal";
1139  throw cms::Exception("DDException") << "HcalNumberingFromDDD: cannot get array " << str;
1140  }
1141  } else {
1142  if (nval < 2) {
1143  edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str
1144  << " bins " << nval << " < 2 ==> illegal"
1145  << " (nmin=" << nmin << ")";
1146  throw cms::Exception("DDException") << "HcalNumberingFromDDD: cannot get array " << str;
1147  }
1148  }
1149  nmin = nval;
1150  return fvec;
1151  } else {
1152  edm::LogError("HCalGeom") << "HcalNumberingFromDDD: cannot get array "
1153  << str;
1154  throw cms::Exception("DDException") << "HcalNumberingFromDDD: cannot get array " << str;
1155  }
1156 }
1157 
1158 int HcalNumberingFromDDD::getShift(HcalSubdetector subdet, int depth) const {
1159 
1160  int shift;
1161  switch(subdet) {
1162  case HcalEndcap:
1163  shift = shiftHE[depth-1];
1164  break;
1165  case HcalForward:
1166  shift = shiftHF[depth-1];
1167  break;
1168  default:
1169  shift = shiftHB[depth-1];
1170  break;
1171  }
1172  return shift;
1173 }
1174 
1175 double HcalNumberingFromDDD::getGain(HcalSubdetector subdet, int depth) const {
1176 
1177  double gain;
1178  switch(subdet) {
1179  case HcalEndcap:
1180  gain = gainHE[depth-1];
1181  break;
1182  case HcalForward:
1183  gain = gainHF[depth-1];
1184  break;
1185  default:
1186  gain = gainHB[depth-1];
1187  break;
1188  }
1189  return gain;
1190 }
1191 
1192 unsigned int HcalNumberingFromDDD::find(int element,
1193  std::vector<int> array) const {
1194 
1195  unsigned int id = array.size();
1196  for (unsigned int i = 0; i < array.size(); i++) {
1197  if (element == array[i]) {
1198  id = i;
1199  break;
1200  }
1201  }
1202  return id;
1203 }
1204 
1205 int HcalNumberingFromDDD::unitPhi(int det, int etaR) const {
1206 
1207  const double fiveDegInRad = 2*M_PI/72;
1208  int units=0;
1209  if (det == static_cast<int>(HcalForward))
1210  units=int(phibin[nEta+etaR-etaMin[2]-1]/fiveDegInRad+0.5);
1211  else
1212  units=int(phibin[etaR-1]/fiveDegInRad+0.5);
1213 
1214  return units;
1215 }
1216 
1217 void HcalNumberingFromDDD::tileHB(int eta, int depth) {
1218 
1219  double etaL = etaTable[eta-1];
1220  double thetaL = 2.*atan(exp(-etaL));
1221  double etaH = etaTable[eta];
1222  double thetaH = 2.*atan(exp(-etaH));
1223  int layL=0, layH=0;
1224  if (depth == 1) {
1225  layH = depth1[eta-1];
1226  } else {
1227  layL = depth1[eta-1];
1228  layH = depth2[eta-1];
1229  }
1230  std::cout << "\ntileHB:: eta|depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << "\n";
1231  for (int lay=layL; lay<layH; ++lay) {
1232  std::vector<double> area(2,0);
1233  int kk=0;
1234  for (unsigned int k=0; k<layb.size(); ++k) {
1235  if (lay == layb[k]) {
1236  double zmin = rhoxb[k]*std::cos(thetaL)/std::sin(thetaL);
1237  double zmax = rhoxb[k]*std::cos(thetaH)/std::sin(thetaH);
1238  double dz = (std::min(zmax,dzxb[k]) - zmin);
1239  if (dz > 0) {
1240  area[kk] = dz*dyxb[k];
1241  kk++;
1242  }
1243  }
1244  }
1245  if (area[0] > 0) std::cout << std::setw(2) << lay << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
1246  }
1247 }
1248 
1249 void HcalNumberingFromDDD::tileHE(int eta, int depth) {
1250 
1251  double etaL = etaTable[eta-1];
1252  double thetaL = 2.*atan(exp(-etaL));
1253  double etaH = etaTable[eta];
1254  double thetaH = 2.*atan(exp(-etaH));
1255  int layL=0, layH=0;
1256  if (eta == 16) {
1257  layH = depth3[eta-1];
1258  } else if (depth == 1) {
1259  layH = depth1[eta-1];
1260  } else if (depth == 2) {
1261  layL = depth1[eta-1];
1262  layH = depth2[eta-1];
1263  } else {
1264  layL = depth2[eta-1];
1265  layH = depth3[eta-1];
1266  }
1267  double phib = phibin[eta-1];
1268  int nphi = 2;
1269  if (phib > 6*CLHEP::deg) nphi = 1;
1270  std::cout << "\ntileHE:: Eta/depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << " phi " << nphi << "\n";
1271  for (int lay=layL; lay<layH; ++lay) {
1272  std::vector<double> area(4,0);
1273  int kk=0;
1274  for (unsigned int k=0; k<laye.size(); ++k) {
1275  if (lay == laye[k]) {
1276  double rmin = zxe[k]*std::tan(thetaH);
1277  double rmax = zxe[k]*std::tan(thetaL);
1278  if ((lay != 0 || eta == 18) &&
1279  (lay != 1 || (eta == 18 && rhoxe[k]-dyxe[k] > 1000) || (eta != 18 && rhoxe[k]-dyxe[k] < 1000)) &&
1280  rmin+30 < rhoxe[k]+dyxe[k] && rmax > rhoxe[k]-dyxe[k]) {
1281  rmin = std::max(rmin,rhoxe[k]-dyxe[k]);
1282  rmax = std::min(rmax,rhoxe[k]+dyxe[k]);
1283  double dx1 = rmin*std::tan(phib);
1284  double dx2 = rmax*std::tan(phib);
1285  double ar1=0, ar2=0;
1286  if (nphi == 1) {
1287  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-4.*dx1e[k]);
1288  } else {
1289  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-2.*dx1e[k]);
1290  ar2 = 0.5*(rmax-rmin)*((rmax+rmin)*tan(10.*CLHEP::deg)-4*dx1e[k])-ar1;
1291  }
1292  area[kk] = ar1;
1293  area[kk+2] = ar2;
1294  kk++;
1295  }
1296  }
1297  }
1298  if (area[0] > 0 && area[1] > 0) {
1299  int lay0 = lay-1;
1300  if (eta == 18) lay0++;
1301  if (nphi == 1) {
1302  std::cout << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
1303  } else {
1304  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";
1305  }
1306  }
1307  }
1308 }
1309 
1310 double HcalNumberingFromDDD::getEtaHO(double& etaR, double& x, double& y,
1311  double& z) const {
1312 
1313  double eta = fabs(etaR);
1314  double r = std::sqrt(x*x+y*y);
1315  if (r > rminHO) {
1316  double zz = fabs(z);
1317  if (zz > zho[3]) {
1318  if (eta <= etaTable[10]) eta = etaTable[10]+0.001;
1319  } else if (zz > zho[1]) {
1320  if (eta <= etaTable[4]) eta = etaTable[4]+0.001;
1321  }
1322  }
1323  eta = (z >= 0. ? eta : -eta);
1324  // std::cout << "R " << r << " Z " << z << " eta " << etaR << ":" << eta <<"\n";
1325  // if (eta != etaR) std::cout << "**** Check *****\n";
1326  return eta;
1327 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
HcalNumberingFromDDD::shiftHE
std::vector< int > shiftHE
Definition: HcalNumberingFromDDD.h:83
HcalNumberingFromDDD::shiftHB
std::vector< int > shiftHB
Definition: HcalNumberingFromDDD.h:81
HcalNumberingFromDDD::dx2e
std::vector< double > dx2e
Definition: HcalNumberingFromDDD.h:98
i
int i
Definition: DBlmapReader.cc:9
DDSolid::parameters
const std::vector< double > & parameters(void) const
Give the parameters of the solid.
Definition: DDSolid.cc:150
create_public_lumi_plots.units
tuple units
Definition: create_public_lumi_plots.py:946
DDTrap::halfZ
double halfZ(void) const
half of the z-Axis
Definition: DDSolid.cc:167
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:134
DDBox::halfZ
double halfZ(void) const
Definition: DDSolid.cc:258
DDFilteredView::addFilter
void addFilter(const DDFilter &, log_op op=AND)
Definition: DDFilteredView.cc:24
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1093
DDTrap::x1
double x1(void) const
Half-length along x of the side at y=-pDy1 of the face at -pDz.
Definition: DDSolid.cc:175
DDBase::name
const N & name() const
Definition: DDBase.h:82
DDBox::halfY
double halfY(void) const
Definition: DDSolid.cc:256
HcalNumberingFromDDD::depth2
std::vector< int > depth2
Definition: HcalNumberingFromDDD.h:77
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:98
HcalNumberingFromDDD::tileHB
void tileHB(int eta, int depth)
Definition: HcalNumberingFromDDD.cc:1217
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
HcalCellType::setMissingPhi
void setMissingPhi(std::vector< int >, std::vector< int >)
Definition: HcalCellType.cc:57
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1108
HcalNumberingFromDDD::zxe
std::vector< double > zxe
Definition: HcalNumberingFromDDD.h:98
HcalNumberingFromDDD::etaTable
std::vector< double > etaTable
Definition: HcalNumberingFromDDD.h:68
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:165
HcalNumberingFromDDD::layb
std::vector< int > layb
Definition: HcalNumberingFromDDD.h:97
abs
#define abs(x)
Definition: mlp_lapack.h:159
HcalNumberingFromDDD::getShift
int getShift(HcalSubdetector subdet, int depth) const
Definition: HcalNumberingFromDDD.cc:1158
HcalNumberingFromDDD::zxb
std::vector< double > zxb
Definition: HcalNumberingFromDDD.h:96
min
#define min(a, b)
Definition: mlp_lapack.h:161
HcalNumberingFromDDD::unitPhi
int unitPhi(int det, int etaR) const
Definition: HcalNumberingFromDDD.cc:1205
HcalNumberingFromDDD::rHB
std::vector< double > rHB
Definition: HcalNumberingFromDDD.h:90
DDLogicalPart::solid
const DDSolid & solid(void) const
Returns a reference object of the solid being the shape of this LogicalPart.
Definition: DDLogicalPart.cc:149
HcalNumberingFromDDD::HcalCellTypes
std::vector< HcalCellType > HcalCellTypes() const
Definition: HcalNumberingFromDDD.cc:391
alcazmumu_cfi.filter
tuple filter
Definition: alcazmumu_cfi.py:34
HcalNumberingFromDDD::drHB
std::vector< double > drHB
Definition: HcalNumberingFromDDD.h:90
eta
T eta() const
Definition: Basic3DVectorLD.h:172
HcalNumberingFromDDD::getGain
double getGain(HcalSubdetector subdet, int depth) const
Definition: HcalNumberingFromDDD.cc:1175
DDCompactView
type of data representation of DDCompactView
Definition: DDCompactView.h:77
HcalNumberingFromDDD::phibin
std::vector< double > phibin
Definition: HcalNumberingFromDDD.h:74
DDfetch
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:102
detailsBasic3DVector::z
double double double z
Definition: newBasic3DVector.h:17
HcalNumberingFromDDD::rhoxe
std::vector< double > rhoxe
Definition: HcalNumberingFromDDD.h:98
DDSolid
A DDSolid represents the shape of a part.
Definition: DDSolid.h:35
HcalNumberingFromDDD::numberOfCells
unsigned int numberOfCells(HcalSubdetector) const
Definition: HcalNumberingFromDDD.cc:373
HcalNumberingFromDDD::zho
std::vector< double > zho
Definition: HcalNumberingFromDDD.h:92
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:21
HcalNumberingFromDDD::dzxb
std::vector< double > dzxb
Definition: HcalNumberingFromDDD.h:96
errorMatrix2Lands_multiChannel.id
tuple id
Definition: errorMatrix2Lands_multiChannel.py:151
HcalNumberingFromDDD::gainHF
std::vector< double > gainHF
Definition: HcalNumberingFromDDD.h:84
HcalNumberingFromDDD::shiftHF
std::vector< int > shiftHF
Definition: HcalNumberingFromDDD.h:85
DDTrap::x4
double x4(void) const
Half-length along x of the side at y=+pDy2 of the face at +pDz.
Definition: DDSolid.cc:185
HcalNumberingFromDDD::dyxb
std::vector< double > dyxb
Definition: HcalNumberingFromDDD.h:96
GetRecoTauVFromDQM_MC_cff.kk
tuple kk
Definition: GetRecoTauVFromDQM_MC_cff.py:83
HcalNumberingFromDDD::getEtaHO
double getEtaHO(double &etaR, double &x, double &y, double &z) const
Definition: HcalNumberingFromDDD.cc:1310
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
HcalNumberingFromDDD::phioff
std::vector< double > phioff
Definition: HcalNumberingFromDDD.h:67
HcalNumberingFromDDD::dx1e
std::vector< double > dx1e
Definition: HcalNumberingFromDDD.h:98
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
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:97
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:77
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
DDTrap::y1
double y1(void) const
Half-length along y of the face at -pDz.
Definition: DDSolid.cc:173
HcalSubdetector
HcalSubdetector
Definition: HcalAssistant.h:32
j
int j
Definition: DBlmapReader.cc:9
DDSolid::shape
DDSolidShape shape(void) const
The type of the solid.
Definition: DDSolid.cc:144
HcalNumberingFromDDD::nOff
std::vector< int > nOff
Definition: HcalNumberingFromDDD.h:89
DDTubs::rOut
double rOut(void) const
Definition: DDSolid.cc:509
HcalNumberingFromDDD::depth3
std::vector< int > depth3
Definition: HcalNumberingFromDDD.h:78
corr
JetCorrectorParameters corr
Definition: classes.h:9
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
HcalNumberingFromDDD::rhoxb
std::vector< double > rhoxb
Definition: HcalNumberingFromDDD.h:96
HcalNumberingFromDDD::find
unsigned find(int element, std::vector< int > array) const
Definition: HcalNumberingFromDDD.cc:1192
HcalNumberingFromDDD::etaMin
std::vector< int > etaMin
Definition: HcalNumberingFromDDD.h:72
HcalNumberingFromDDD::etaMax
std::vector< int > etaMax
Definition: HcalNumberingFromDDD.h:73
DDBox::halfX
double halfX(void) const
Definition: DDSolid.cc:254
M_PI
#define M_PI
Definition: BFit3D.cc:3
HcalNumberingFromDDD::loadGeometry
void loadGeometry(DDFilteredView)
Definition: HcalNumberingFromDDD.cc:817
HcalNumberingFromDDD::deltaEta
double deltaEta(int det, int eta, int depth) const
Definition: HcalNumberingFromDDD.cc:571
HcalNumberingFromDDD::getDDDArray
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &) const
Definition: HcalNumberingFromDDD.cc:1120
DDBox
Interface to a Box.
Definition: DDSolid.h:162
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:80
evf::utils::num
long long int num
Definition: procUtils.cc:71
HcalNumberingFromDDD::gainHE
std::vector< double > gainHE
Definition: HcalNumberingFromDDD.h:82
HcalNumberingFromDDD::depth1
std::vector< int > depth1
Definition: HcalNumberingFromDDD.h:76
DDFilteredView::mergedSpecifics
DDsvalues_type mergedSpecifics() const
Definition: DDFilteredView.cc:66
HcalNumberingFromDDD::zHE
std::vector< double > zHE
Definition: HcalNumberingFromDDD.h:91
HcalNumberingFromDDD::tileHE
void tileHE(int eta, int depth)
Definition: HcalNumberingFromDDD.cc:1249
DDTrap::x2
double x2(void) const
Half-length along x of the side at y=+pDy1 of the face at -pDz.
Definition: DDSolid.cc:177
DDTubs::zhalf
double zhalf(void) const
Definition: DDSolid.cc:505
HcalNumberingFromDDD::dzHE
std::vector< double > dzHE
Definition: HcalNumberingFromDDD.h:91
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
DDFilteredView::firstChild
bool firstChild()
set the current node to the first child ...
Definition: DDFilteredView.cc:131
DDTrap::y2
double y2(void) const
Half-length along y of the face at +pDz.
Definition: DDSolid.cc:181
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:121
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:31
x
Definition: DDAxes.h:10
DDTubs::rIn
double rIn(void) const
Definition: DDSolid.cc:507
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
DDTrap::x3
double x3(void) const
Half-length along x of the side at y=-pDy2 of the face at +pDz.
Definition: DDSolid.cc:183
HcalNumberingFromDDD::rTable
std::vector< double > rTable
Definition: HcalNumberingFromDDD.h:70
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