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  const 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(const 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(const DDFilteredView& _fv) {
818 
819  DDFilteredView fv = _fv;
820  bool dodet=true, hf=false;
821  std::vector<double> rb(20,0.0), ze(20,0.0), thkb(20,-1.0), thke(20,-1.0);
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 #ifdef DebugLog
932  LogDebug("HCalGeom") << "Detector " << idet << " Lay " << lay << " fi " << ifi << " " << ich << " z " << z1 << " " << z2;
933 #endif
934  }
935  }
936  } else if (idet == 4) {
937  // HE
938 #ifdef DebugLog
939  LogDebug("HCalGeom") << "HE " << sol.name() << " Shape " << sol.shape()
940  << " Layer " << lay << " Z " << t.z();
941 #endif
942  if (lay >=0 && lay < 20) {
943  ie[lay]++;
944  ze[lay] += std::abs(t.z());
945  if (thke[lay] <= 0) thke[lay] = dz;
946  bool found = false;
947  for (unsigned int k=0; k<zxe.size(); k++) {
948  if (std::abs(zxe[k]-std::abs(t.z())) < 0.01) {
949  found = true;
950  break;
951  }
952  }
953  if (!found) {
954  zxe.push_back(std::abs(t.z()));
955  rhoxe.push_back(t.Rho()*std::cos(t.phi()));
956  dyxe.push_back(dy*std::cos(t.phi()));
957  dx1 -= 0.5*(t.rho()-dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
958  dx2 -= 0.5*(t.rho()+dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
959  dx1e.push_back(-dx1);
960  dx2e.push_back(-dx2);
961  laye.push_back(lay);
962  }
963  }
964  if (copy[nsiz-1] == 21) {
965  int iz = copy[nsiz-7];
966  int fi = copy[nsiz-5];
967  unsigned int it1 = find(iz, ize);
968  if (it1 == ize.size()) ize.push_back(iz);
969  unsigned int it2 = find(fi, phie);
970  if (it2 == phie.size()) phie.push_back(fi);
971  }
972  } else if (idet == 5) {
973  // HF
974  if (!hf) {
975  const std::vector<double> & paras = sol.parameters();
976 #ifdef DebugLog
977  LogDebug("HCalGeom") << "HF " << sol.name() << " Shape " << sol.shape()
978  << " Z " << t.z() << " with " << paras.size()
979  << " Parameters";
980  for (unsigned j=0; j<paras.size(); j++)
981  LogDebug("HCalGeom") << "HF Parameter[" << j << "] = " << paras[j];
982 #endif
983  zf = fabs(t.z());
984  if (sol.shape() == ddpolycone_rrz) {
985  int nz = (int)(paras.size())-3;
986  zf += paras[3];
987  dzVcal = 0.5*(paras[nz]-paras[3]);
988  hf = true;
989  } else if (sol.shape() == ddtubs || sol.shape() == ddcons) {
990  dzVcal = paras[0];
991  zf -= paras[0];
992  hf = true;
993  }
994  }
995 #ifdef DebugLog
996  } else {
997  LogDebug("HCalGeom") << "Unknown Detector " << idet << " for "
998  << sol.name() << " Shape " << sol.shape() << " R "
999  << t.Rho() << " Z " << t.z();
1000 #endif
1001  }
1002  dodet = fv.next();
1003  }
1004 
1005  int ibmx = 0, iemx = 0;
1006  for (int i = 0; i < 20; i++) {
1007  if (ib[i]>0) {
1008  rb[i] /= (double)(ib[i]);
1009  ibmx = i+1;
1010  }
1011  if (ie[i]>0) {
1012  ze[i] /= (double)(ie[i]);
1013  iemx = i+1;
1014  }
1015 #ifdef DebugLog
1016  LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
1017  << rb[i] << " Endcap " << ie[i] << " " << ze[i];
1018 #endif
1019  }
1020  for (int i = 4; i >= 0; i--) {
1021  if (ib[i] == 0) {rb[i] = rb[i+1]; thkb[i] = thkb[i+1];}
1022  if (ie[i] == 0) {ze[i] = ze[i+1]; thke[i] = thke[i+1];}
1023 #ifdef DebugLog
1024  if (ib[i] == 0 || ie[i] == 0)
1025  LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
1026  << rb[i] << " Endcap " << ie[i] << " " << ze[i];
1027 #endif
1028  }
1029 
1030 #ifdef DebugLog
1031  for (unsigned int k=0; k<layb.size(); ++k)
1032  std::cout << "HB: " << layb[k] << " R " << rxb[k] << " " << rhoxb[k] << " Z " << zxb[k] << " DY " << dyxb[k] << " DZ " << dzxb[k] << "\n";
1033  for (unsigned int k=0; k<laye.size(); ++k)
1034  std::cout << "HE: " << laye[k] << " R " << rhoxe[k] << " Z " << zxe[k] << " X1|X2 " << dx1e[k] << "|" << dx2e[k] << " DY " << dyxe[k] << "\n";
1035 
1036  printTile();
1037  LogDebug("HCalGeom") << "HcalNumberingFromDDD: Maximum Layer for HB "
1038  << ibmx << " for HE " << iemx << " Z for HF " << zf
1039  << " extent " << dzVcal;
1040 #endif
1041 
1042  if (ibmx > 0) {
1043  rHB.resize(ibmx);
1044  drHB.resize(ibmx);
1045  for (int i=0; i<ibmx; i++) {
1046  rHB[i] = rb[i];
1047  drHB[i] = thkb[i];
1048 #ifdef DebugLog
1049  LogDebug("HCalGeom") << "HcalNumberingFromDDD: rHB[" << i << "] = "
1050  << rHB[i] << " drHB[" << i << "] = " << drHB[i];
1051 #endif
1052  }
1053  }
1054  if (iemx > 0) {
1055  zHE.resize(iemx);
1056  dzHE.resize(iemx);
1057  for (int i=0; i<iemx; i++) {
1058  zHE[i] = ze[i];
1059  dzHE[i] = thke[i];
1060 #ifdef DebugLog
1061  LogDebug("HCalGeom") << "HcalNumberingFromDDD: zHE[" << i << "] = "
1062  << zHE[i] << " dzHE[" << i << "] = " << dzHE[i];
1063 #endif
1064  }
1065  }
1066 
1067  nzHB = (int)(izb.size());
1068  nmodHB = (int)(phib.size());
1069 #ifdef DebugLog
1070  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHB
1071  << " barrel half-sectors";
1072  for (int i=0; i<nzHB; i++)
1073  LogDebug("HCalGeom") << "Section " << i << " Copy number " << izb[i];
1074  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHB
1075  << " barrel modules";
1076  for (int i=0; i<nmodHB; i++)
1077  LogDebug("HCalGeom") << "Module " << i << " Copy number " << phib[i];
1078 #endif
1079 
1080  nzHE = (int)(ize.size());
1081  nmodHE = (int)(phie.size());
1082 #ifdef DebugLog
1083  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHE
1084  << " endcap half-sectors";
1085  for (int i=0; i<nzHE; i++)
1086  LogDebug("HCalGeom") << "Section " << i << " Copy number " << ize[i];
1087  LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHE
1088  << " endcap modules";
1089  for (int i=0; i<nmodHE; i++)
1090  LogDebug("HCalGeom") << "Module " << i << " Copy number " << phie[i];
1091 #endif
1092 
1093 #ifdef DebugLog
1094  LogDebug("HCalGeom") << "HO has Z of size " << zho.size();
1095  for (unsigned int kk=0; kk<zho.size(); kk++)
1096  LogDebug("HCalGeom") << "ZHO[" << kk << "] = " << zho[kk];
1097 #endif
1098  if (ibmx > 17 && zho.size() > 4) {
1099  rminHO = rHB[17]-100.0;
1100  etaHO[0] = getEta(0.5*(rHB[17]+rHB[18]), zho[1]);
1101  etaHO[1] = getEta(rHB[18]+drHB[18], zho[2]);
1102  etaHO[2] = getEta(rHB[18]-drHB[18], zho[3]);
1103  etaHO[3] = getEta(rHB[18]+drHB[18], zho[4]);
1104  } else {
1105  rminHO =-1.0;
1106  etaHO[0] = etaTable[4];
1107  etaHO[1] = etaTable[4];
1108  etaHO[2] = etaTable[10];
1109  etaHO[3] = etaTable[10];
1110  }
1111 #ifdef DebugLog
1112  LogDebug("HCalGeom") << "HO Eta boundaries " << etaHO[0] << " " << etaHO[1]
1113  << " " << etaHO[2] << " " << etaHO[3];
1114  std::cout << "HO Parameters " << rminHO << " " << zho.size();
1115  for (int i=0; i<4; ++i) std::cout << " eta[" << i << "] = " << etaHO[i];
1116  for (unsigned int i=0; i<zho.size(); ++i) std::cout << " zho[" << i << "] = " << zho[i];
1117  std::cout << std::endl;
1118 #endif
1119 }
1120 
1121 std::vector<double> HcalNumberingFromDDD::getDDDArray(const std::string & str,
1122  const DDsvalues_type & sv,
1123  int & nmin) const {
1124 #ifdef DebugLog
1125  LogDebug("HCalGeom") << "HcalNumberingFromDDD:getDDDArray called for "
1126  << str << " with nMin " << nmin;
1127 #endif
1128  DDValue value(str);
1129  if (DDfetch(&sv,value)) {
1130 #ifdef DebugLog
1131  LogDebug("HCalGeom") << "HcalNumberingFromDDD: " << value;
1132 #endif
1133  const std::vector<double> & fvec = value.doubles();
1134  int nval = fvec.size();
1135  if (nmin > 0) {
1136  if (nval < nmin) {
1137  edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str
1138  << " bins " << nval << " < " << nmin
1139  << " ==> illegal";
1140  throw cms::Exception("DDException") << "HcalNumberingFromDDD: cannot get array " << str;
1141  }
1142  } else {
1143  if (nval < 2) {
1144  edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str
1145  << " bins " << nval << " < 2 ==> illegal"
1146  << " (nmin=" << nmin << ")";
1147  throw cms::Exception("DDException") << "HcalNumberingFromDDD: cannot get array " << str;
1148  }
1149  }
1150  nmin = nval;
1151  return fvec;
1152  } else {
1153  edm::LogError("HCalGeom") << "HcalNumberingFromDDD: cannot get array "
1154  << str;
1155  throw cms::Exception("DDException") << "HcalNumberingFromDDD: cannot get array " << str;
1156  }
1157 }
1158 
1159 int HcalNumberingFromDDD::getShift(HcalSubdetector subdet, int depth) const {
1160 
1161  int shift;
1162  switch(subdet) {
1163  case HcalEndcap:
1164  shift = shiftHE[depth-1];
1165  break;
1166  case HcalForward:
1167  shift = shiftHF[depth-1];
1168  break;
1169  default:
1170  shift = shiftHB[depth-1];
1171  break;
1172  }
1173  return shift;
1174 }
1175 
1176 double HcalNumberingFromDDD::getGain(HcalSubdetector subdet, int depth) const {
1177 
1178  double gain;
1179  switch(subdet) {
1180  case HcalEndcap:
1181  gain = gainHE[depth-1];
1182  break;
1183  case HcalForward:
1184  gain = gainHF[depth-1];
1185  break;
1186  default:
1187  gain = gainHB[depth-1];
1188  break;
1189  }
1190  return gain;
1191 }
1192 
1193 unsigned int HcalNumberingFromDDD::find(int element,
1194  std::vector<int>& array) const {
1195 
1196  unsigned int id = array.size();
1197  for (unsigned int i = 0; i < array.size(); i++) {
1198  if (element == array[i]) {
1199  id = i;
1200  break;
1201  }
1202  }
1203  return id;
1204 }
1205 
1206 int HcalNumberingFromDDD::unitPhi(int det, int etaR) const {
1207 
1208  const double fiveDegInRad = 2*M_PI/72;
1209  int units=0;
1210  if (det == static_cast<int>(HcalForward))
1211  units=int(phibin[nEta+etaR-etaMin[2]-1]/fiveDegInRad+0.5);
1212  else
1213  units=int(phibin[etaR-1]/fiveDegInRad+0.5);
1214 
1215  return units;
1216 }
1217 
1218 void HcalNumberingFromDDD::tileHB(int eta, int depth) {
1219 
1220  double etaL = etaTable[eta-1];
1221  double thetaL = 2.*atan(exp(-etaL));
1222  double etaH = etaTable[eta];
1223  double thetaH = 2.*atan(exp(-etaH));
1224  int layL=0, layH=0;
1225  if (depth == 1) {
1226  layH = depth1[eta-1];
1227  } else {
1228  layL = depth1[eta-1];
1229  layH = depth2[eta-1];
1230  }
1231  std::cout << "\ntileHB:: eta|depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << "\n";
1232  for (int lay=layL; lay<layH; ++lay) {
1233  std::vector<double> area(2,0);
1234  int kk=0;
1235  for (unsigned int k=0; k<layb.size(); ++k) {
1236  if (lay == layb[k]) {
1237  double zmin = rhoxb[k]*std::cos(thetaL)/std::sin(thetaL);
1238  double zmax = rhoxb[k]*std::cos(thetaH)/std::sin(thetaH);
1239  double dz = (std::min(zmax,dzxb[k]) - zmin);
1240  if (dz > 0) {
1241  area[kk] = dz*dyxb[k];
1242  kk++;
1243  }
1244  }
1245  }
1246  if (area[0] > 0) std::cout << std::setw(2) << lay << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
1247  }
1248 }
1249 
1250 void HcalNumberingFromDDD::tileHE(int eta, int depth) {
1251 
1252  double etaL = etaTable[eta-1];
1253  double thetaL = 2.*atan(exp(-etaL));
1254  double etaH = etaTable[eta];
1255  double thetaH = 2.*atan(exp(-etaH));
1256  int layL=0, layH=0;
1257  if (eta == 16) {
1258  layH = depth3[eta-1];
1259  } else if (depth == 1) {
1260  layH = depth1[eta-1];
1261  } else if (depth == 2) {
1262  layL = depth1[eta-1];
1263  layH = depth2[eta-1];
1264  } else {
1265  layL = depth2[eta-1];
1266  layH = depth3[eta-1];
1267  }
1268  double phib = phibin[eta-1];
1269  int nphi = 2;
1270  if (phib > 6*CLHEP::deg) nphi = 1;
1271  std::cout << "\ntileHE:: Eta/depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << " phi " << nphi << "\n";
1272  for (int lay=layL; lay<layH; ++lay) {
1273  std::vector<double> area(4,0);
1274  int kk=0;
1275  for (unsigned int k=0; k<laye.size(); ++k) {
1276  if (lay == laye[k]) {
1277  double rmin = zxe[k]*std::tan(thetaH);
1278  double rmax = zxe[k]*std::tan(thetaL);
1279  if ((lay != 0 || eta == 18) &&
1280  (lay != 1 || (eta == 18 && rhoxe[k]-dyxe[k] > 1000) || (eta != 18 && rhoxe[k]-dyxe[k] < 1000)) &&
1281  rmin+30 < rhoxe[k]+dyxe[k] && rmax > rhoxe[k]-dyxe[k]) {
1282  rmin = std::max(rmin,rhoxe[k]-dyxe[k]);
1283  rmax = std::min(rmax,rhoxe[k]+dyxe[k]);
1284  double dx1 = rmin*std::tan(phib);
1285  double dx2 = rmax*std::tan(phib);
1286  double ar1=0, ar2=0;
1287  if (nphi == 1) {
1288  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-4.*dx1e[k]);
1289  } else {
1290  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-2.*dx1e[k]);
1291  ar2 = 0.5*(rmax-rmin)*((rmax+rmin)*tan(10.*CLHEP::deg)-4*dx1e[k])-ar1;
1292  }
1293  area[kk] = ar1;
1294  area[kk+2] = ar2;
1295  kk++;
1296  }
1297  }
1298  }
1299  if (area[0] > 0 && area[1] > 0) {
1300  int lay0 = lay-1;
1301  if (eta == 18) lay0++;
1302  if (nphi == 1) {
1303  std::cout << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
1304  } else {
1305  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";
1306  }
1307  }
1308  }
1309 }
1310 
1311 double HcalNumberingFromDDD::getEtaHO(double& etaR, double& x, double& y,
1312  double& z) const {
1313 
1314  double eta = fabs(etaR);
1315  double r = std::sqrt(x*x+y*y);
1316  if (r > rminHO) {
1317  double zz = fabs(z);
1318  if (zz > zho[3]) {
1319  if (eta <= etaTable[10]) eta = etaTable[10]+0.001;
1320  } else if (zz > zho[1]) {
1321  if (eta <= etaTable[4]) eta = etaTable[4]+0.001;
1322  }
1323  }
1324  eta = (z >= 0. ? eta : -eta);
1325  // std::cout << "R " << r << " Z " << z << " eta " << etaR << ":" << eta <<"\n";
1326  // if (eta != etaR) std::cout << "**** Check *****\n";
1327  return eta;
1328 }
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
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:139
HcalCellType::setMissingPhi
void setMissingPhi(std::vector< int > &, std::vector< int > &)
Definition: HcalCellType.cc:57
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:1087
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
cuy.ib
int ib
Definition: cuy.py:660
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:1218
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1102
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:170
HcalNumberingFromDDD::layb
std::vector< int > layb
Definition: HcalNumberingFromDDD.h:97
HcalNumberingFromDDD::getShift
int getShift(HcalSubdetector subdet, int depth) const
Definition: HcalNumberingFromDDD.cc:1159
HcalNumberingFromDDD::zxb
std::vector< double > zxb
Definition: HcalNumberingFromDDD.h:96
HcalNumberingFromDDD::unitPhi
int unitPhi(int det, int etaR) const
Definition: HcalNumberingFromDDD.cc:1206
HcalNumberingFromDDD::rHB
std::vector< double > rHB
Definition: HcalNumberingFromDDD.h:90
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:210
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:177
HcalNumberingFromDDD::getGain
double getGain(HcalSubdetector subdet, int depth) const
Definition: HcalNumberingFromDDD.cc:1176
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:80
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:15
HcalNumberingFromDDD::rhoxe
std::vector< double > rhoxe
Definition: HcalNumberingFromDDD.h:98
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
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
HcalNumberingFromDDD::gainHF
std::vector< double > gainHF
Definition: HcalNumberingFromDDD.h:84
roll_playback.k
tuple k
Definition: roll_playback.py:143
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:1311
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:48
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
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:31
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
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
min
T min(T a, T b)
Definition: MathUtil.h:58
DDTubs::rOut
double rOut(void) const
Definition: DDSolid.cc:509
M_PI
#define M_PI
Definition: readTestVector.cc:25
HcalNumberingFromDDD::depth3
std::vector< int > depth3
Definition: HcalNumberingFromDDD.h:78
corr
JetCorrectorParameters corr
Definition: classes.h:5
HcalNumberingFromDDD::loadSpecPars
void loadSpecPars(const DDFilteredView &)
Definition: HcalNumberingFromDDD.cc:645
HcalNumberingFromDDD::rhoxb
std::vector< double > rhoxb
Definition: HcalNumberingFromDDD.h:96
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
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:1121
SiStripMonitorClusterAlca_cfi.zmin
tuple zmin
Definition: SiStripMonitorClusterAlca_cfi.py:182
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
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:1250
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
units
TString units(TString variable, Char_t axis)
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
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
HcalNumberingFromDDD::find
unsigned find(int element, std::vector< int > &array) const
Definition: HcalNumberingFromDDD.cc:1193
HcalNumberingFromDDD::loadGeometry
void loadGeometry(const DDFilteredView &)
Definition: HcalNumberingFromDDD.cc:817
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
HcalNumberingFromDDD::getEtaTable
std::vector< double > getEtaTable() const
Definition: HcalNumberingFromDDD.cc:367
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
x
Definition: DDAxes.h:10
DDTubs::rIn
double rIn(void) const
Definition: DDSolid.cc:507
HcalNumberingFromDDD::unitID
HcalID unitID(int det, const CLHEP::Hep3Vector &pos, int depth, int lay=-1) const
Definition: HcalNumberingFromDDD.cc:31
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
SiStripMonitorClusterAlca_cfi.zmax
tuple zmax
Definition: SiStripMonitorClusterAlca_cfi.py: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