Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-17 23:36:22

0001 #include "Geometry/HcalCommonData/interface/HcalDDDSimConstants.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "DataFormats/Math/interface/GeantUnits.h"
0006 
0007 //#define EDM_ML_DEBUG
0008 using namespace geant_units::operators;
0009 
0010 HcalDDDSimConstants::HcalDDDSimConstants(const HcalParameters* hp) : hpar(hp) {
0011 #ifdef EDM_ML_DEBUG
0012   edm::LogInfo("HCalGeom") << "HcalDDDSimConstants::HcalDDDSimConstants (const HcalParameters* hp) constructor\n";
0013 #endif
0014 
0015   initialize();
0016 #ifdef EDM_ML_DEBUG
0017   std::vector<HcalCellType> cellTypes = HcalCellTypes();
0018   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: " << cellTypes.size() << " cells of type HCal (All)\n";
0019 #endif
0020 }
0021 
0022 HcalDDDSimConstants::~HcalDDDSimConstants() {
0023 #ifdef EDM_ML_DEBUG
0024   edm::LogInfo("HCalGeom") << "HcalDDDSimConstants::destructed!!!\n";
0025 #endif
0026 }
0027 
0028 HcalCellType::HcalCell HcalDDDSimConstants::cell(
0029     const int& idet, const int& zside, const int& depth, const int& etaR, const int& iphi) const {
0030   double etaMn = hpar->etaMin[0];
0031   double etaMx = hpar->etaMax[0];
0032   if (idet == static_cast<int>(HcalEndcap)) {
0033     etaMn = hpar->etaMin[1];
0034     etaMx = hpar->etaMax[1];
0035   } else if (idet == static_cast<int>(HcalForward)) {
0036     etaMn = hpar->etaMin[2];
0037     etaMx = hpar->etaMax[2];
0038   }
0039   double eta = 0, deta = 0, phi = 0, dphi = 0, rz = 0, drz = 0;
0040   bool ok = false, flagrz = true;
0041   if ((idet == static_cast<int>(HcalBarrel) || idet == static_cast<int>(HcalEndcap) ||
0042        idet == static_cast<int>(HcalOuter) || idet == static_cast<int>(HcalForward)) &&
0043       etaR >= etaMn && etaR <= etaMx && depth > 0)
0044     ok = true;
0045   if (idet == static_cast<int>(HcalEndcap) && depth > (int)(hpar->zHE.size()))
0046     ok = false;
0047   else if (idet == static_cast<int>(HcalBarrel) && depth > maxLayerHB_ + 1)
0048     ok = false;
0049   else if (idet == static_cast<int>(HcalOuter) && depth != 4)
0050     ok = false;
0051   else if (idet == static_cast<int>(HcalForward) && depth > maxDepth[2])
0052     ok = false;
0053   if (ok) {
0054     eta = getEta(idet, etaR, zside, depth);
0055     deta = deltaEta(idet, etaR, depth);
0056     double fibin, fioff;
0057     if (idet == static_cast<int>(HcalBarrel) || idet == static_cast<int>(HcalOuter)) {
0058       fioff = hpar->phioff[0];
0059       fibin = hpar->phibin[etaR - 1];
0060     } else if (idet == static_cast<int>(HcalEndcap)) {
0061       fioff = hpar->phioff[1];
0062       fibin = hpar->phibin[etaR - 1];
0063     } else {
0064       fioff = hpar->phioff[2];
0065       fibin = hpar->phitable[etaR - hpar->etaMin[2]];
0066       if (unitPhi(fibin) > 2)
0067         fioff = hpar->phioff[4];
0068     }
0069     phi = -fioff + (iphi - 0.5) * fibin;
0070     dphi = 0.5 * fibin;
0071     if (idet == static_cast<int>(HcalForward)) {
0072       int ir = nR + hpar->etaMin[2] - etaR - 1;
0073       if (ir > 0 && ir < nR) {
0074         rz = 0.5 * (hpar->rTable[ir] + hpar->rTable[ir - 1]);
0075         drz = 0.5 * (hpar->rTable[ir] - hpar->rTable[ir - 1]);
0076       } else {
0077         ok = false;
0078 #ifdef EDM_ML_DEBUG
0079         edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: wrong eta " << etaR << " (" << ir << "/" << nR
0080                                  << ") Detector " << idet << std::endl;
0081 #endif
0082       }
0083     } else if (etaR <= nEta) {
0084       int laymin(depth), laymax(depth);
0085       if (idet == static_cast<int>(HcalOuter)) {
0086         laymin = (etaR > hpar->noff[2]) ? ((int)(hpar->zHE.size())) : ((int)(hpar->zHE.size())) - 1;
0087         laymax = ((int)(hpar->zHE.size()));
0088       }
0089       double d1 = 0, d2 = 0;
0090       if (idet == static_cast<int>(HcalEndcap)) {
0091         flagrz = false;
0092         d1 = hpar->zHE[laymin - 1] - hpar->dzHE[laymin - 1];
0093         d2 = hpar->zHE[laymax - 1] + hpar->dzHE[laymax - 1];
0094       } else {
0095         d1 = hpar->rHB[laymin - 1] - hpar->drHB[laymin - 1];
0096         d2 = hpar->rHB[laymax - 1] + hpar->drHB[laymax - 1];
0097       }
0098       rz = 0.5 * (d2 + d1);
0099       drz = 0.5 * (d2 - d1);
0100     } else {
0101       ok = false;
0102       edm::LogWarning("HCalGeom") << "HcalDDDSimConstants: wrong depth " << depth << " or etaR " << etaR
0103                                   << " for detector " << idet;
0104     }
0105   } else {
0106     ok = false;
0107     edm::LogWarning("HCalGeom") << "HcalDDDSimConstants: wrong depth " << depth << " det " << idet;
0108   }
0109   HcalCellType::HcalCell tmp(ok, eta, deta, phi, dphi, rz, drz, flagrz);
0110 
0111 #ifdef EDM_ML_DEBUG
0112   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: det/side/depth/etaR/"
0113                                << "phi " << idet << "/" << zside << "/" << depth << "/" << etaR << "/" << iphi
0114                                << " Cell Flag " << tmp.ok << " " << tmp.eta << " " << tmp.deta << " phi " << tmp.phi
0115                                << " " << tmp.dphi << " r(z) " << tmp.rz << " " << tmp.drz << " " << tmp.flagrz;
0116 #endif
0117   return tmp;
0118 }
0119 
0120 int HcalDDDSimConstants::findDepth(
0121     const int& det, const int& eta, const int& phi, const int& zside, const int& lay) const {
0122   int depth = (ldmap_.isValid(det, phi, zside)) ? ldmap_.getDepth(det, eta, phi, zside, lay) : -1;
0123   return depth;
0124 }
0125 
0126 unsigned int HcalDDDSimConstants::findLayer(const int& layer,
0127                                             const std::vector<HcalParameters::LayerItem>& layerGroup) const {
0128   unsigned int id = layerGroup.size();
0129   for (unsigned int i = 0; i < layerGroup.size(); i++) {
0130     if (layer == (int)(layerGroup[i].layer)) {
0131       id = i;
0132       break;
0133     }
0134   }
0135   return id;
0136 }
0137 
0138 std::vector<std::pair<double, double> > HcalDDDSimConstants::getConstHBHE(const int& type) const {
0139   std::vector<std::pair<double, double> > gcons;
0140   if (type == 0) {
0141     for (unsigned int i = 0; i < hpar->rHB.size(); ++i) {
0142       gcons.emplace_back(std::pair<double, double>(hpar->rHB[i], hpar->drHB[i]));
0143     }
0144   } else {
0145     for (unsigned int i = 0; i < hpar->zHE.size(); ++i) {
0146       gcons.emplace_back(std::pair<double, double>(hpar->zHE[i], hpar->dzHE[i]));
0147     }
0148   }
0149   return gcons;
0150 }
0151 
0152 int HcalDDDSimConstants::getDepthEta16(const int& det, const int& phi, const int& zside) const {
0153   int depth = ldmap_.getDepth16(det, phi, zside);
0154   if (depth < 0)
0155     depth = (det == 2) ? depthEta16[1] : depthEta16[0];
0156 #ifdef EDM_ML_DEBUG
0157   edm::LogVerbatim("HCalGeom") << "getDepthEta16: " << det << ":" << depth;
0158 #endif
0159   return depth;
0160 }
0161 
0162 int HcalDDDSimConstants::getDepthEta16M(const int& det) const {
0163   int depth = (det == 2) ? depthEta16[1] : depthEta16[0];
0164   std::vector<int> phis;
0165   int detsp = ldmap_.validDet(phis);
0166   if (detsp == det) {
0167     int zside = (phis[0] > 0) ? 1 : -1;
0168     int iphi = (phis[0] > 0) ? phis[0] : -phis[0];
0169     int depthsp = ldmap_.getDepth16(det, iphi, zside);
0170     if (det == 1 && depthsp > depth)
0171       depth = depthsp;
0172     if (det == 2 && depthsp < depth)
0173       depth = depthsp;
0174   }
0175   return depth;
0176 }
0177 
0178 int HcalDDDSimConstants::getDepthEta29(const int& phi, const int& zside, const int& i) const {
0179   int depth = (i == 0) ? ldmap_.getMaxDepthLastHE(2, phi, zside) : -1;
0180   if (depth < 0)
0181     depth = (i == 1) ? depthEta29[1] : depthEta29[0];
0182   return depth;
0183 }
0184 
0185 int HcalDDDSimConstants::getDepthEta29M(const int& i, const bool& planOne) const {
0186   int depth = (i == 1) ? depthEta29[1] : depthEta29[0];
0187   if (i == 0 && planOne) {
0188     std::vector<int> phis;
0189     int detsp = ldmap_.validDet(phis);
0190     if (detsp == 2) {
0191       int zside = (phis[0] > 0) ? 1 : -1;
0192       int iphi = (phis[0] > 0) ? phis[0] : -phis[0];
0193       int depthsp = ldmap_.getMaxDepthLastHE(2, iphi, zside);
0194       if (depthsp > depth)
0195         depth = depthsp;
0196     }
0197   }
0198   return depth;
0199 }
0200 
0201 std::pair<int, double> HcalDDDSimConstants::getDetEta(const double& eta, const int& depth) const {
0202   int hsubdet(0), ieta(0);
0203   double etaR(0);
0204   double heta = fabs(eta);
0205   for (int i = 0; i < nEta; i++)
0206     if (heta > hpar->etaTable[i])
0207       ieta = i + 1;
0208   if (heta <= hpar->etaRange[1]) {
0209     if (((ieta == hpar->etaMin[1] && depth == depthEta16[1]) || (ieta > hpar->etaMax[0])) &&
0210         (ieta <= hpar->etaMax[1])) {
0211       hsubdet = static_cast<int>(HcalEndcap);
0212     } else {
0213       hsubdet = static_cast<int>(HcalBarrel);
0214     }
0215     etaR = eta;
0216   } else {
0217     hsubdet = static_cast<int>(HcalForward);
0218     double theta = 2. * atan(exp(-heta));
0219     double hR = zVcal * tan(theta);
0220     etaR = (eta >= 0. ? hR : -hR);
0221   }
0222   return std::pair<int, double>(hsubdet, etaR);
0223 }
0224 
0225 int HcalDDDSimConstants::getEta(const int& det, const int& lay, const double& hetaR) const {
0226   int ieta(0);
0227   if (det == static_cast<int>(HcalForward)) {  // Forward HCal
0228     ieta = hpar->etaMax[2];
0229     for (int i = nR - 1; i > 0; i--)
0230       if (hetaR < hpar->rTable[i])
0231         ieta = hpar->etaMin[2] + nR - i - 1;
0232   } else {  // Barrel or Endcap
0233     ieta = 1;
0234     for (int i = 0; i < nEta - 1; i++)
0235       if (hetaR > hpar->etaTable[i])
0236         ieta = i + 1;
0237     if (det == static_cast<int>(HcalBarrel)) {
0238       if (ieta > hpar->etaMax[0])
0239         ieta = hpar->etaMax[0];
0240       if (lay == maxLayer_) {
0241         if (hetaR > etaHO[1] && ieta == hpar->noff[2])
0242           ieta++;
0243       }
0244     } else if (det == static_cast<int>(HcalEndcap)) {
0245       if (ieta <= hpar->etaMin[1])
0246         ieta = hpar->etaMin[1];
0247     }
0248   }
0249   return ieta;
0250 }
0251 
0252 std::pair<int, int> HcalDDDSimConstants::getEtaDepth(
0253     const int& det, int etaR, const int& phi, const int& zside, int depth, const int& lay) const {
0254 #ifdef EDM_ML_DEBUG
0255   edm::LogVerbatim("HCalGeom") << "HcalDDDEsimConstants:getEtaDepth: I/P " << det << ":" << etaR << ":" << phi << ":"
0256                                << zside << ":" << depth << ":" << lay;
0257 #endif
0258   //Modify the depth index
0259   if ((det == static_cast<int>(HcalEndcap)) && (etaR == 17) && (lay == 1))
0260     etaR = 18;
0261   if (det == static_cast<int>(HcalForward)) {  // Forward HCal
0262   } else if (det == static_cast<int>(HcalOuter)) {
0263     depth = 4;
0264   } else {
0265     if (lay >= 0) {
0266       depth = layerGroup(det, etaR, phi, zside, lay - 1);
0267       if (etaR == hpar->noff[0] && lay > 1) {
0268         int kphi = phi + int((hpar->phioff[3] + 0.1) / hpar->phibin[etaR - 1]);
0269         kphi = (kphi - 1) % 4 + 1;
0270         if (kphi == 2 || kphi == 3)
0271           depth = layerGroup(det, etaR, phi, zside, lay - 2);
0272       }
0273     } else if (det == static_cast<int>(HcalBarrel)) {
0274       if (depth > getMaxDepth(det, etaR, phi, zside, false))
0275         depth = getMaxDepth(det, etaR, phi, zside, false);
0276     }
0277     if (etaR >= hpar->noff[1] && depth > getDepthEta29(phi, zside, 0)) {
0278       etaR -= getDepthEta29(phi, zside, 1);
0279     } else if (etaR == hpar->etaMin[1]) {
0280       if (det == static_cast<int>(HcalBarrel)) {
0281         if (depth > getDepthEta16(det, phi, zside))
0282           depth = getDepthEta16(det, phi, zside);
0283       } else {
0284         if (depth < getDepthEta16(det, phi, zside))
0285           depth = getDepthEta16(det, phi, zside);
0286       }
0287     }
0288   }
0289 #ifdef EDM_ML_DEBUG
0290   edm::LogVerbatim("HCalGeom") << "HcalDDDEsimConstants:getEtaDepth: O/P " << etaR << ":" << depth;
0291 #endif
0292   return std::pair<int, int>(etaR, depth);
0293 }
0294 
0295 double HcalDDDSimConstants::getEtaHO(const double& etaR, const double& x, const double& y, const double& z) const {
0296   if (hpar->zHO.size() > 4) {
0297     double eta = fabs(etaR);
0298     double r = std::sqrt(x * x + y * y);
0299     if (r > rminHO) {
0300       double zz = fabs(z);
0301       if (zz > hpar->zHO[3]) {
0302         if (eta <= hpar->etaTable[10])
0303           eta = hpar->etaTable[10] + 0.001;
0304       } else if (zz > hpar->zHO[1]) {
0305         if (eta <= hpar->etaTable[4])
0306           eta = hpar->etaTable[4] + 0.001;
0307       }
0308     }
0309     eta = (z >= 0. ? eta : -eta);
0310 #ifdef EDM_ML_DEBUG
0311     edm::LogVerbatim("HCalGeom") << "R " << r << " Z " << z << " eta " << etaR << ":" << eta;
0312     if (eta != etaR)
0313       edm::LogInfo("HCalGeom") << "**** Check *****";
0314 #endif
0315     return eta;
0316   } else {
0317     return etaR;
0318   }
0319 }
0320 
0321 int HcalDDDSimConstants::getFrontLayer(const int& det, const int& eta) const {
0322   int lay = 0;
0323   if (det == 1) {
0324     if (std::abs(eta) == 16)
0325       lay = layFHB[1];
0326     else
0327       lay = layFHB[0];
0328   } else {
0329     if (std::abs(eta) == 16)
0330       lay = layFHE[1];
0331     else if (std::abs(eta) == 18)
0332       lay = layFHE[2];
0333     else
0334       lay = layFHE[0];
0335   }
0336   return lay;
0337 }
0338 
0339 int HcalDDDSimConstants::getLastLayer(const int& det, const int& eta) const {
0340   int lay = 0;
0341   if (det == 1) {
0342     if (std::abs(eta) == 15)
0343       lay = layBHB[1];
0344     else if (std::abs(eta) == 16)
0345       lay = layBHB[2];
0346     else
0347       lay = layBHB[0];
0348   } else {
0349     if (std::abs(eta) == 16)
0350       lay = layBHE[1];
0351     else if (std::abs(eta) == 17)
0352       lay = layBHE[2];
0353     else if (std::abs(eta) == 18)
0354       lay = layBHE[3];
0355     else
0356       lay = layBHE[0];
0357   }
0358   return lay;
0359 }
0360 
0361 double HcalDDDSimConstants::getLayer0Wt(const int& det, const int& phi, const int& zside) const {
0362   double wt = ldmap_.getLayer0Wt(det, phi, zside);
0363   if (wt < 0)
0364     wt = (det == 2) ? hpar->Layer0Wt[1] : hpar->Layer0Wt[0];
0365   return wt;
0366 }
0367 
0368 int HcalDDDSimConstants::getLayerFront(
0369     const int& det, const int& eta, const int& phi, const int& zside, const int& depth) const {
0370   int layer = ldmap_.getLayerFront(det, eta, phi, zside, depth);
0371   if (layer < 0) {
0372     if (det == 1 || det == 2) {
0373       layer = 1;
0374       for (int l = 0; l < getLayerMax(eta, depth); ++l) {
0375         if ((int)(layerGroup(eta - 1, l)) == depth) {
0376           layer = l + 1;
0377           break;
0378         }
0379       }
0380     } else {
0381       layer = (eta > hpar->noff[2]) ? maxLayerHB_ + 1 : maxLayer_;
0382     }
0383   }
0384   return layer;
0385 }
0386 
0387 int HcalDDDSimConstants::getLayerBack(
0388     const int& det, const int& eta, const int& phi, const int& zside, const int& depth) const {
0389   int layer = ldmap_.getLayerBack(det, eta, phi, zside, depth);
0390   if (layer < 0) {
0391     if (det == 1 || det == 2) {
0392       layer = depths[depth - 1][eta - 1];
0393     } else {
0394       layer = maxLayer_;
0395     }
0396   }
0397   return layer;
0398 }
0399 
0400 int HcalDDDSimConstants::getLayerMax(const int& eta, const int& depth) const {
0401   int layermx = ((eta < hpar->etaMin[1]) && depth - 1 < maxDepth[0]) ? maxLayerHB_ + 1 : (int)layerGroupSize(eta - 1);
0402   return layermx;
0403 }
0404 
0405 int HcalDDDSimConstants::getMaxDepth(
0406     const int& det, const int& eta, const int& phi, const int& zside, const bool& partialDetOnly) const {
0407   int dmax(-1);
0408   if (partialDetOnly) {
0409     if (ldmap_.isValid(det, phi, zside)) {
0410       dmax = ldmap_.getDepths(eta).second;
0411     }
0412   } else if (det == 1 || det == 2) {
0413     if (ldmap_.isValid(det, phi, zside))
0414       dmax = ldmap_.getDepths(eta).second;
0415     else if (det == 2)
0416       dmax = (maxDepth[1] > 0) ? layerGroup(eta - 1, maxLayer_) : 0;
0417     else if (eta == hpar->etaMax[0])
0418       dmax = getDepthEta16(det, phi, zside);
0419     else
0420       dmax = layerGroup(eta - 1, maxLayerHB_);
0421   } else if (det == 3) {  // HF
0422     dmax = maxHFDepth(zside * eta, phi);
0423   } else if (det == 4) {  // HO
0424     dmax = maxDepth[3];
0425   } else {
0426     dmax = -1;
0427   }
0428   return dmax;
0429 }
0430 
0431 int HcalDDDSimConstants::getMinDepth(
0432     const int& det, const int& eta, const int& phi, const int& zside, const bool& partialDetOnly) const {
0433   int lmin(-1);
0434   if (partialDetOnly) {
0435     if (ldmap_.isValid(det, phi, zside)) {
0436       lmin = ldmap_.getDepths(eta).first;
0437     }
0438   } else if (det == 3) {  // HF
0439     lmin = 1;
0440   } else if (det == 4) {  // HO
0441     lmin = maxDepth[3];
0442   } else {
0443     if (ldmap_.isValid(det, phi, zside)) {
0444       lmin = ldmap_.getDepths(eta).first;
0445     } else if (layerGroupSize(eta - 1) > 0) {
0446       lmin = (int)(layerGroup(eta - 1, 0));
0447       unsigned int type = (det == 1) ? 0 : 1;
0448       if (type == 1 && eta == hpar->etaMin[1])
0449         lmin = getDepthEta16(det, phi, zside);
0450     } else {
0451       lmin = 1;
0452     }
0453   }
0454   return lmin;
0455 }
0456 
0457 std::pair<int, int> HcalDDDSimConstants::getModHalfHBHE(const int& type) const {
0458   if (type == 0) {
0459     return std::pair<int, int>(nmodHB, nzHB);
0460   } else {
0461     return std::pair<int, int>(nmodHE, nzHE);
0462   }
0463 }
0464 
0465 std::pair<double, double> HcalDDDSimConstants::getPhiCons(const int& det, const int& ieta) const {
0466   double fioff(0), fibin(0);
0467   if (det == static_cast<int>(HcalForward)) {  // Forward HCal
0468     fioff = hpar->phioff[2];
0469     fibin = hpar->phitable[ieta - hpar->etaMin[2]];
0470     if (unitPhi(fibin) > 2) {  // HF double-phi
0471       fioff = hpar->phioff[4];
0472     }
0473   } else {  // Barrel or Endcap
0474     if (det == static_cast<int>(HcalEndcap)) {
0475       fioff = hpar->phioff[1];
0476     } else {
0477       fioff = hpar->phioff[0];
0478     }
0479     fibin = hpar->phibin[ieta - 1];
0480   }
0481   return std::pair<double, double>(fioff, fibin);
0482 }
0483 
0484 std::vector<std::pair<int, double> > HcalDDDSimConstants::getPhis(const int& subdet, const int& ieta) const {
0485   std::vector<std::pair<int, double> > phis;
0486   int ietaAbs = (ieta > 0) ? ieta : -ieta;
0487   std::pair<double, double> ficons = getPhiCons(subdet, ietaAbs);
0488   int nphi = int((2._pi + 0.1 * ficons.second) / ficons.second);
0489   int units = unitPhi(subdet, ietaAbs);
0490   for (int ifi = 0; ifi < nphi; ++ifi) {
0491     double phi = -ficons.first + (ifi + 0.5) * ficons.second;
0492     int iphi = phiNumber(ifi + 1, units);
0493     phis.emplace_back(std::pair<int, double>(iphi, phi));
0494   }
0495 #ifdef EDM_ML_DEBUG
0496   edm::LogVerbatim("HCalGeom") << "getPhis: subdet|ieta|iphi " << subdet << "|" << ieta << " with " << phis.size()
0497                                << " phi bins";
0498   for (unsigned int k = 0; k < phis.size(); ++k)
0499     edm::LogVerbatim("HCalGeom") << "[" << k << "] iphi " << phis[k].first << " phi "
0500                                  << convertRadToDeg(phis[k].second);
0501 #endif
0502   return phis;
0503 }
0504 
0505 std::vector<HcalCellType> HcalDDDSimConstants::HcalCellTypes() const {
0506   std::vector<HcalCellType> cellTypes = HcalCellTypes(HcalBarrel);
0507 #ifdef EDM_ML_DEBUG
0508   edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: " << cellTypes.size() << " cells of type HCal Barrel\n";
0509   for (unsigned int i = 0; i < cellTypes.size(); i++)
0510     edm::LogInfo("HCalGeom") << "Cell " << i << " " << cellTypes[i] << "\n";
0511 #endif
0512 
0513   std::vector<HcalCellType> hoCells = HcalCellTypes(HcalOuter);
0514 #ifdef EDM_ML_DEBUG
0515   edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: " << hoCells.size() << " cells of type HCal Outer\n";
0516   for (unsigned int i = 0; i < hoCells.size(); i++)
0517     edm::LogInfo("HCalGeom") << "Cell " << i << " " << hoCells[i] << "\n";
0518 #endif
0519   cellTypes.insert(cellTypes.end(), hoCells.begin(), hoCells.end());
0520 
0521   std::vector<HcalCellType> heCells = HcalCellTypes(HcalEndcap);
0522 #ifdef EDM_ML_DEBUG
0523   edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: " << heCells.size() << " cells of type HCal Endcap\n";
0524   for (unsigned int i = 0; i < heCells.size(); i++)
0525     edm::LogInfo("HCalGeom") << "Cell " << i << " " << heCells[i] << "\n";
0526 #endif
0527   cellTypes.insert(cellTypes.end(), heCells.begin(), heCells.end());
0528 
0529   std::vector<HcalCellType> hfCells = HcalCellTypes(HcalForward);
0530 #ifdef EDM_ML_DEBUG
0531   edm::LogInfo("HCalGeom") << "HcalDDDSimConstants: " << hfCells.size() << " cells of type HCal Forward\n";
0532   for (unsigned int i = 0; i < hfCells.size(); i++)
0533     edm::LogInfo("HCalGeom") << "Cell " << i << " " << hfCells[i] << "\n";
0534 #endif
0535   cellTypes.insert(cellTypes.end(), hfCells.begin(), hfCells.end());
0536 
0537   return cellTypes;
0538 }
0539 
0540 std::vector<HcalCellType> HcalDDDSimConstants::HcalCellTypes(const HcalSubdetector& subdet,
0541                                                              int ieta,
0542                                                              int depthl) const {
0543   std::vector<HcalCellType> cellTypes;
0544   if (subdet == HcalForward) {
0545     if (dzVcal < 0)
0546       return cellTypes;
0547   }
0548 
0549   int dmin, dmax, indx, nz;
0550   double hsize = 0;
0551   switch (subdet) {
0552     case HcalEndcap:
0553       dmin = 1;
0554       dmax = (maxDepth[1] > 0) ? maxLayer_ + 1 : 0;
0555       indx = 1;
0556       nz = nzHE;
0557       break;
0558     case HcalForward:
0559       dmin = 1;
0560       dmax = (!idHF2QIE.empty()) ? 2 : maxDepth[2];
0561       indx = 2;
0562       nz = 2;
0563       break;
0564     case HcalOuter:
0565       dmin = 4;
0566       dmax = 4;
0567       indx = 0;
0568       nz = nzHB;
0569       break;
0570     default:
0571       dmin = 1;
0572       dmax = maxLayerHB_ + 1;
0573       indx = 0;
0574       nz = nzHB;
0575       break;
0576   }
0577   if (depthl > 0)
0578     dmin = dmax = depthl;
0579   int ietamin = (ieta > 0) ? ieta : hpar->etaMin[indx];
0580   int ietamax = (ieta > 0) ? ieta : hpar->etaMax[indx];
0581   int phi = (indx == 2) ? 3 : 1;
0582 
0583   // Get the Cells
0584   int subdet0 = static_cast<int>(subdet);
0585   for (int depth = dmin; depth <= dmax; depth++) {
0586     int shift = getShift(subdet, depth);
0587     double gain = getGain(subdet, depth);
0588     if (subdet == HcalForward) {
0589       if (depth % 2 == 1)
0590         hsize = dzVcal;
0591       else
0592         hsize = dzVcal - 0.5 * dlShort;
0593     }
0594     for (int eta = ietamin; eta <= ietamax; eta++) {
0595       for (int iz = 0; iz < nz; ++iz) {
0596         int zside = -2 * iz + 1;
0597         HcalCellType::HcalCell temp1 = cell(subdet0, zside, depth, eta, phi);
0598         if (temp1.ok) {
0599           std::vector<std::pair<int, double> > phis = getPhis(subdet0, eta);
0600           HcalCellType temp2(subdet, eta, zside, depth, temp1, shift, gain, hsize);
0601           double dphi, fioff;
0602           std::vector<int> phiMiss2;
0603           if ((subdet == HcalBarrel) || (subdet == HcalOuter)) {
0604             fioff = hpar->phioff[0];
0605             dphi = hpar->phibin[eta - 1];
0606             if (subdet == HcalOuter) {
0607               if (eta == hpar->noff[4]) {
0608                 int kk = (iz == 0) ? 7 : (7 + hpar->noff[5]);
0609                 for (int miss = 0; miss < hpar->noff[5 + iz]; miss++) {
0610                   phiMiss2.emplace_back(hpar->noff[kk]);
0611                   kk++;
0612                 }
0613               }
0614             }
0615           } else if (subdet == HcalEndcap) {
0616             fioff = hpar->phioff[1];
0617             dphi = hpar->phibin[eta - 1];
0618           } else {
0619             fioff = hpar->phioff[2];
0620             dphi = hpar->phitable[eta - hpar->etaMin[2]];
0621             if (unitPhi(dphi) > 2)
0622               fioff = hpar->phioff[4];
0623           }
0624           int unit = unitPhi(dphi);
0625           temp2.setPhi(phis, phiMiss2, fioff, dphi, unit);
0626           cellTypes.emplace_back(temp2);
0627           // For HF look at extra cells
0628           if ((subdet == HcalForward) && (!idHF2QIE.empty())) {
0629             HcalCellType temp3(subdet, eta, zside + 2, depth, temp1, shift, gain, hsize);
0630             std::vector<int> phiMiss3;
0631             for (auto& phi : phis) {
0632               bool ok(false);
0633               for (auto l : idHF2QIE) {
0634                 if ((eta * zside == l.ieta()) && (phi.first == l.iphi())) {
0635                   ok = true;
0636                   break;
0637                 }
0638               }
0639               if (!ok)
0640                 phiMiss3.emplace_back(phi.first);
0641             }
0642             dphi = hpar->phitable[eta - hpar->etaMin[2]];
0643             unit = unitPhi(dphi);
0644             fioff = (unit > 2) ? hpar->phioff[4] : hpar->phioff[2];
0645             temp3.setPhi(phis, phiMiss2, fioff, dphi, unit);
0646             cellTypes.emplace_back(temp3);
0647           }
0648         }
0649       }
0650     }
0651   }
0652   return cellTypes;
0653 }
0654 
0655 int HcalDDDSimConstants::maxHFDepth(const int& eta, const int& iphi) const {
0656   int mxdepth = maxDepth[2];
0657   if (!idHF2QIE.empty()) {
0658     bool ok(false);
0659     for (auto k : idHF2QIE) {
0660       if ((eta == k.ieta()) && (iphi == k.iphi())) {
0661         ok = true;
0662         break;
0663       }
0664     }
0665     if (!ok)
0666       mxdepth = 2;
0667   }
0668   return mxdepth;
0669 }
0670 
0671 unsigned int HcalDDDSimConstants::numberOfCells(const HcalSubdetector& subdet) const {
0672   unsigned int num = 0;
0673   std::vector<HcalCellType> cellTypes = HcalCellTypes(subdet);
0674   for (auto& cellType : cellTypes) {
0675     num += (unsigned int)(cellType.nPhiBins());
0676   }
0677 #ifdef EDM_ML_DEBUG
0678   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants:numberOfCells " << cellTypes.size() << " " << num
0679                                << " for subdetector " << subdet;
0680 #endif
0681   return num;
0682 }
0683 
0684 int HcalDDDSimConstants::phiNumber(const int& phi, const int& units) const {
0685   int iphi_skip = phi;
0686   if (units == 2)
0687     iphi_skip = (phi - 1) * 2 + 1;
0688   else if (units == 4)
0689     iphi_skip = (phi - 1) * 4 - 1;
0690   if (iphi_skip < 0)
0691     iphi_skip += 72;
0692   return iphi_skip;
0693 }
0694 
0695 void HcalDDDSimConstants::printTiles() const {
0696   std::vector<int> phis;
0697   int detsp = ldmap_.validDet(phis);
0698   int kphi = (detsp > 0) ? phis[0] : 1;
0699   int zside = (kphi > 0) ? 1 : -1;
0700   int iphi = (kphi > 0) ? kphi : -kphi;
0701   edm::LogVerbatim("HCalGeom") << "Tile Information for HB from " << hpar->etaMin[0] << " to " << hpar->etaMax[0]
0702                                << "\n";
0703   for (int eta = hpar->etaMin[0]; eta <= hpar->etaMax[0]; eta++) {
0704     int dmax = getMaxDepth(1, eta, iphi, -zside, false);
0705     for (int depth = 1; depth <= dmax; depth++)
0706       printTileHB(eta, iphi, -zside, depth);
0707     if (detsp == 1) {
0708       int dmax = getMaxDepth(1, eta, iphi, zside, false);
0709       for (int depth = 1; depth <= dmax; depth++)
0710         printTileHB(eta, iphi, zside, depth);
0711     }
0712   }
0713 
0714   edm::LogVerbatim("HCalGeom") << "\nTile Information for HE from " << hpar->etaMin[1] << " to " << hpar->etaMax[1]
0715                                << "\n";
0716   for (int eta = hpar->etaMin[1]; eta <= hpar->etaMax[1]; eta++) {
0717     int dmin = (eta == hpar->etaMin[1]) ? getDepthEta16(2, iphi, -zside) : 1;
0718     int dmax = getMaxDepth(2, eta, iphi, -zside, false);
0719     for (int depth = dmin; depth <= dmax; depth++)
0720       printTileHE(eta, iphi, -zside, depth);
0721     if (detsp == 2) {
0722       int dmax = getMaxDepth(2, eta, iphi, zside, false);
0723       for (int depth = 1; depth <= dmax; depth++)
0724         printTileHE(eta, iphi, zside, depth);
0725     }
0726   }
0727 }
0728 
0729 int HcalDDDSimConstants::unitPhi(const int& det, const int& etaR) const {
0730   double dphi =
0731       (det == static_cast<int>(HcalForward)) ? hpar->phitable[etaR - hpar->etaMin[2]] : hpar->phibin[etaR - 1];
0732   return unitPhi(dphi);
0733 }
0734 
0735 int HcalDDDSimConstants::unitPhi(const double& dphi) const {
0736   const double fiveDegInRad = 2 * M_PI / 72;
0737   int units = int(dphi / fiveDegInRad + 0.5);
0738   if (units < 1)
0739     units = 1;
0740   return units;
0741 }
0742 
0743 void HcalDDDSimConstants::initialize(void) {
0744   nEta = hpar->etaTable.size();
0745   nR = hpar->rTable.size();
0746   nPhiF = nR - 1;
0747   isBH_ = false;
0748 
0749 #ifdef EDM_ML_DEBUG
0750   for (int i = 0; i < nEta - 1; ++i) {
0751     edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants:Read LayerGroup" << i << ":";
0752     for (unsigned int k = 0; k < layerGroupSize(i); k++)
0753       edm::LogVerbatim("HCalGeom") << " [" << k << "] = " << layerGroup(i, k);
0754   }
0755 #endif
0756 
0757   // Geometry parameters for HF
0758   dlShort = hpar->gparHF[0];
0759   zVcal = hpar->gparHF[4];
0760   dzVcal = hpar->dzVcal;
0761 #ifdef EDM_ML_DEBUG
0762   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: dlShort " << dlShort << " zVcal " << zVcal << " and dzVcal "
0763                                << dzVcal;
0764 #endif
0765 
0766   //Transform some of the parameters
0767   maxDepth = hpar->maxDepth;
0768   maxDepth[0] = maxDepth[1] = 0;
0769   for (int i = 0; i < nEta - 1; ++i) {
0770     unsigned int imx = layerGroupSize(i);
0771     int laymax = (imx > 0) ? layerGroup(i, imx - 1) : 0;
0772     if (i < hpar->etaMax[0]) {
0773       int laymax0 = (imx > 16) ? layerGroup(i, 16) : laymax;
0774       if (i + 1 == hpar->etaMax[0] && laymax0 > 2)
0775         laymax0 = 2;
0776       if (maxDepth[0] < laymax0)
0777         maxDepth[0] = laymax0;
0778     }
0779     if (i >= hpar->etaMin[1] - 1 && i < hpar->etaMax[1]) {
0780       if (maxDepth[1] < laymax)
0781         maxDepth[1] = laymax;
0782     }
0783   }
0784 #ifdef EDM_ML_DEBUG
0785   for (int i = 0; i < 4; ++i)
0786     edm::LogVerbatim("HCalGeom") << "Detector Type [" << i << "] iEta " << hpar->etaMin[i] << ":" << hpar->etaMax[i]
0787                                  << " MaxDepth " << maxDepth[i];
0788 #endif
0789 
0790   int maxdepth = (maxDepth[1] > maxDepth[0]) ? maxDepth[1] : maxDepth[0];
0791   for (int i = 0; i < maxdepth; ++i) {
0792     for (int k = 0; k < nEta - 1; ++k) {
0793       int layermx = getLayerMax(k + 1, i + 1);
0794       int ll = layermx;
0795       for (int l = layermx - 1; l >= 0; --l) {
0796         if ((int)layerGroup(k, l) == i + 1) {
0797           ll = l + 1;
0798           break;
0799         }
0800       }
0801       depths[i].emplace_back(ll);
0802     }
0803 
0804 #ifdef EDM_ML_DEBUG
0805     edm::LogVerbatim("HCalGeom") << "Depth " << i + 1 << " with " << depths[i].size() << " etas:";
0806     for (int k = 0; k < nEta - 1; ++k)
0807       edm::LogVerbatim("HCalGeom") << " [" << k << "] " << depths[i][k];
0808 #endif
0809   }
0810 
0811   nzHB = hpar->modHB[1];
0812   nmodHB = hpar->modHB[0];
0813   nzHE = hpar->modHE[1];
0814   nmodHE = hpar->modHE[0];
0815 #ifdef EDM_ML_DEBUG
0816   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants:: " << nzHB << ":" << nmodHB << " barrel and " << nzHE << ":"
0817                                << nmodHE << " endcap half-sectors";
0818 #endif
0819 
0820   if (hpar->rHB.size() > maxLayerHB_ + 1 && hpar->zHO.size() > 4) {
0821     rminHO = hpar->rHO[0];
0822     for (int k = 0; k < 4; ++k)
0823       etaHO[k] = hpar->rHO[k + 1];
0824   } else {
0825     rminHO = -1.0;
0826     etaHO[0] = hpar->etaTable[4];
0827     etaHO[1] = hpar->etaTable[4];
0828     etaHO[2] = hpar->etaTable[10];
0829     etaHO[3] = hpar->etaTable[10];
0830   }
0831 #ifdef EDM_ML_DEBUG
0832   edm::LogVerbatim("HCalGeom") << "HO Eta boundaries " << etaHO[0] << " " << etaHO[1] << " " << etaHO[2] << " "
0833                                << etaHO[3];
0834   edm::LogVerbatim("HCalGeom") << "HO Parameters " << rminHO << " " << hpar->zHO.size();
0835   for (int i = 0; i < 4; ++i)
0836     edm::LogVerbatim("HCalGeom") << " eta[" << i << "] = " << etaHO[i];
0837   for (unsigned int i = 0; i < hpar->zHO.size(); ++i)
0838     edm::LogVerbatim("HCalGeom") << " zHO[" << i << "] = " << hpar->zHO[i];
0839 #endif
0840 
0841   int noffsize = 7 + hpar->noff[5] + hpar->noff[6];
0842   int noffl(noffsize + 5);
0843   if ((int)(hpar->noff.size()) > (noffsize + 3)) {
0844     depthEta16[0] = hpar->noff[noffsize];
0845     depthEta16[1] = hpar->noff[noffsize + 1];
0846     depthEta29[0] = hpar->noff[noffsize + 2];
0847     depthEta29[1] = hpar->noff[noffsize + 3];
0848     if ((int)(hpar->noff.size()) > (noffsize + 4)) {
0849       noffl += (2 * hpar->noff[noffsize + 4]);
0850       if ((int)(hpar->noff.size()) > noffl)
0851         isBH_ = (hpar->noff[noffl] > 0);
0852     }
0853   } else {
0854     depthEta16[0] = 2;
0855     depthEta16[1] = 3;
0856     depthEta29[0] = 2;
0857     depthEta29[1] = 1;
0858   }
0859 #ifdef EDM_ML_DEBUG
0860   edm::LogVerbatim("HCalGeom") << "isBH_ " << hpar->noff.size() << ":" << noffsize << ":" << noffl << ":" << isBH_;
0861   edm::LogVerbatim("HCalGeom") << "Depth index at ieta = 16 for HB (max) " << depthEta16[0] << " HE (min) "
0862                                << depthEta16[1] << "; max depth for itea = 29 : (" << depthEta29[0] << ":"
0863                                << depthEta29[1] << ")";
0864 #endif
0865 
0866   if ((int)(hpar->noff.size()) > (noffsize + 4)) {
0867     int npair = hpar->noff[noffsize + 4];
0868     int kk = noffsize + 4;
0869     for (int k = 0; k < npair; ++k) {
0870       idHF2QIE.emplace_back(HcalDetId(HcalForward, hpar->noff[kk + 1], hpar->noff[kk + 2], 1));
0871       kk += 2;
0872     }
0873   }
0874 #ifdef EDM_ML_DEBUG
0875   edm::LogVerbatim("HCalGeom") << idHF2QIE.size() << " detector channels having 2 QIE cards:";
0876   for (unsigned int k = 0; k < idHF2QIE.size(); ++k)
0877     edm::LogVerbatim("HCalGeom") << " [" << k << "] " << idHF2QIE[k];
0878 #endif
0879 
0880   layFHB[0] = 0;
0881   layFHB[1] = 1;
0882   layBHB[0] = 16;
0883   layBHB[1] = 15;
0884   layBHB[2] = 8;
0885   if (maxDepth[1] == 0) {
0886     layFHE[0] = layFHE[1] = layFHE[2] = 0;
0887     layBHE[0] = layBHE[1] = layBHE[2] = layBHE[3] = 0;
0888   } else {
0889     layFHE[0] = 1;
0890     layFHE[1] = 4;
0891     layFHE[2] = 0;
0892     layBHE[0] = 18;
0893     layBHE[1] = 9;
0894     layBHE[2] = 14;
0895     layBHE[3] = 16;
0896   }
0897   depthMaxSp_ = std::pair<int, int>(0, 0);
0898   int noffk(noffsize + 5);
0899   if ((int)(hpar->noff.size()) > (noffsize + 5)) {
0900     noffk += (2 * hpar->noff[noffsize + 4]);
0901     if ((int)(hpar->noff.size()) >= noffk + 7) {
0902       int dtype = hpar->noff[noffk + 1];
0903       int nphi = hpar->noff[noffk + 2];
0904       int ndeps = hpar->noff[noffk + 3];
0905       int ndp16 = hpar->noff[noffk + 4];
0906       int ndp29 = hpar->noff[noffk + 5];
0907       double wt = 0.1 * (hpar->noff[noffk + 6]);
0908       if ((int)(hpar->noff.size()) >= (noffk + 7 + nphi + 3 * ndeps)) {
0909         if (dtype == 1 || dtype == 2) {
0910           std::vector<int> ifi, iet, ily, idp;
0911           for (int i = 0; i < nphi; ++i)
0912             ifi.emplace_back(hpar->noff[noffk + 7 + i]);
0913           for (int i = 0; i < ndeps; ++i) {
0914             iet.emplace_back(hpar->noff[noffk + 7 + nphi + 3 * i]);
0915             ily.emplace_back(hpar->noff[noffk + 7 + nphi + 3 * i + 1]);
0916             idp.emplace_back(hpar->noff[noffk + 7 + nphi + 3 * i + 2]);
0917           }
0918 #ifdef EDM_ML_DEBUG
0919           edm::LogVerbatim("HCalGeom") << "Initialize HcalLayerDepthMap for "
0920                                        << "Detector " << dtype << " etaMax " << hpar->etaMax[dtype] << " with " << nphi
0921                                        << " sectors";
0922           for (int i = 0; i < nphi; ++i)
0923             edm::LogVerbatim("HCalGeom") << " [" << i << "] " << ifi[i];
0924           edm::LogVerbatim("HCalGeom") << "And " << ndeps << " depth sections";
0925           for (int i = 0; i < ndeps; ++i)
0926             edm::LogVerbatim("HCalGeom") << " [" << i << "] " << iet[i] << "  " << ily[i] << "  " << idp[i];
0927           edm::LogVerbatim("HCalGeom") << "Maximum depth for last HE Eta tower " << depthEta29[0] << ":" << ndp16 << ":"
0928                                        << ndp29 << " L0 Wt " << hpar->Layer0Wt[dtype - 1] << ":" << wt;
0929 #endif
0930           ldmap_.initialize(dtype, hpar->etaMax[dtype - 1], ndp16, ndp29, wt, ifi, iet, ily, idp);
0931           int zside = (ifi[0] > 0) ? 1 : -1;
0932           int iphi = (ifi[0] > 0) ? ifi[0] : -ifi[0];
0933           depthMaxSp_ = std::pair<int, int>(dtype, ldmap_.getDepthMax(dtype, iphi, zside));
0934         }
0935       }
0936       int noffm = (noffk + 7 + nphi + 3 * ndeps);
0937       if ((int)(hpar->noff.size()) > noffm) {
0938         int ndnext = hpar->noff[noffm];
0939         if (ndnext > 4 && (int)(hpar->noff.size()) >= noffm + ndnext) {
0940           for (int i = 0; i < 2; ++i)
0941             layFHB[i] = hpar->noff[noffm + i + 1];
0942           for (int i = 0; i < 3; ++i)
0943             layFHE[i] = hpar->noff[noffm + i + 3];
0944         }
0945         if (ndnext > 11 && (int)(hpar->noff.size()) >= noffm + ndnext) {
0946           for (int i = 0; i < 3; ++i)
0947             layBHB[i] = hpar->noff[noffm + i + 6];
0948           for (int i = 0; i < 4; ++i)
0949             layBHE[i] = hpar->noff[noffm + i + 9];
0950         }
0951       }
0952     }
0953   }
0954 #ifdef EDM_ML_DEBUG
0955   edm::LogVerbatim("HCalGeom") << "Front Layer Definition for HB: " << layFHB[0] << ":" << layFHB[1]
0956                                << " and for HE: " << layFHE[0] << ":" << layFHE[1] << ":" << layFHE[2];
0957   edm::LogVerbatim("HCalGeom") << "Last Layer Definition for HB: " << layBHB[0] << ":" << layBHB[1] << ":" << layBHB[2]
0958                                << " and for HE: " << layBHE[0] << ":" << layBHE[1] << ":" << layBHE[2] << ":"
0959                                << layBHE[3];
0960 #endif
0961   if (depthMaxSp_.first == 0) {
0962     depthMaxSp_ = depthMaxDf_ = std::pair<int, int>(2, maxDepth[1]);
0963   } else if (depthMaxSp_.first == 1) {
0964     depthMaxDf_ = std::pair<int, int>(1, maxDepth[0]);
0965     if (depthMaxSp_.second > maxDepth[0])
0966       maxDepth[0] = depthMaxSp_.second;
0967   } else {
0968     depthMaxDf_ = std::pair<int, int>(2, maxDepth[1]);
0969     if (depthMaxSp_.second > maxDepth[1])
0970       maxDepth[1] = depthMaxSp_.second;
0971   }
0972 #ifdef EDM_ML_DEBUG
0973   edm::LogVerbatim("HCalGeom") << "Detector type and maximum depth for all RBX " << depthMaxDf_.first << ":"
0974                                << depthMaxDf_.second << " and for special RBX " << depthMaxSp_.first << ":"
0975                                << depthMaxSp_.second;
0976 #endif
0977 }
0978 
0979 double HcalDDDSimConstants::deltaEta(const int& det, const int& etaR, const int& depth) const {
0980   double tmp = 0;
0981   if (det == static_cast<int>(HcalForward)) {
0982     int ir = nR + hpar->etaMin[2] - etaR - 1;
0983     if (ir > 0 && ir < nR) {
0984       double z = zVcal;
0985       if (depth % 2 != 1)
0986         z += dlShort;
0987       tmp = 0.5 * (getEta(hpar->rTable[ir - 1], z) - getEta(hpar->rTable[ir], z));
0988     }
0989   } else {
0990     if (etaR > 0 && etaR < nEta) {
0991       if (etaR == hpar->noff[1] - 1 && depth > depthEta29[0]) {
0992         tmp = 0.5 * (hpar->etaTable[etaR + 1] - hpar->etaTable[etaR - 1]);
0993       } else if (det == static_cast<int>(HcalOuter)) {
0994         if (etaR == hpar->noff[2]) {
0995           tmp = 0.5 * (etaHO[0] - hpar->etaTable[etaR - 1]);
0996         } else if (etaR == hpar->noff[2] + 1) {
0997           tmp = 0.5 * (hpar->etaTable[etaR] - etaHO[1]);
0998         } else if (etaR == hpar->noff[3]) {
0999           tmp = 0.5 * (etaHO[2] - hpar->etaTable[etaR - 1]);
1000         } else if (etaR == hpar->noff[3] + 1) {
1001           tmp = 0.5 * (hpar->etaTable[etaR] - etaHO[3]);
1002         } else {
1003           tmp = 0.5 * (hpar->etaTable[etaR] - hpar->etaTable[etaR - 1]);
1004         }
1005       } else {
1006         tmp = 0.5 * (hpar->etaTable[etaR] - hpar->etaTable[etaR - 1]);
1007       }
1008     }
1009   }
1010 #ifdef EDM_ML_DEBUG
1011   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants::deltaEta " << etaR << " " << depth << " ==> " << tmp;
1012 #endif
1013   return tmp;
1014 }
1015 
1016 double HcalDDDSimConstants::getEta(const int& det, const int& etaR, const int& zside, int depth) const {
1017   double tmp = 0;
1018   if (det == static_cast<int>(HcalForward)) {
1019     int ir = nR + hpar->etaMin[2] - etaR - 1;
1020     if (ir > 0 && ir < nR) {
1021       double z = zVcal;
1022       if (depth % 2 != 1)
1023         z += dlShort;
1024       tmp = 0.5 * (getEta(hpar->rTable[ir - 1], z) + getEta(hpar->rTable[ir], z));
1025     }
1026   } else {
1027     if (etaR > 0 && etaR < nEta) {
1028       if (etaR == hpar->noff[1] - 1 && depth > depthEta29[0]) {
1029         tmp = 0.5 * (hpar->etaTable[etaR + 1] + hpar->etaTable[etaR - 1]);
1030       } else if (det == static_cast<int>(HcalOuter)) {
1031         if (etaR == hpar->noff[2]) {
1032           tmp = 0.5 * (etaHO[0] + hpar->etaTable[etaR - 1]);
1033         } else if (etaR == hpar->noff[2] + 1) {
1034           tmp = 0.5 * (hpar->etaTable[etaR] + etaHO[1]);
1035         } else if (etaR == hpar->noff[3]) {
1036           tmp = 0.5 * (etaHO[2] + hpar->etaTable[etaR - 1]);
1037         } else if (etaR == hpar->noff[3] + 1) {
1038           tmp = 0.5 * (hpar->etaTable[etaR] + etaHO[3]);
1039         } else {
1040           tmp = 0.5 * (hpar->etaTable[etaR] + hpar->etaTable[etaR - 1]);
1041         }
1042       } else {
1043         tmp = 0.5 * (hpar->etaTable[etaR] + hpar->etaTable[etaR - 1]);
1044       }
1045     }
1046   }
1047   if (zside == 0)
1048     tmp = -tmp;
1049 #ifdef EDM_ML_DEBUG
1050   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants::getEta " << etaR << " " << zside << " " << depth << " ==> "
1051                                << tmp;
1052 #endif
1053   return tmp;
1054 }
1055 
1056 double HcalDDDSimConstants::getEta(const double& r, const double& z) const {
1057   double tmp = 0;
1058   if (z != 0)
1059     tmp = -log(tan(0.5 * atan(r / z)));
1060 #ifdef EDM_ML_DEBUG
1061   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants::getEta " << r << " " << z << " ==> " << tmp;
1062 #endif
1063   return tmp;
1064 }
1065 
1066 int HcalDDDSimConstants::getShift(const HcalSubdetector& subdet, const int& depth) const {
1067   int shift;
1068   switch (subdet) {
1069     case HcalEndcap:
1070       shift = hpar->HEShift[0];
1071       break;
1072     case HcalForward:
1073       shift = hpar->HFShift[(depth - 1) % 2];
1074       break;
1075     case HcalOuter:
1076       shift = hpar->HBShift[3];
1077       break;
1078     default:
1079       shift = hpar->HBShift[0];
1080       break;
1081   }
1082   return shift;
1083 }
1084 
1085 double HcalDDDSimConstants::getGain(const HcalSubdetector& subdet, const int& depth) const {
1086   double gain;
1087   switch (subdet) {
1088     case HcalEndcap:
1089       gain = hpar->HEGains[0];
1090       break;
1091     case HcalForward:
1092       gain = hpar->HFGains[(depth - 1) % 2];
1093       break;
1094     case HcalOuter:
1095       gain = hpar->HBGains[3];
1096       break;
1097     default:
1098       gain = hpar->HBGains[0];
1099       break;
1100   }
1101   return gain;
1102 }
1103 
1104 void HcalDDDSimConstants::printTileHB(const int& eta, const int& phi, const int& zside, const int& depth) const {
1105   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants::printTileHB for eta " << eta << " and depth " << depth;
1106 
1107   double etaL = hpar->etaTable.at(eta - 1);
1108   double thetaL = 2. * atan(exp(-etaL));
1109   double etaH = hpar->etaTable.at(eta);
1110   double thetaH = 2. * atan(exp(-etaH));
1111   int layL = getLayerFront(1, eta, phi, zside, depth);
1112   int layH = getLayerBack(1, eta, phi, zside, depth);
1113   edm::LogVerbatim("HCalGeom") << "\ntileHB:: eta|depth " << zside * eta << "|" << depth << " theta "
1114                                << convertRadToDeg(thetaH) << ":" << convertRadToDeg(thetaL) << " Layer " << layL - 1
1115                                << ":" << layH - 1;
1116   for (int lay = layL - 1; lay < layH; ++lay) {
1117     std::vector<double> area(2, 0);
1118     int kk(0);
1119     double mean(0);
1120     for (unsigned int k = 0; k < hpar->layHB.size(); ++k) {
1121       if (lay == hpar->layHB[k]) {
1122         double zmin = hpar->rhoxHB[k] * std::cos(thetaL) / std::sin(thetaL);
1123         double zmax = hpar->rhoxHB[k] * std::cos(thetaH) / std::sin(thetaH);
1124         double dz = (std::min(zmax, hpar->dxHB[k]) - zmin);
1125         if (dz > 0) {
1126           area[kk] = dz * hpar->dyHB[k];
1127           mean += area[kk];
1128           kk++;
1129         }
1130       }
1131     }
1132     if (area[0] > 0) {
1133       mean /= (kk * 100);
1134       edm::LogVerbatim("HCalGeom") << std::setw(2) << lay << " Area " << std::setw(8) << area[0] << " " << std::setw(8)
1135                                    << area[1] << " Mean " << mean;
1136     }
1137   }
1138 }
1139 
1140 void HcalDDDSimConstants::printTileHE(const int& eta, const int& phi, const int& zside, const int& depth) const {
1141   double etaL = hpar->etaTable[eta - 1];
1142   double thetaL = 2. * atan(exp(-etaL));
1143   double etaH = hpar->etaTable[eta];
1144   double thetaH = 2. * atan(exp(-etaH));
1145   int layL = getLayerFront(2, eta, phi, zside, depth);
1146   int layH = getLayerBack(2, eta, phi, zside, depth);
1147   double phib = hpar->phibin[eta - 1];
1148   int nphi = 2;
1149   if (phib > 6._deg)
1150     nphi = 1;
1151   edm::LogVerbatim("HCalGeom") << "\ntileHE:: Eta/depth " << zside * eta << "|" << depth << " theta "
1152                                << convertRadToDeg(thetaH) << ":" << convertRadToDeg(thetaL) << " Layer " << layL - 1
1153                                << ":" << layH - 1 << " phi " << nphi;
1154   for (int lay = layL - 1; lay < layH; ++lay) {
1155     std::vector<double> area(4, 0);
1156     int kk(0);
1157     double mean(0);
1158     for (unsigned int k = 0; k < hpar->layHE.size(); ++k) {
1159       if (lay == hpar->layHE[k]) {
1160         double rmin = hpar->zxHE[k] * std::tan(thetaH);
1161         double rmax = hpar->zxHE[k] * std::tan(thetaL);
1162         if ((lay != 0 || eta == 18) &&
1163             (lay != 1 || (eta == 18 && hpar->rhoxHE[k] - hpar->dyHE[k] > 1000) ||
1164              (eta != 18 && hpar->rhoxHE[k] - hpar->dyHE[k] < 1000)) &&
1165             rmin + 30 < hpar->rhoxHE[k] + hpar->dyHE[k] && rmax > hpar->rhoxHE[k] - hpar->dyHE[k]) {
1166           rmin = std::max(rmin, hpar->rhoxHE[k] - hpar->dyHE[k]);
1167           rmax = std::min(rmax, hpar->rhoxHE[k] + hpar->dyHE[k]);
1168           double dx1 = rmin * std::tan(phib);
1169           double dx2 = rmax * std::tan(phib);
1170           double ar1 = 0, ar2 = 0;
1171           if (nphi == 1) {
1172             ar1 = 0.5 * (rmax - rmin) * (dx1 + dx2 - 4. * hpar->dx1HE[k]);
1173             mean += ar1;
1174           } else {
1175             ar1 = 0.5 * (rmax - rmin) * (dx1 + dx2 - 2. * hpar->dx1HE[k]);
1176             ar2 = 0.5 * (rmax - rmin) * ((rmax + rmin) * tan(10._deg) - 4 * hpar->dx1HE[k]) - ar1;
1177             mean += (ar1 + ar2);
1178           }
1179           area[kk] = ar1;
1180           area[kk + 2] = ar2;
1181           kk++;
1182         }
1183       }
1184     }
1185     if (area[0] > 0 && area[1] > 0) {
1186       int lay0 = lay - 1;
1187       if (eta == 18)
1188         lay0++;
1189       if (nphi == 1) {
1190         mean /= (kk * 100);
1191         edm::LogVerbatim("HCalGeom") << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " "
1192                                      << std::setw(8) << area[1] << " Mean " << mean;
1193       } else {
1194         mean /= (kk * 200);
1195         edm::LogVerbatim("HCalGeom") << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " "
1196                                      << std::setw(8) << area[1] << ":" << std::setw(8) << area[2] << " " << std::setw(8)
1197                                      << area[3] << " Mean " << mean;
1198       }
1199     }
1200   }
1201 }
1202 
1203 unsigned int HcalDDDSimConstants::layerGroupSize(int eta) const {
1204   unsigned int k = 0;
1205   for (auto const& it : hpar->layerGroupEtaSim) {
1206     if (it.layer == (unsigned int)(eta + 1)) {
1207       return it.layerGroup.size();
1208     }
1209     if (it.layer > (unsigned int)(eta + 1))
1210       break;
1211     k = it.layerGroup.size();
1212   }
1213   return k;
1214 }
1215 
1216 unsigned int HcalDDDSimConstants::layerGroup(int eta, int i) const {
1217   unsigned int k = 0;
1218   for (auto const& it : hpar->layerGroupEtaSim) {
1219     if (it.layer == (unsigned int)(eta + 1)) {
1220       return it.layerGroup.at(i);
1221     }
1222     if (it.layer > (unsigned int)(eta + 1))
1223       break;
1224     k = it.layerGroup.at(i);
1225   }
1226   return k;
1227 }
1228 
1229 unsigned int HcalDDDSimConstants::layerGroup(int det, int eta, int phi, int zside, int lay) const {
1230   int depth0 = findDepth(det, eta, phi, zside, lay);
1231   unsigned int depth = (depth0 > 0) ? (unsigned int)(depth0) : layerGroup(eta - 1, lay);
1232   return depth;
1233 }