Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:49:32

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::LogVerbatim("HCalGeom") << "HcalDDDSimConstants::HcalDDDSimConstants (const HcalParameters* hp) constructor";
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)";
0019 #endif
0020 }
0021 
0022 HcalDDDSimConstants::~HcalDDDSimConstants() {
0023 #ifdef EDM_ML_DEBUG
0024   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants::destructed!!!";
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 > static_cast<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::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: wrong eta " << etaR << " (" << ir << "/" << nR
0080                                      << ") Detector " << idet;
0081 #endif
0082       }
0083     } else if (etaR <= nEta) {
0084       int laymin(depth), laymax(depth);
0085       if (idet == static_cast<int>(HcalOuter)) {
0086         laymin =
0087             (etaR > hpar->noff[2]) ? (static_cast<int>(hpar->zHE.size())) : (static_cast<int>(hpar->zHE.size())) - 1;
0088         laymax = (static_cast<int>(hpar->zHE.size()));
0089       }
0090       double d1 = 0, d2 = 0;
0091       if (idet == static_cast<int>(HcalEndcap)) {
0092         flagrz = false;
0093         d1 = hpar->zHE[laymin - 1] - hpar->dzHE[laymin - 1];
0094         d2 = hpar->zHE[laymax - 1] + hpar->dzHE[laymax - 1];
0095       } else {
0096         d1 = hpar->rHB[laymin - 1] - hpar->drHB[laymin - 1];
0097         d2 = hpar->rHB[laymax - 1] + hpar->drHB[laymax - 1];
0098       }
0099       rz = 0.5 * (d2 + d1);
0100       drz = 0.5 * (d2 - d1);
0101     } else {
0102       ok = false;
0103       edm::LogWarning("HCalGeom") << "HcalDDDSimConstants: wrong depth " << depth << " or etaR " << etaR
0104                                   << " for detector " << idet;
0105     }
0106   } else {
0107     ok = false;
0108     edm::LogWarning("HCalGeom") << "HcalDDDSimConstants: wrong depth " << depth << " det " << idet;
0109   }
0110   HcalCellType::HcalCell tmp(ok, eta, deta, phi, dphi, rz, drz, flagrz);
0111 
0112 #ifdef EDM_ML_DEBUG
0113   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: det/side/depth/etaR/"
0114                                << "phi " << idet << "/" << zside << "/" << depth << "/" << etaR << "/" << iphi
0115                                << " Cell Flag " << tmp.ok << " " << tmp.eta << " " << tmp.deta << " phi " << tmp.phi
0116                                << " " << tmp.dphi << " r(z) " << tmp.rz << " " << tmp.drz << " " << tmp.flagrz;
0117 #endif
0118   return tmp;
0119 }
0120 
0121 int HcalDDDSimConstants::findDepth(
0122     const int& det, const int& eta, const int& phi, const int& zside, const int& lay) const {
0123   int depth = (ldmap_.isValid(det, phi, zside)) ? ldmap_.getDepth(det, eta, phi, zside, lay) : -1;
0124   return depth;
0125 }
0126 
0127 unsigned int HcalDDDSimConstants::findLayer(const int& layer,
0128                                             const std::vector<HcalParameters::LayerItem>& layerGroup) const {
0129   unsigned int id = layerGroup.size();
0130   for (unsigned int i = 0; i < layerGroup.size(); i++) {
0131     if (layer == static_cast<int>(layerGroup[i].layer)) {
0132       id = i;
0133       break;
0134     }
0135   }
0136   return id;
0137 }
0138 
0139 std::vector<std::pair<double, double> > HcalDDDSimConstants::getConstHBHE(const int& type) const {
0140   std::vector<std::pair<double, double> > gcons;
0141   if (type == 0) {
0142     for (unsigned int i = 0; i < hpar->rHB.size(); ++i) {
0143       gcons.emplace_back(std::pair<double, double>(hpar->rHB[i], hpar->drHB[i]));
0144     }
0145   } else {
0146     for (unsigned int i = 0; i < hpar->zHE.size(); ++i) {
0147       gcons.emplace_back(std::pair<double, double>(hpar->zHE[i], hpar->dzHE[i]));
0148     }
0149   }
0150   return gcons;
0151 }
0152 
0153 int HcalDDDSimConstants::getDepthEta16(const int& det, const int& phi, const int& zside) const {
0154   int depth = ldmap_.getDepth16(det, phi, zside);
0155   if (depth < 0)
0156     depth = (det == 2) ? depthEta16[1] : depthEta16[0];
0157 #ifdef EDM_ML_DEBUG
0158   edm::LogVerbatim("HCalGeom") << "getDepthEta16: " << det << ":" << depth;
0159 #endif
0160   return depth;
0161 }
0162 
0163 int HcalDDDSimConstants::getDepthEta16M(const int& det) const {
0164   int depth = (det == 2) ? depthEta16[1] : depthEta16[0];
0165   std::vector<int> phis;
0166   int detsp = ldmap_.validDet(phis);
0167   if (detsp == det) {
0168     int zside = (phis[0] > 0) ? 1 : -1;
0169     int iphi = (phis[0] > 0) ? phis[0] : -phis[0];
0170     int depthsp = ldmap_.getDepth16(det, iphi, zside);
0171     if (det == 1 && depthsp > depth)
0172       depth = depthsp;
0173     if (det == 2 && depthsp < depth)
0174       depth = depthsp;
0175   }
0176   return depth;
0177 }
0178 
0179 int HcalDDDSimConstants::getDepthEta29(const int& phi, const int& zside, const int& i) const {
0180   int depth = (i == 0) ? ldmap_.getMaxDepthLastHE(2, phi, zside) : -1;
0181   if (depth < 0)
0182     depth = (i == 1) ? depthEta29[1] : depthEta29[0];
0183   return depth;
0184 }
0185 
0186 int HcalDDDSimConstants::getDepthEta29M(const int& i, const bool& planOne) const {
0187   int depth = (i == 1) ? depthEta29[1] : depthEta29[0];
0188   if (i == 0 && planOne) {
0189     std::vector<int> phis;
0190     int detsp = ldmap_.validDet(phis);
0191     if (detsp == 2) {
0192       int zside = (phis[0] > 0) ? 1 : -1;
0193       int iphi = (phis[0] > 0) ? phis[0] : -phis[0];
0194       int depthsp = ldmap_.getMaxDepthLastHE(2, iphi, zside);
0195       if (depthsp > depth)
0196         depth = depthsp;
0197     }
0198   }
0199   return depth;
0200 }
0201 
0202 std::pair<int, double> HcalDDDSimConstants::getDetEta(const double& eta, const int& depth) const {
0203   int hsubdet(0), ieta(0);
0204   double etaR(0);
0205   double heta = fabs(eta);
0206   for (int i = 0; i < nEta; i++)
0207     if (heta > hpar->etaTable[i])
0208       ieta = i + 1;
0209   if (heta <= hpar->etaRange[1]) {
0210     if (((ieta == hpar->etaMin[1] && depth == depthEta16[1]) || (ieta > hpar->etaMax[0])) &&
0211         (ieta <= hpar->etaMax[1])) {
0212       hsubdet = static_cast<int>(HcalEndcap);
0213     } else {
0214       hsubdet = static_cast<int>(HcalBarrel);
0215     }
0216     etaR = eta;
0217   } else {
0218     hsubdet = static_cast<int>(HcalForward);
0219     double theta = 2. * atan(exp(-heta));
0220     double hR = zVcal * tan(theta);
0221     etaR = (eta >= 0. ? hR : -hR);
0222   }
0223   return std::pair<int, double>(hsubdet, etaR);
0224 }
0225 
0226 int HcalDDDSimConstants::getEta(const int& det, const int& lay, const double& hetaR) const {
0227   int ieta(0);
0228   if (det == static_cast<int>(HcalForward)) {  // Forward HCal
0229     ieta = hpar->etaMax[2];
0230     for (int i = nR - 1; i > 0; i--)
0231       if (hetaR < hpar->rTable[i])
0232         ieta = hpar->etaMin[2] + nR - i - 1;
0233   } else {  // Barrel or Endcap
0234     ieta = 1;
0235     for (int i = 0; i < nEta - 1; i++)
0236       if (hetaR > hpar->etaTable[i])
0237         ieta = i + 1;
0238     if (det == static_cast<int>(HcalBarrel)) {
0239       if (ieta > hpar->etaMax[0])
0240         ieta = hpar->etaMax[0];
0241       if (lay == maxLayer_) {
0242         if (hetaR > etaHO[1] && ieta == hpar->noff[2])
0243           ieta++;
0244       }
0245     } else if (det == static_cast<int>(HcalEndcap)) {
0246       if (ieta <= hpar->etaMin[1])
0247         ieta = hpar->etaMin[1];
0248     }
0249   }
0250   return ieta;
0251 }
0252 
0253 std::pair<int, int> HcalDDDSimConstants::getEtaDepth(
0254     const int& det, int etaR, const int& phi, const int& zside, int depth, const int& lay) const {
0255 #ifdef EDM_ML_DEBUG
0256   edm::LogVerbatim("HCalGeom") << "HcalDDDEsimConstants:getEtaDepth: I/P " << det << ":" << etaR << ":" << phi << ":"
0257                                << zside << ":" << depth << ":" << lay;
0258 #endif
0259   //Modify the depth index
0260   if ((det == static_cast<int>(HcalEndcap)) && (etaR == 17) && (lay == 1))
0261     etaR = 18;
0262   if (det == static_cast<int>(HcalForward)) {  // Forward HCal
0263   } else if (det == static_cast<int>(HcalOuter)) {
0264     depth = 4;
0265   } else {
0266     if (lay >= 0) {
0267       depth = layerGroup(det, etaR, phi, zside, lay - 1);
0268       if (((det == 2) && (etaR == 18)) || ((det == 1) && (etaR == 15)))
0269         if (etaR == hpar->noff[0] && lay > 1) {
0270           int kphi = phi + static_cast<int>((hpar->phioff[3] + 0.1) / hpar->phibin[etaR - 1]);
0271           kphi = (kphi - 1) % 4 + 1;
0272           if (kphi == 2 || kphi == 3)
0273             depth = layerGroup(det, etaR, phi, zside, lay - 2);
0274         }
0275     } else if (det == static_cast<int>(HcalBarrel)) {
0276       if (depth > getMaxDepth(det, etaR, phi, zside, false))
0277         depth = getMaxDepth(det, etaR, phi, zside, false);
0278     }
0279     if (etaR >= hpar->noff[1] && depth > getDepthEta29(phi, zside, 0)) {
0280       etaR -= getDepthEta29(phi, zside, 1);
0281     } else if (etaR == hpar->etaMin[1]) {
0282       if (det == static_cast<int>(HcalBarrel)) {
0283         if (depth > getDepthEta16(det, phi, zside))
0284           depth = getDepthEta16(det, phi, zside);
0285       } else {
0286         if (depth < getDepthEta16(det, phi, zside))
0287           depth = getDepthEta16(det, phi, zside);
0288       }
0289     }
0290   }
0291 #ifdef EDM_ML_DEBUG
0292   edm::LogVerbatim("HCalGeom") << "HcalDDDEsimConstants:getEtaDepth: O/P " << etaR << ":" << depth;
0293 #endif
0294   return std::pair<int, int>(etaR, depth);
0295 }
0296 
0297 double HcalDDDSimConstants::getEtaHO(const double& etaR, const double& x, const double& y, const double& z) const {
0298   if (hpar->zHO.size() > 4) {
0299     double eta = fabs(etaR);
0300     double r = std::sqrt(x * x + y * y);
0301     if (r > rminHO) {
0302       double zz = fabs(z);
0303       if (zz > hpar->zHO[3]) {
0304         if (eta <= hpar->etaTable[10])
0305           eta = hpar->etaTable[10] + 0.001;
0306       } else if (zz > hpar->zHO[1]) {
0307         if (eta <= hpar->etaTable[4])
0308           eta = hpar->etaTable[4] + 0.001;
0309       }
0310     }
0311     eta = (z >= 0. ? eta : -eta);
0312 #ifdef EDM_ML_DEBUG
0313     std::string chk = (eta != etaR) ? " **** Check *****" : " ";
0314     edm::LogVerbatim("HCalGeom") << "R " << r << " Z " << z << " eta " << etaR << ":" << eta << chk;
0315 #endif
0316     return eta;
0317   } else {
0318     return etaR;
0319   }
0320 }
0321 
0322 int HcalDDDSimConstants::getFrontLayer(const int& det, const int& eta) const {
0323   int lay = 0;
0324   if (det == 1) {
0325     if (std::abs(eta) == 16)
0326       lay = layFHB[1];
0327     else
0328       lay = layFHB[0];
0329   } else {
0330     if (std::abs(eta) == 16)
0331       lay = layFHE[1];
0332     else if (std::abs(eta) == 18)
0333       lay = layFHE[2];
0334     else
0335       lay = layFHE[0];
0336   }
0337   return lay;
0338 }
0339 
0340 int HcalDDDSimConstants::getLastLayer(const int& det, const int& eta) const {
0341   int lay = 0;
0342   if (det == 1) {
0343     if (std::abs(eta) == 15)
0344       lay = layBHB[1];
0345     else if (std::abs(eta) == 16)
0346       lay = layBHB[2];
0347     else
0348       lay = layBHB[0];
0349   } else {
0350     if (std::abs(eta) == 16)
0351       lay = layBHE[1];
0352     else if (std::abs(eta) == 17)
0353       lay = layBHE[2];
0354     else if (std::abs(eta) == 18)
0355       lay = layBHE[3];
0356     else
0357       lay = layBHE[0];
0358   }
0359   return lay;
0360 }
0361 
0362 double HcalDDDSimConstants::getLayer0Wt(const int& det, const int& phi, const int& zside) const {
0363   double wt = ldmap_.getLayer0Wt(det, phi, zside);
0364   if (wt < 0)
0365     wt = (det == 2) ? hpar->Layer0Wt[1] : hpar->Layer0Wt[0];
0366   return wt;
0367 }
0368 
0369 int HcalDDDSimConstants::getLayerFront(
0370     const int& det, const int& eta, const int& phi, const int& zside, const int& depth) const {
0371   int layer = ldmap_.getLayerFront(det, eta, phi, zside, depth);
0372   if (layer < 0) {
0373     if (det == 1 || det == 2) {
0374       layer = 1;
0375       for (int l = 0; l < getLayerMax(eta, depth); ++l) {
0376         if (static_cast<int>(layerGroup(eta - 1, l)) == depth) {
0377           layer = l + 1;
0378           break;
0379         }
0380       }
0381     } else {
0382       layer = (eta > hpar->noff[2]) ? maxLayerHB_ + 1 : maxLayer_;
0383     }
0384   }
0385   return layer;
0386 }
0387 
0388 int HcalDDDSimConstants::getLayerBack(
0389     const int& det, const int& eta, const int& phi, const int& zside, const int& depth) const {
0390   int layer = ldmap_.getLayerBack(det, eta, phi, zside, depth);
0391   if (layer < 0) {
0392     if (det == 1 || det == 2) {
0393       layer = depths[depth - 1][eta - 1];
0394     } else {
0395       layer = maxLayer_;
0396     }
0397   }
0398   return layer;
0399 }
0400 
0401 int HcalDDDSimConstants::getLayerMax(const int& eta, const int& depth) const {
0402   int layermx = ((eta < hpar->etaMin[1]) && depth - 1 < maxDepth[0]) ? maxLayerHB_ + 1
0403                                                                      : static_cast<int>(layerGroupSize(eta - 1));
0404   return layermx;
0405 }
0406 
0407 int HcalDDDSimConstants::getMaxDepth(
0408     const int& det, const int& eta, const int& phi, const int& zside, const bool& partialDetOnly) const {
0409   int dmax(-1);
0410   if (partialDetOnly) {
0411     if (ldmap_.isValid(det, phi, zside)) {
0412       dmax = ldmap_.getDepths(eta).second;
0413     }
0414   } else if (det == 1 || det == 2) {
0415     if (ldmap_.isValid(det, phi, zside))
0416       dmax = ldmap_.getDepths(eta).second;
0417     else if (det == 2)
0418       dmax = (maxDepth[1] > 0) ? layerGroup(eta - 1, maxLayer_) : 0;
0419     else if (eta == hpar->etaMax[0])
0420       dmax = getDepthEta16(det, phi, zside);
0421     else
0422       dmax = layerGroup(eta - 1, maxLayerHB_);
0423   } else if (det == 3) {  // HF
0424     dmax = maxHFDepth(zside * eta, phi);
0425   } else if (det == 4) {  // HO
0426     dmax = maxDepth[3];
0427   } else {
0428     dmax = -1;
0429   }
0430   return dmax;
0431 }
0432 
0433 int HcalDDDSimConstants::getMinDepth(
0434     const int& det, const int& eta, const int& phi, const int& zside, const bool& partialDetOnly) const {
0435   int lmin(-1);
0436   if (partialDetOnly) {
0437     if (ldmap_.isValid(det, phi, zside)) {
0438       lmin = ldmap_.getDepths(eta).first;
0439     }
0440   } else if (det == 3) {  // HF
0441     lmin = 1;
0442   } else if (det == 4) {  // HO
0443     lmin = maxDepth[3];
0444   } else {
0445     if (ldmap_.isValid(det, phi, zside)) {
0446       lmin = ldmap_.getDepths(eta).first;
0447     } else if (layerGroupSize(eta - 1) > 0) {
0448       lmin = static_cast<int>(layerGroup(eta - 1, 0));
0449       unsigned int type = (det == 1) ? 0 : 1;
0450       if (type == 1 && eta == hpar->etaMin[1])
0451         lmin = getDepthEta16(det, phi, zside);
0452     } else {
0453       lmin = 1;
0454     }
0455   }
0456   return lmin;
0457 }
0458 
0459 std::pair<int, int> HcalDDDSimConstants::getModHalfHBHE(const int& type) const {
0460   if (type == 0) {
0461     return std::pair<int, int>(nmodHB, nzHB);
0462   } else {
0463     return std::pair<int, int>(nmodHE, nzHE);
0464   }
0465 }
0466 
0467 std::pair<double, double> HcalDDDSimConstants::getPhiCons(const int& det, const int& ieta) const {
0468   double fioff(0), fibin(0);
0469   if (det == static_cast<int>(HcalForward)) {  // Forward HCal
0470     fioff = hpar->phioff[2];
0471     fibin = hpar->phitable[ieta - hpar->etaMin[2]];
0472     if (unitPhi(fibin) > 2) {  // HF double-phi
0473       fioff = hpar->phioff[4];
0474     }
0475   } else {  // Barrel or Endcap
0476     if (det == static_cast<int>(HcalEndcap)) {
0477       fioff = hpar->phioff[1];
0478     } else {
0479       fioff = hpar->phioff[0];
0480     }
0481     fibin = hpar->phibin[ieta - 1];
0482   }
0483   return std::pair<double, double>(fioff, fibin);
0484 }
0485 
0486 std::vector<std::pair<int, double> > HcalDDDSimConstants::getPhis(const int& subdet, const int& ieta) const {
0487   std::vector<std::pair<int, double> > phis;
0488   int ietaAbs = (ieta > 0) ? ieta : -ieta;
0489   std::pair<double, double> ficons = getPhiCons(subdet, ietaAbs);
0490   int nphi = int((2._pi + 0.1 * ficons.second) / ficons.second);
0491   int units = unitPhi(subdet, ietaAbs);
0492   for (int ifi = 0; ifi < nphi; ++ifi) {
0493     double phi = -ficons.first + (ifi + 0.5) * ficons.second;
0494     int iphi = phiNumber(ifi + 1, units);
0495     phis.emplace_back(std::pair<int, double>(iphi, phi));
0496   }
0497 #ifdef EDM_ML_DEBUG
0498   edm::LogVerbatim("HCalGeom") << "getPhis: subdet|ieta|iphi " << subdet << "|" << ieta << " with " << phis.size()
0499                                << " phi bins";
0500   for (unsigned int k = 0; k < phis.size(); ++k)
0501     edm::LogVerbatim("HCalGeom") << "[" << k << "] iphi " << phis[k].first << " phi "
0502                                  << convertRadToDeg(phis[k].second);
0503 #endif
0504   return phis;
0505 }
0506 
0507 std::vector<HcalCellType> HcalDDDSimConstants::HcalCellTypes() const {
0508   std::vector<HcalCellType> cellTypes = HcalCellTypes(HcalBarrel);
0509 #ifdef EDM_ML_DEBUG
0510   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: " << cellTypes.size() << " cells of type HCal Barrel";
0511   for (unsigned int i = 0; i < cellTypes.size(); i++)
0512     edm::LogVerbatim("HCalGeom") << "Cell " << i << " " << cellTypes[i];
0513 #endif
0514 
0515   std::vector<HcalCellType> hoCells = HcalCellTypes(HcalOuter);
0516 #ifdef EDM_ML_DEBUG
0517   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: " << hoCells.size() << " cells of type HCal Outer";
0518   for (unsigned int i = 0; i < hoCells.size(); i++)
0519     edm::LogVerbatim("HCalGeom") << "Cell " << i << " " << hoCells[i];
0520 #endif
0521   cellTypes.insert(cellTypes.end(), hoCells.begin(), hoCells.end());
0522 
0523   std::vector<HcalCellType> heCells = HcalCellTypes(HcalEndcap);
0524 #ifdef EDM_ML_DEBUG
0525   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: " << heCells.size() << " cells of type HCal Endcap";
0526   for (unsigned int i = 0; i < heCells.size(); i++)
0527     edm::LogVerbatim("HCalGeom") << "Cell " << i << " " << heCells[i];
0528 #endif
0529   cellTypes.insert(cellTypes.end(), heCells.begin(), heCells.end());
0530 
0531   std::vector<HcalCellType> hfCells = HcalCellTypes(HcalForward);
0532 #ifdef EDM_ML_DEBUG
0533   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants: " << hfCells.size() << " cells of type HCal Forward";
0534   for (unsigned int i = 0; i < hfCells.size(); i++)
0535     edm::LogVerbatim("HCalGeom") << "Cell " << i << " " << hfCells[i];
0536 #endif
0537   cellTypes.insert(cellTypes.end(), hfCells.begin(), hfCells.end());
0538 
0539   return cellTypes;
0540 }
0541 
0542 std::vector<HcalCellType> HcalDDDSimConstants::HcalCellTypes(const HcalSubdetector& subdet,
0543                                                              int ieta,
0544                                                              int depthl) const {
0545   std::vector<HcalCellType> cellTypes;
0546   if (subdet == HcalForward) {
0547     if (dzVcal < 0)
0548       return cellTypes;
0549   }
0550 
0551   int dmin, dmax, indx, nz;
0552   double hsize = 0;
0553   switch (subdet) {
0554     case HcalEndcap:
0555       dmin = 1;
0556       dmax = (maxDepth[1] > 0) ? maxLayer_ + 1 : 0;
0557       indx = 1;
0558       nz = nzHE;
0559       break;
0560     case HcalForward:
0561       dmin = 1;
0562       dmax = (!idHF2QIE.empty()) ? 2 : maxDepth[2];
0563       indx = 2;
0564       nz = 2;
0565       break;
0566     case HcalOuter:
0567       dmin = 4;
0568       dmax = 4;
0569       indx = 0;
0570       nz = nzHB;
0571       break;
0572     default:
0573       dmin = 1;
0574       dmax = maxLayerHB_ + 1;
0575       indx = 0;
0576       nz = nzHB;
0577       break;
0578   }
0579   if (depthl > 0)
0580     dmin = dmax = depthl;
0581   int ietamin = (ieta > 0) ? ieta : hpar->etaMin[indx];
0582   int ietamax = (ieta > 0) ? ieta : hpar->etaMax[indx];
0583   int phi = (indx == 2) ? 3 : 1;
0584 
0585   // Get the Cells
0586   int subdet0 = static_cast<int>(subdet);
0587   for (int depth = dmin; depth <= dmax; depth++) {
0588     int shift = getShift(subdet, depth);
0589     double gain = getGain(subdet, depth);
0590     if (subdet == HcalForward) {
0591       if (depth % 2 == 1)
0592         hsize = dzVcal;
0593       else
0594         hsize = dzVcal - 0.5 * dlShort;
0595     }
0596     for (int eta = ietamin; eta <= ietamax; eta++) {
0597       for (int iz = 0; iz < nz; ++iz) {
0598         int zside = -2 * iz + 1;
0599         HcalCellType::HcalCell temp1 = cell(subdet0, zside, depth, eta, phi);
0600         if (temp1.ok) {
0601           std::vector<std::pair<int, double> > phis = getPhis(subdet0, eta);
0602           HcalCellType temp2(subdet, eta, zside, depth, temp1, shift, gain, hsize);
0603           double dphi, fioff;
0604           std::vector<int> phiMiss2;
0605           if ((subdet == HcalBarrel) || (subdet == HcalOuter)) {
0606             fioff = hpar->phioff[0];
0607             dphi = hpar->phibin[eta - 1];
0608             if (subdet == HcalOuter) {
0609               if (eta == hpar->noff[4]) {
0610                 int kk = (iz == 0) ? 7 : (7 + hpar->noff[5]);
0611                 for (int miss = 0; miss < hpar->noff[5 + iz]; miss++) {
0612                   phiMiss2.emplace_back(hpar->noff[kk]);
0613                   kk++;
0614                 }
0615               }
0616             }
0617           } else if (subdet == HcalEndcap) {
0618             fioff = hpar->phioff[1];
0619             dphi = hpar->phibin[eta - 1];
0620           } else {
0621             fioff = hpar->phioff[2];
0622             dphi = hpar->phitable[eta - hpar->etaMin[2]];
0623             if (unitPhi(dphi) > 2)
0624               fioff = hpar->phioff[4];
0625           }
0626           int unit = unitPhi(dphi);
0627           temp2.setPhi(phis, phiMiss2, fioff, dphi, unit);
0628           cellTypes.emplace_back(temp2);
0629           // For HF look at extra cells
0630           if ((subdet == HcalForward) && (!idHF2QIE.empty())) {
0631             HcalCellType temp3(subdet, eta, zside + 2, depth, temp1, shift, gain, hsize);
0632             std::vector<int> phiMiss3;
0633             for (auto& phi : phis) {
0634               bool ok(false);
0635               for (auto l : idHF2QIE) {
0636                 if ((eta * zside == l.ieta()) && (phi.first == l.iphi())) {
0637                   ok = true;
0638                   break;
0639                 }
0640               }
0641               if (!ok)
0642                 phiMiss3.emplace_back(phi.first);
0643             }
0644             dphi = hpar->phitable[eta - hpar->etaMin[2]];
0645             unit = unitPhi(dphi);
0646             fioff = (unit > 2) ? hpar->phioff[4] : hpar->phioff[2];
0647             temp3.setPhi(phis, phiMiss2, fioff, dphi, unit);
0648             cellTypes.emplace_back(temp3);
0649           }
0650         }
0651       }
0652     }
0653   }
0654   return cellTypes;
0655 }
0656 
0657 int HcalDDDSimConstants::maxHFDepth(const int& eta, const int& iphi) const {
0658   int mxdepth = maxDepth[2];
0659   if (!idHF2QIE.empty()) {
0660     bool ok(false);
0661     for (auto k : idHF2QIE) {
0662       if ((eta == k.ieta()) && (iphi == k.iphi())) {
0663         ok = true;
0664         break;
0665       }
0666     }
0667     if (!ok)
0668       mxdepth = 2;
0669   }
0670   return mxdepth;
0671 }
0672 
0673 unsigned int HcalDDDSimConstants::numberOfCells(const HcalSubdetector& subdet) const {
0674   unsigned int num = 0;
0675   std::vector<HcalCellType> cellTypes = HcalCellTypes(subdet);
0676   for (auto& cellType : cellTypes) {
0677     num += static_cast<unsigned int>(cellType.nPhiBins());
0678   }
0679 #ifdef EDM_ML_DEBUG
0680   edm::LogVerbatim("HCalGeom") << "HcalDDDSimConstants:numberOfCells " << cellTypes.size() << " " << num
0681                                << " for subdetector " << subdet;
0682 #endif
0683   return num;
0684 }
0685 
0686 int HcalDDDSimConstants::phiNumber(const int& phi, const int& units) const {
0687   int iphi_skip = phi;
0688   if (units == 2)
0689     iphi_skip = (phi - 1) * 2 + 1;
0690   else if (units == 4)
0691     iphi_skip = (phi - 1) * 4 - 1;
0692   if (iphi_skip < 0)
0693     iphi_skip += 72;
0694   return iphi_skip;
0695 }
0696 
0697 void HcalDDDSimConstants::printTiles() const {
0698   std::vector<int> phis;
0699   int detsp = ldmap_.validDet(phis);
0700   int kphi = (detsp > 0) ? phis[0] : 1;
0701   int zside = (kphi > 0) ? 1 : -1;
0702   int iphi = (kphi > 0) ? kphi : -kphi;
0703   edm::LogVerbatim("HCalGeom") << "Tile Information for HB from " << hpar->etaMin[0] << " to " << hpar->etaMax[0];
0704   for (int eta = hpar->etaMin[0]; eta <= hpar->etaMax[0]; eta++) {
0705     int dmax = getMaxDepth(1, eta, iphi, -zside, false);
0706     for (int depth = 1; depth <= dmax; depth++)
0707       printTileHB(eta, iphi, -zside, depth);
0708     if (detsp == 1) {
0709       int dmax = getMaxDepth(1, eta, iphi, zside, false);
0710       for (int depth = 1; depth <= dmax; depth++)
0711         printTileHB(eta, iphi, zside, depth);
0712     }
0713   }
0714 
0715   edm::LogVerbatim("HCalGeom") << "\nTile Information for HE from " << hpar->etaMin[1] << " to " << hpar->etaMax[1];
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 (static_cast<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 (static_cast<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 (static_cast<int>(hpar->noff.size()) > (noffsize + 4)) {
0849       noffl += (2 * hpar->noff[noffsize + 4]);
0850       if (static_cast<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 (static_cast<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 (static_cast<int>(hpar->noff.size()) > (noffsize + 5)) {
0900     noffk += (2 * hpar->noff[noffsize + 4]);
0901     if (static_cast<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 (static_cast<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 (static_cast<int>(hpar->noff.size()) > noffm) {
0938         int ndnext = hpar->noff[noffm];
0939         if (ndnext > 4 && static_cast<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 && static_cast<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 == static_cast<unsigned int>(eta + 1)) {
1207       return it.layerGroup.size();
1208     }
1209     if (it.layer > static_cast<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 == static_cast<unsigned int>(eta + 1)) {
1220       return it.layerGroup.at(i);
1221     }
1222     if (it.layer > static_cast<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) ? static_cast<unsigned int>(depth0) : layerGroup(eta - 1, lay);
1232   return depth;
1233 }