test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalDDDSimConstants.cc
Go to the documentation of this file.
2 
5 
6 #include "CLHEP/Units/GlobalPhysicalConstants.h"
7 #include "CLHEP/Units/GlobalSystemOfUnits.h"
8 
9 //#define DebugLog
10 
12 
13 #ifdef DebugLog
14  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants::HcalDDDSimConstants (const HcalParameter* hp) constructor";
15 #endif
16 
17  initialize();
18 #ifdef DebugLog
19  std::vector<HcalCellType> cellTypes = HcalCellTypes();
20  edm::LogInfo ("HCalGeom") << "HcalDDDSimConstants: " << cellTypes.size()
21  << " cells of type HCal (All)";
22 #endif
23 }
24 
25 
27 #ifdef DebugLog
28  edm::LogInfo ("HCalGeom") << "HcalDDDSimConstants::destructed!!!";
29 #endif
30 }
31 
33  int depth, int etaR,
34  int iphi) const {
35 
36  double etaMn = hpar->etaMin[0];
37  double etaMx = hpar->etaMax[0];
38  if (idet==static_cast<int>(HcalEndcap)) {
39  etaMn = hpar->etaMin[1]; etaMx = hpar->etaMax[1];
40  } else if (idet==static_cast<int>(HcalForward)) {
41  etaMn = hpar->etaMin[2]; etaMx = hpar->etaMax[2];
42  }
43  double eta = 0, deta = 0, phi = 0, dphi = 0, rz = 0, drz = 0;
44  bool ok = false, flagrz = true;
45  if ((idet==static_cast<int>(HcalBarrel)||idet==static_cast<int>(HcalEndcap)||
46  idet==static_cast<int>(HcalOuter)||idet==static_cast<int>(HcalForward))
47  && etaR >=etaMn && etaR <= etaMx && depth > 0) ok = true;
48  if (idet == static_cast<int>(HcalEndcap) && depth>(int)(hpar->zHE.size()))ok=false;
49  else if (idet == static_cast<int>(HcalBarrel) && depth > 17) ok=false;
50  else if (idet == static_cast<int>(HcalOuter) && depth != 4) ok=false;
51  else if (idet == static_cast<int>(HcalForward) && depth > maxDepth[2]) ok=false;
52  if (ok) {
53  eta = getEta(idet, etaR, zside, depth);
54  deta = deltaEta(idet, etaR, depth);
55  double fibin, fioff;
56  if (idet == static_cast<int>(HcalBarrel)||
57  idet == static_cast<int>(HcalOuter)) {
58  fioff = hpar->phioff[0];
59  fibin = hpar->phibin[etaR-1];
60  } else if (idet == static_cast<int>(HcalEndcap)) {
61  fioff = hpar->phioff[1];
62  fibin = hpar->phibin[etaR-1];
63  } else {
64  fioff = hpar->phioff[2];
65  fibin = hpar->phitable[etaR-hpar->etaMin[2]];
66  if (unitPhi(fibin) > 2) fioff = hpar->phioff[4];
67  }
68  phi = fioff + (iphi - 0.5)*fibin;
69  dphi = 0.5*fibin;
70  if (idet == static_cast<int>(HcalForward)) {
71  int ir = nR + hpar->etaMin[2] - etaR - 1;
72  if (ir > 0 && ir < nR) {
73  rz = 0.5*(hpar->rTable[ir]+hpar->rTable[ir-1]);
74  drz = 0.5*(hpar->rTable[ir]-hpar->rTable[ir-1]);
75  } else {
76  ok = false;
77 #ifdef DebugLog
78  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: wrong eta " << etaR
79  << " (" << ir << "/" << nR << ") Detector "
80  << idet;
81 #endif
82  }
83  } else if (etaR <= nEta) {
84  int laymin(depth), laymax(depth);
85  if (idet == static_cast<int>(HcalOuter)) {
86  laymin = (etaR > hpar->noff[2]) ? ((int)(hpar->zHE.size())) : ((int)(hpar->zHE.size()))-1;
87  laymax = ((int)(hpar->zHE.size()));
88  }
89  double d1=0, d2=0;
90  if (idet == static_cast<int>(HcalEndcap)) {
91  flagrz = false;
92  d1 = hpar->zHE[laymin-1] - hpar->dzHE[laymin-1];
93  d2 = hpar->zHE[laymax-1] + hpar->dzHE[laymax-1];
94  } else {
95  d1 = hpar->rHB[laymin-1] - hpar->drHB[laymin-1];
96  d2 = hpar->rHB[laymax-1] + hpar->drHB[laymax-1];
97  }
98  rz = 0.5*(d2+d1);
99  drz = 0.5*(d2-d1);
100  } else {
101  ok = false;
102 #ifdef DebugLog
103  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: wrong depth " << depth
104  << " or etaR " << etaR << " for detector "
105  << idet;
106 #endif
107  }
108  } else {
109  ok = false;
110 #ifdef DebugLog
111  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: wrong depth " << depth
112  << " det " << idet;
113 #endif
114  }
115  HcalCellType::HcalCell tmp(ok,eta,deta,phi,dphi,rz,drz,flagrz);
116 
117 #ifdef DebugLog
118  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: det/side/depth/etaR/phi "
119  << idet << "/" << zside << "/" << depth << "/"
120  << etaR << "/" << iphi << " Cell Flag " << tmp.ok
121  << " " << tmp.eta << " " << tmp.deta << " phi "
122  << tmp.phi << " " << tmp.dphi << " r(z) " << tmp.rz
123  << " " << tmp.drz << " " << tmp.flagrz;
124 #endif
125  return tmp;
126 }
127 
128 std::vector<std::pair<double,double> > HcalDDDSimConstants::getConstHBHE(const int type) const {
129 
130  std::vector<std::pair<double,double> > gcons;
131  if (type == 0) {
132  for (unsigned int i=0; i<hpar->rHB.size(); ++i) {
133  gcons.push_back(std::pair<double,double>(hpar->rHB[i],hpar->drHB[i]));
134  }
135  } else {
136  for (unsigned int i=0; i<hpar->zHE.size(); ++i) {
137  gcons.push_back(std::pair<double,double>(hpar->zHE[i],hpar->dzHE[i]));
138  }
139  }
140  return gcons;
141 }
142 
143 
144 std::pair<int,double> HcalDDDSimConstants::getDetEta(double eta, int depth) {
145 
146  int hsubdet(0), ieta(0);
147  double etaR(0);
148  double heta = fabs(eta);
149  for (int i = 0; i < nEta; i++)
150  if (heta > hpar->etaTable[i]) ieta = i + 1;
151  if (heta <= hpar->etaRange[1]) {
152  if ((ieta <= hpar->etaMin[1] && depth==3) || ieta > hpar->etaMax[0]) {
153  hsubdet = static_cast<int>(HcalEndcap);
154  } else {
155  hsubdet = static_cast<int>(HcalBarrel);
156  }
157  etaR = eta;
158  } else {
159  hsubdet = static_cast<int>(HcalForward);
160  double theta = 2.*atan(exp(-heta));
161  double hR = zVcal*tan(theta);
162  etaR = (eta >= 0. ? hR : -hR);
163  }
164  return std::pair<int,double>(hsubdet,etaR);
165 }
166 
167 int HcalDDDSimConstants::getEta(int det, int lay, double hetaR) {
168 
169  int ieta(0);
170  if (det == static_cast<int>(HcalForward)) { // Forward HCal
171  ieta = hpar->etaMax[2];
172  for (int i = nR-1; i > 0; i--)
173  if (hetaR < hpar->rTable[i]) ieta = hpar->etaMin[2] + nR - i - 1;
174  } else { // Barrel or Endcap
175  ieta = 1;
176  for (int i = 0; i < nEta-1; i++)
177  if (hetaR > hpar->etaTable[i]) ieta = i + 1;
178  if (det == static_cast<int>(HcalBarrel)) {
179  if (ieta > hpar->etaMax[0]) ieta = hpar->etaMax[0];
180  if (lay == 18) {
181  if (hetaR > etaHO[1] && ieta == hpar->noff[2]) ieta++;
182  }
183  } else if (det == static_cast<int>(HcalEndcap)) {
184  if (ieta <= hpar->etaMin[1]) ieta = hpar->etaMin[1];
185  }
186  }
187  return ieta;
188 }
189 
190 std::pair<int,int> HcalDDDSimConstants::getEtaDepth(int det, int etaR, int phi,
191  int depth, int lay) {
192 
193  //Modify the depth index
194  if (det == static_cast<int>(HcalForward)) { // Forward HCal
195  } else if (det == static_cast<int>(HcalOuter)) {
196  depth = 4;
197  } else {
198  if (lay >= 0) {
199  depth= layerGroup( etaR-1, lay-1 );
200  if (etaR == hpar->noff[0] && lay > 1) {
201  int kphi = phi + int((hpar->phioff[3]+0.1)/hpar->phibin[etaR-1]);
202  kphi = (kphi-1)%4 + 1;
203  if (kphi == 2 || kphi == 3) depth = layerGroup( etaR-1, lay-2 );
204  }
205  } else if (det == static_cast<int>(HcalBarrel)) {
206  if (depth>maxDepth[0]) depth = maxDepth[0];
207  }
208  if (etaR >= hpar->noff[1] && depth > depthEta29[0]) {
209  etaR -= depthEta29[1];
210  } else if (etaR == hpar->etaMin[1]) {
211  if (det == static_cast<int>(HcalBarrel)) {
212  if (depth > depthEta16[0]) depth = depthEta16[0];
213  } else {
214  if (depth < depthEta16[1]) depth = depthEta16[1];
215  }
216  }
217  }
218  return std::pair<int,int>(etaR,depth);
219 }
220 
221 double HcalDDDSimConstants::getEtaHO(double& etaR, double& x, double& y,
222  double& z) const {
223 
224  if (hpar->zHO.size() > 4) {
225  double eta = fabs(etaR);
226  double r = std::sqrt(x*x+y*y);
227  if (r > rminHO) {
228  double zz = fabs(z);
229  if (zz > hpar->zHO[3]) {
230  if (eta <= hpar->etaTable[10]) eta = hpar->etaTable[10]+0.001;
231  } else if (zz > hpar->zHO[1]) {
232  if (eta <= hpar->etaTable[4]) eta = hpar->etaTable[4]+0.001;
233  }
234  }
235  eta = (z >= 0. ? eta : -eta);
236 #ifdef DebugLog
237  edm::LogInfo ("HCalGeom") << "R " << r << " Z " << z << " eta " << etaR
238  << ":" << eta;
239  if (eta != etaR) edm::LogInfo ("HCalGeom") << "**** Check *****";
240 #endif
241  return eta;
242  } else {
243  return etaR;
244  }
245 }
246 
247 unsigned int HcalDDDSimConstants::findLayer(int layer, const std::vector<HcalParameters::LayerItem>& layerGroup) const {
248 
249  unsigned int id = layerGroup.size();
250  for (unsigned int i = 0; i < layerGroup.size(); i++) {
251  if (layer == (int)(layerGroup[i].layer)) {
252  id = i;
253  break;
254  }
255  }
256  return id;
257 }
258 
259 std::pair<int,int> HcalDDDSimConstants::getModHalfHBHE(const int type) const {
260 
261  if (type == 0) {
262  return std::pair<int,int>(nmodHB,nzHB);
263  } else {
264  return std::pair<int,int>(nmodHE,nzHE);
265  }
266 }
267 
268 std::pair<double,double> HcalDDDSimConstants::getPhiCons(int det, int ieta) {
269 
270  double fioff(0), fibin(0);
271  if (det == static_cast<int>(HcalForward)) { // Forward HCal
272  fioff = hpar->phioff[2];
273  fibin = hpar->phitable[ieta-hpar->etaMin[2]];
274  if (unitPhi(fibin) > 2) { // HF double-phi
275  fioff = hpar->phioff[4];
276  }
277  } else { // Barrel or Endcap
278  if (det == static_cast<int>(HcalBarrel)) {
279  fioff = hpar->phioff[0];
280  } else {
281  fioff = hpar->phioff[1];
282  }
283  fibin = hpar->phibin[ieta-1];
284  }
285  return std::pair<double,double>(fioff,fibin);
286 }
287 
288 std::vector<HcalCellType> HcalDDDSimConstants::HcalCellTypes() const{
289 
290  std::vector<HcalCellType> cellTypes = HcalCellTypes(HcalBarrel);
291 #ifdef DebugLog
292  edm::LogInfo ("HCalGeom") << "HcalDDDSimConstants: " << cellTypes.size()
293  << " cells of type HCal Barrel";
294  for (unsigned int i=0; i<cellTypes.size(); i++)
295  edm::LogInfo ("HCalGeom") << "Cell " << i << " " << cellTypes[i];
296 #endif
297 
298  std::vector<HcalCellType> hoCells = HcalCellTypes(HcalOuter);
299 #ifdef DebugLog
300  edm::LogInfo ("HCalGeom") << "HcalDDDSimConstants: " << hoCells.size()
301  << " cells of type HCal Outer";
302  for (unsigned int i=0; i<hoCells.size(); i++)
303  edm::LogInfo ("HCalGeom") << "Cell " << i << " " << hoCells[i];
304 #endif
305  cellTypes.insert(cellTypes.end(), hoCells.begin(), hoCells.end());
306 
307  std::vector<HcalCellType> heCells = HcalCellTypes(HcalEndcap);
308 #ifdef DebugLog
309  edm::LogInfo ("HCalGeom") << "HcalDDDSimConstants: " << heCells.size()
310  << " cells of type HCal Endcap";
311  for (unsigned int i=0; i<heCells.size(); i++)
312  edm::LogInfo ("HCalGeom") << "Cell " << i << " " << heCells[i];
313 #endif
314  cellTypes.insert(cellTypes.end(), heCells.begin(), heCells.end());
315 
316  std::vector<HcalCellType> hfCells = HcalCellTypes(HcalForward);
317 #ifdef DebugLog
318  edm::LogInfo ("HCalGeom") << "HcalDDDSimConstants: " << hfCells.size()
319  << " cells of type HCal Forward";
320  for (unsigned int i=0; i<hfCells.size(); i++)
321  edm::LogInfo ("HCalGeom") << "Cell " << i << " " << hfCells[i];
322 #endif
323  cellTypes.insert(cellTypes.end(), hfCells.begin(), hfCells.end());
324 
325  return cellTypes;
326 }
327 
328 std::vector<HcalCellType> HcalDDDSimConstants::HcalCellTypes(HcalSubdetector subdet,
329  int ieta, int depthl) const {
330 
331  std::vector<HcalCellType> cellTypes;
332  if (subdet == HcalForward) {
333  if (dzVcal < 0) return cellTypes;
334  }
335 
336  int dmin, dmax, indx, nz, nmod;
337  double hsize = 0;
338  switch(subdet) {
339  case HcalEndcap:
340  dmin = 1; dmax = 19; indx = 1; nz = nzHE; nmod = nmodHE;
341  break;
342  case HcalForward:
343  dmin = 1; dmax = maxDepth[2]; indx = 2; nz = 2; nmod = 18;
344  break;
345  case HcalOuter:
346  dmin = 4; dmax = 4; indx = 0; nz = nzHB; nmod = nmodHB;
347  break;
348  default:
349  dmin = 1; dmax = 17; indx = 0; nz = nzHB; nmod = nmodHB;
350  break;
351  }
352  if (depthl > 0) dmin = dmax = depthl;
353  int ietamin = (ieta>0) ? ieta : hpar->etaMin[indx];
354  int ietamax = (ieta>0) ? ieta : hpar->etaMax[indx];
355  int phi = 1, zside = 1;
356 
357  // Get the Cells
358  int subdet0 = static_cast<int>(subdet);
359  for (int depth=dmin; depth<=dmax; depth++) {
360  int shift = getShift(subdet, depth);
361  double gain = getGain (subdet, depth);
362  if (subdet == HcalForward) {
363  if (depth%2 == 1) hsize = dzVcal;
364  else hsize = dzVcal-0.5*dlShort;
365  }
366  for (int eta=ietamin; eta<= ietamax; eta++) {
367  HcalCellType::HcalCell temp1 = cell(subdet0,zside,depth,eta,phi);
368  if (temp1.ok) {
369  int units = unitPhi (subdet0, eta);
370  HcalCellType temp2(subdet, eta, phi, depth, temp1,
371  shift, gain, nz, nmod, hsize, units);
372  if (subdet == HcalOuter) {
373  if (eta == hpar->noff[4]) {
374  std::vector<int> missPlus, missMinus;
375  int kk = 7;
376  for (int miss=0; miss<hpar->noff[5]; miss++) {
377  missPlus.push_back(hpar->noff[kk]);
378  kk++;
379  }
380  for (int miss=0; miss<hpar->noff[6]; miss++) {
381  missMinus.push_back(hpar->noff[kk]);
382  kk++;
383  }
384  temp2.setMissingPhi(missPlus, missMinus);
385  }
386  } else if (subdet == HcalForward) {
387  if (idHF2QIE.size() > 0 && depth > 2) {
388  int phiInc = 4/temp2.nPhiModule();
389  int iphi = 1;
390  if (temp2.unitPhi() == 4) iphi = 3;
391  std::vector<int> missPlus, missMinus;
392  for (int k=0; k<temp2.nPhiBins(); ++k) {
393  bool okp(true), okn(true);
394  for (unsigned int l=0; l<idHF2QIE.size(); ++l) {
395  if ((eta == idHF2QIE[l].ietaAbs()) &&
396  (iphi == idHF2QIE[l].iphi())) {
397  if (idHF2QIE[l].zside() == 1) okp = false;
398  else okn = false;
399  }
400  }
401  if (okp) missPlus.push_back(iphi);
402  if (okn) missMinus.push_back(iphi);
403  iphi += phiInc;
404  }
405  temp2.setMissingPhi(missPlus, missMinus);
406  }
407  }
408  cellTypes.push_back(temp2);
409  }
410  }
411  }
412  return cellTypes;
413 }
414 
415 int HcalDDDSimConstants::maxHFDepth(int eta, int iphi) const {
416  int mxdepth = maxDepth[2];
417  if (idHF2QIE.size() > 0) {
418  bool ok(false);
419  for (unsigned int k=0; k<idHF2QIE.size(); ++k) {
420  if ((eta == idHF2QIE[k].ieta()) && (iphi == idHF2QIE[k].iphi())) {
421  ok = true; break;
422  }
423  }
424  if (!ok) mxdepth = 2;
425  }
426  return mxdepth;
427 }
428 
430 
431  unsigned int num = 0;
432  std::vector<HcalCellType> cellTypes = HcalCellTypes(subdet);
433  for (unsigned int i=0; i<cellTypes.size(); i++) {
434  int ncur = (cellTypes[i].nHalves() > 1) ? 2*cellTypes[i].nPhiBins() :
435  cellTypes[i].nPhiBins();
436  num += (unsigned int)(ncur);
437  num -= (unsigned int)(cellTypes[i].nPhiMissingBins());
438  }
439 #ifdef DebugLog
440  edm::LogInfo ("HCalGeom") << "HcalDDDSimConstants:numberOfCells "
441  << cellTypes.size() << " " << num
442  << " for subdetector " << subdet;
443 #endif
444  return num;
445 }
446 
448 
449  int iphi_skip = phi;
450  if (units==2) iphi_skip = (phi-1)*2+1;
451  else if (units==4) iphi_skip = (phi-1)*4-1;
452  if (iphi_skip < 0) iphi_skip += 72;
453  return iphi_skip;
454 }
455 
457 
458  std::cout << "Tile Information for HB from " << hpar->etaMin[0] << " to " << hpar->etaMax[0] << "\n\n";
459  for (int eta=hpar->etaMin[0]; eta<= hpar->etaMax[0]; eta++) {
460  int dmax = 1;
461  if (depths[0][eta-1] < 17) dmax = 2;
462  if (maxDepth[0] > 2) {
463  if (eta == hpar->etaMin[1]) dmax = depthEta16[0];
464  else dmax = maxDepth[0];
465  }
466  for (int depth=1; depth<=dmax; depth++)
468  }
469 
470  std::cout << "\nTile Information for HE from " << hpar->etaMin[1] << " to " << hpar->etaMax[1] << "\n\n";
471  for (int eta=hpar->etaMin[1]; eta<= hpar->etaMax[1]; eta++) {
472  int dmin=1, dmax=3;
473  if (eta == hpar->etaMin[1]) {
474  dmax = dmin = depthEta16[1];
475  } else if (depths[0][eta-1] > 18) {
476  dmax = 1;
477  } else if (depths[1][eta-1] > 18) {
478  dmax = 2;
479  }
480  for (int depth=dmin; depth<=dmax; depth++)
482  }
483 }
484 
485 int HcalDDDSimConstants::unitPhi(int det, int etaR) const {
486 
487  double dphi = (det == static_cast<int>(HcalForward)) ? hpar->phitable[etaR-hpar->etaMin[2]] : hpar->phibin[etaR-1];
488  return unitPhi(dphi);
489 }
490 
491 int HcalDDDSimConstants::unitPhi(double dphi) const {
492 
493  const double fiveDegInRad = 2*M_PI/72;
494  int units = int(dphi/fiveDegInRad+0.5);
495  return units;
496 }
497 
499 
500  nEta = hpar->etaTable.size();
501  nR = hpar->rTable.size();
502  nPhiF = nR - 1;
503 
504 #ifdef DebugLog
505  for (int i=0; i<nEta-1; ++i) {
506  std::cout << "HcalDDDSimConstants:Read LayerGroup" << i << ":";
507  for (unsigned int k=0; k<layerGroupSize( i ); k++)
508  std::cout << " [" << k << "] = " << layerGroup( i, k );
509  std::cout << std::endl;
510  }
511 #endif
512 
513  // Geometry parameters for HF
514  dlShort = hpar->gparHF[0];
515  zVcal = hpar->gparHF[4];
516  dzVcal = hpar->dzVcal;
517 #ifdef DebugLog
518  std::cout << "HcalDDDSimConstants: dlShort " << dlShort << " zVcal " << zVcal
519  << " and dzVcal " << dzVcal << std::endl;
520 #endif
521 
522  //Transform some of the parameters
524  maxDepth[0] = maxDepth[1] = 0;
525  for (int i=0; i<nEta-1; ++i) {
526  unsigned int imx = layerGroupSize( i );
527  int laymax = (imx > 0) ? layerGroup( i, imx-1 ) : 0;
528  if (i < hpar->etaMax[0]) {
529  int laymax0 = (imx > 16) ? layerGroup( i, 16 ) : laymax;
530  if (i+1 == hpar->etaMax[0] && laymax0 > 2) laymax0 = 2;
531  if (maxDepth[0] < laymax0) maxDepth[0] = laymax0;
532  }
533  if (i >= hpar->etaMin[1]-1 && i < hpar->etaMax[1]) {
534  if (maxDepth[1] < laymax) maxDepth[1] = laymax;
535  }
536  }
537 #ifdef DebugLog
538  for (int i=0; i<4; ++i)
539  std::cout << "Detector Type [" << i << "] iEta " << hpar->etaMin[i] << ":"
540  << hpar->etaMax[i] << " MaxDepth " << maxDepth[i] << std::endl;
541 #endif
542 
543  int maxdepth = (maxDepth[1]>maxDepth[0]) ? maxDepth[1] : maxDepth[0];
544  for (int i=0; i<maxdepth; ++i) {
545  for (int k=0; k<nEta-1; ++k) {
546  int layermx = ((k+1 < hpar->etaMin[1]) && i < maxDepth[0]) ? 17 : (int)layerGroupSize( k );
547  int ll = layermx;
548  for (int l=layermx-1; l >= 0; --l) {
549  if ((int)layerGroup( k, l ) == i+1) {
550  ll = l+1; break;
551  }
552  }
553  depths[i].push_back(ll);
554  }
555 
556 #ifdef DebugLog
557  std::cout << "Depth " << i << " with " << depths[i].size() << " etas:";
558  for (int k=0; k<nEta-1; ++k) std::cout << " " << depths[i][k];
559  std::cout << std::endl;
560 #endif
561  }
562 
563  nzHB = hpar->modHB[0];
564  nmodHB = hpar->modHB[1];
565  nzHE = hpar->modHE[0];
566  nmodHE = hpar->modHE[1];
567 #ifdef DebugLog
568  std::cout << "HcalDDDSimConstants:: " << nzHB << ":" << nmodHB
569  << " barrel and " << nzHE << ":" << nmodHE
570  << " endcap half-sectors" << std::endl;
571 #endif
572 
573  if (hpar->rHB.size() > 17 && hpar->zHO.size() > 4) {
574  rminHO = hpar->rHO[0];
575  for (int k=0; k<4; ++k) etaHO[k] = hpar->rHO[k+1];
576  } else {
577  rminHO =-1.0;
578  etaHO[0] = hpar->etaTable[4];
579  etaHO[1] = hpar->etaTable[4];
580  etaHO[2] = hpar->etaTable[10];
581  etaHO[3] = hpar->etaTable[10];
582  }
583 #ifdef DebugLog
584  std::cout << "HO Eta boundaries " << etaHO[0] << " " << etaHO[1]
585  << " " << etaHO[2] << " " << etaHO[3] << std::endl;
586  std::cout << "HO Parameters " << rminHO << " " << hpar->zHO.size();
587  for (int i=0; i<4; ++i) std::cout << " eta[" << i << "] = " << etaHO[i];
588  for (unsigned int i=0; i<hpar->zHO.size(); ++i)
589  std::cout << " zHO[" << i << "] = " << hpar->zHO[i];
590  std::cout << std::endl;
591 #endif
592 
593  int noffsize = 7 + hpar->noff[5] + hpar->noff[6];
594  if ((int)(hpar->noff.size()) > (noffsize+3)) {
595  depthEta16[0] = hpar->noff[noffsize];
596  depthEta16[1] = hpar->noff[noffsize+1];
597  depthEta29[0] = hpar->noff[noffsize+2];
598  depthEta29[1] = hpar->noff[noffsize+3];
599  } else {
600  depthEta16[0] = 2;
601  depthEta16[1] = 3;
602  depthEta29[0] = 2;
603  depthEta29[1] = 1;
604  }
605 #ifdef DebugLog
606  std::cout << "Depth index at ieta = 16 for HB (max) " << depthEta16[0]
607  << " HE (min) " << depthEta16[1] << "; max depth for itea = 29 : ("
608  << depthEta29[0] << ":" << depthEta29[1] << ")" << std::endl;
609 #endif
610 
611  if ((int)(hpar->noff.size()) > (noffsize+4)) {
612  int npair = hpar->noff[noffsize+4];
613  int kk = noffsize+4;
614  for (int k=0; k<npair; ++k) {
615  idHF2QIE.push_back(HcalDetId(HcalForward,hpar->noff[kk+1],hpar->noff[kk+2],1));
616  kk += 2;
617  }
618  }
619 #ifdef DebugLog
620  std::cout << idHF2QIE.size() << " detector channels having 2 QIE cards:"
621  << std::endl;
622  for (unsigned int k=0; k<idHF2QIE.size(); ++k)
623  std::cout << " [" << k << "] " << idHF2QIE[k] << std::endl;
624 #endif
625 }
626 
627 double HcalDDDSimConstants::deltaEta(int det, int etaR, int depth) const {
628 
629  double tmp = 0;
630  if (det == static_cast<int>(HcalForward)) {
631  int ir = nR + hpar->etaMin[2] - etaR - 1;
632  if (ir > 0 && ir < nR) {
633  double z = zVcal;
634  if (depth%2 != 1) z += dlShort;
635  tmp = 0.5*(getEta(hpar->rTable[ir-1],z)-getEta(hpar->rTable[ir],z));
636  }
637  } else {
638  if (etaR > 0 && etaR < nEta) {
639  if (etaR == hpar->noff[1]-1 && depth > 2) {
640  tmp = 0.5*(hpar->etaTable[etaR+1]-hpar->etaTable[etaR-1]);
641  } else if (det == static_cast<int>(HcalOuter)) {
642  if (etaR == hpar->noff[2]) {
643  tmp = 0.5*(etaHO[0]-hpar->etaTable[etaR-1]);
644  } else if (etaR == hpar->noff[2]+1) {
645  tmp = 0.5*(hpar->etaTable[etaR]-etaHO[1]);
646  } else if (etaR == hpar->noff[3]) {
647  tmp = 0.5*(etaHO[2]-hpar->etaTable[etaR-1]);
648  } else if (etaR == hpar->noff[3]+1) {
649  tmp = 0.5*(hpar->etaTable[etaR]-etaHO[3]);
650  } else {
651  tmp = 0.5*(hpar->etaTable[etaR]-hpar->etaTable[etaR-1]);
652  }
653  } else {
654  tmp = 0.5*(hpar->etaTable[etaR]-hpar->etaTable[etaR-1]);
655  }
656  }
657  }
658 #ifdef DebugLog
659  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants::deltaEta " << etaR << " "
660  << depth << " ==> " << tmp;
661 #endif
662  return tmp;
663 }
664 
665 double HcalDDDSimConstants::getEta(int det, int etaR, int zside,
666  int depth) const {
667 
668  double tmp = 0;
669  if (det == static_cast<int>(HcalForward)) {
670  int ir = nR + hpar->etaMin[2] - etaR - 1;
671  if (ir > 0 && ir < nR) {
672  double z = zVcal;
673  if (depth%2 != 1) z += dlShort;
674  tmp = 0.5*(getEta(hpar->rTable[ir-1],z)+getEta(hpar->rTable[ir],z));
675  }
676  } else {
677  if (etaR > 0 && etaR < nEta) {
678  if (etaR == hpar->noff[1]-1 && depth > 2) {
679  tmp = 0.5*(hpar->etaTable[etaR+1]+hpar->etaTable[etaR-1]);
680  } else if (det == static_cast<int>(HcalOuter)) {
681  if (etaR == hpar->noff[2]) {
682  tmp = 0.5*(etaHO[0]+hpar->etaTable[etaR-1]);
683  } else if (etaR == hpar->noff[2]+1) {
684  tmp = 0.5*(hpar->etaTable[etaR]+etaHO[1]);
685  } else if (etaR == hpar->noff[3]) {
686  tmp = 0.5*(etaHO[2]+hpar->etaTable[etaR-1]);
687  } else if (etaR == hpar->noff[3]+1) {
688  tmp = 0.5*(hpar->etaTable[etaR]+etaHO[3]);
689  } else {
690  tmp = 0.5*(hpar->etaTable[etaR]+hpar->etaTable[etaR-1]);
691  }
692  } else {
693  tmp = 0.5*(hpar->etaTable[etaR]+hpar->etaTable[etaR-1]);
694  }
695  }
696  }
697  if (zside == 0) tmp = -tmp;
698 #ifdef DebugLog
699  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants::getEta " << etaR << " "
700  << zside << " " << depth << " ==> " << tmp;
701 #endif
702  return tmp;
703 }
704 
705 double HcalDDDSimConstants::getEta(double r, double z) const {
706 
707  double tmp = 0;
708  if (z != 0) tmp = -log(tan(0.5*atan(r/z)));
709 #ifdef DebugLog
710  edm::LogInfo("HCalGeom") << "HcalDDDSimConstants::getEta " << r << " " << z
711  << " ==> " << tmp;
712 #endif
713  return tmp;
714 }
715 
717 
718  int shift;
719  switch(subdet) {
720  case HcalEndcap:
721  shift = hpar->HEShift[0];
722  break;
723  case HcalForward:
724  shift = hpar->HFShift[(depth-1)%2];
725  break;
726  case HcalOuter:
727  shift = hpar->HBShift[3];
728  break;
729  default:
730  shift = hpar->HBShift[0];
731  break;
732  }
733  return shift;
734 }
735 
737 
738  double gain;
739  switch(subdet) {
740  case HcalEndcap:
741  gain = hpar->HEGains[0];
742  break;
743  case HcalForward:
744  gain = hpar->HFGains[(depth-1)%2];
745  break;
746  case HcalOuter:
747  gain = hpar->HBGains[3];
748  break;
749  default:
750  gain = hpar->HBGains[0];
751  break;
752  }
753  return gain;
754 }
755 
757  std::cout << "HcalDDDSimConstants::printTileHB for eta " << eta << " and depth " << depth << "\n";
758 
759  double etaL = hpar->etaTable.at(eta-1);
760  double thetaL = 2.*atan(exp(-etaL));
761  double etaH = hpar->etaTable.at(eta);
762  double thetaH = 2.*atan(exp(-etaH));
763  int layL=0, layH=0;
764  if (depth == 1) {
765  layH = depths[0][eta-1];
766  } else {
767  layL = depths[0][eta-1];
768  layH = depths[1][eta-1];
769  }
770  std::cout << "\ntileHB:: eta|depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << "\n";
771  for (int lay=layL; lay<layH; ++lay) {
772  std::vector<double> area(2,0);
773  int kk=0;
774  for (unsigned int k=0; k<hpar->layHB.size(); ++k) {
775  if (lay == hpar->layHB[k]) {
776  double zmin = hpar->rhoxHB[k]*std::cos(thetaL)/std::sin(thetaL);
777  double zmax = hpar->rhoxHB[k]*std::cos(thetaH)/std::sin(thetaH);
778  double dz = (std::min(zmax,hpar->dxHB[k]) - zmin);
779  if (dz > 0) {
780  area[kk] = dz*hpar->dyHB[k];
781  kk++;
782  }
783  }
784  }
785  if (area[0] > 0) std::cout << std::setw(2) << lay << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
786  }
787 }
788 
790 
791  double etaL = hpar->etaTable[eta-1];
792  double thetaL = 2.*atan(exp(-etaL));
793  double etaH = hpar->etaTable[eta];
794  double thetaH = 2.*atan(exp(-etaH));
795  int layL=0, layH=0;
796  if (eta == 16) {
797  layH = depths[2][eta-1];
798  } else if (depth == 1) {
799  layH = depths[0][eta-1];
800  } else if (depth == 2) {
801  layL = depths[0][eta-1];
802  layH = depths[1][eta-1];
803  } else {
804  layL = depths[1][eta-1];
805  layH = depths[2][eta-1];
806  }
807  double phib = hpar->phibin[eta-1];
808  int nphi = 2;
809  if (phib > 6*CLHEP::deg) nphi = 1;
810  std::cout << "\ntileHE:: Eta/depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << " phi " << nphi << "\n";
811  for (int lay=layL; lay<layH; ++lay) {
812  std::vector<double> area(4,0);
813  int kk=0;
814  for (unsigned int k=0; k<hpar->layHE.size(); ++k) {
815  if (lay == hpar->layHE[k]) {
816  double rmin = hpar->zxHE[k]*std::tan(thetaH);
817  double rmax = hpar->zxHE[k]*std::tan(thetaL);
818  if ((lay != 0 || eta == 18) &&
819  (lay != 1 || (eta == 18 && hpar->rhoxHE[k]-hpar->dyHE[k] > 1000) ||
820  (eta != 18 && hpar->rhoxHE[k]-hpar->dyHE[k] < 1000)) &&
821  rmin+30 < hpar->rhoxHE[k]+hpar->dyHE[k] &&
822  rmax > hpar->rhoxHE[k]-hpar->dyHE[k]) {
823  rmin = std::max(rmin,hpar->rhoxHE[k]-hpar->dyHE[k]);
824  rmax = std::min(rmax,hpar->rhoxHE[k]+hpar->dyHE[k]);
825  double dx1 = rmin*std::tan(phib);
826  double dx2 = rmax*std::tan(phib);
827  double ar1=0, ar2=0;
828  if (nphi == 1) {
829  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-4.*hpar->dx1HE[k]);
830  } else {
831  ar1 = 0.5*(rmax-rmin)*(dx1+dx2-2.*hpar->dx1HE[k]);
832  ar2 = 0.5*(rmax-rmin)*((rmax+rmin)*tan(10.*CLHEP::deg)-4*hpar->dx1HE[k])-ar1;
833  }
834  area[kk] = ar1;
835  area[kk+2] = ar2;
836  kk++;
837  }
838  }
839  }
840  if (area[0] > 0 && area[1] > 0) {
841  int lay0 = lay-1;
842  if (eta == 18) lay0++;
843  if (nphi == 1) {
844  std::cout << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
845  } else {
846  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";
847  }
848  }
849  }
850 }
851 
852 unsigned int HcalDDDSimConstants::layerGroupSize(unsigned int eta) const {
853  unsigned int k = 0;
854  for( auto const & it : hpar->layerGroupEtaSim ) {
855  if( it.layer == eta + 1 ) {
856  return it.layerGroup.size();
857  }
858  if( it.layer > eta + 1 ) break;
859  k = it.layerGroup.size();
860  }
861  return k;
862 }
863 
864 unsigned int HcalDDDSimConstants::layerGroup(unsigned int eta,
865  unsigned int i) const {
866  unsigned int k = 0;
867  for( auto const & it : hpar->layerGroupEtaSim ) {
868  if( it.layer == eta + 1 ) {
869  return it.layerGroup.at( i );
870  }
871  if( it.layer > eta + 1 ) break;
872  k = it.layerGroup.at( i );
873  }
874  return k;
875 }
std::pair< int, int > getEtaDepth(int det, int etaR, int phi, int depth, int lay)
type
Definition: HCALResponse.h:21
std::vector< double > zHO
int i
Definition: DBlmapReader.cc:9
std::vector< double > etaTable
void setMissingPhi(std::vector< int > &, std::vector< int > &)
Definition: HcalCellType.cc:81
int unitPhi(int det, int etaR) const
void printTileHB(int eta, int depth) const
std::vector< int > layHE
std::vector< int > maxDepth
std::vector< double > rHB
int getEta(int det, int lay, double hetaR)
std::vector< double > rhoxHB
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::vector< double > HBGains
Geom::Theta< T > theta() const
std::vector< int > HEShift
int zside(DetId const &)
std::vector< std::pair< double, double > > getConstHBHE(const int type) const
HcalDDDSimConstants(const HcalParameters *hp)
unsigned int numberOfCells(HcalSubdetector) const
float float float z
std::vector< int > etaMax
std::vector< int > depths[nDepthMax]
std::vector< int > modHB
std::vector< double > rhoxHE
std::vector< double > zxHE
T x() const
Cartesian x coordinate.
std::vector< double > zHE
double getEtaHO(double &etaR, double &x, double &y, double &z) const
std::vector< LayerItem > layerGroupEtaSim
void printTileHE(int eta, int depth) const
std::vector< double > dyHE
std::vector< double > dx1HE
std::vector< double > rHO
susybsm::HSCParticleRefProd hp
Definition: classes.h:27
T sqrt(T t)
Definition: SSEVec.h:18
std::vector< double > dzHE
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
int getShift(HcalSubdetector subdet, int depth) const
double getGain(HcalSubdetector subdet, int depth) const
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
HcalSubdetector
Definition: HcalAssistant.h:31
std::vector< HcalCellType > HcalCellTypes() const
std::vector< int > modHE
std::vector< HcalDetId > idHF2QIE
T min(T a, T b)
Definition: MathUtil.h:58
std::pair< int, int > getModHalfHBHE(const int type) const
std::vector< int > HFShift
unsigned int findLayer(int layer, const std::vector< HcalParameters::LayerItem > &layerGroup) const
#define M_PI
std::vector< double > HEGains
std::pair< int, double > getDetEta(double eta, int depth)
std::vector< double > dyHB
std::vector< double > phioff
std::pair< double, double > getPhiCons(int det, int ieta)
std::vector< double > gparHF
double deltaEta(int det, int eta, int depth) const
Geom::Phi< T > phi() const
std::vector< double > dxHB
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
std::vector< double > rTable
std::vector< double > phitable
TString units(TString variable, Char_t axis)
int nPhiModule() const
Definition: HcalCellType.h:52
std::vector< double > phibin
if(dp >Float(M_PI)) dp-
HcalCellType::HcalCell cell(int det, int zside, int depth, int etaR, int iphi) const
int phiNumber(int phi, int unit) const
std::vector< double > drHB
std::vector< double > HFGains
static unsigned int const shift
std::vector< int > noff
tuple cout
Definition: gather_cfg.py:145
std::vector< int > HBShift
int nPhiBins() const
the number of these cells in a ring
Definition: HcalCellType.h:51
const HcalParameters * hpar
std::vector< int > maxDepth
int unitPhi() const
Definition: HcalCellType.h:59
std::vector< int > layHB
unsigned int layerGroup(unsigned int eta, unsigned int i) const
std::vector< int > etaMin
int maxHFDepth(int ieta, int iphi) const
unsigned int layerGroupSize(unsigned int eta) const