File indexing completed on 2024-04-06 12:14:47
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
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)) {
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 {
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
0260 if ((det == static_cast<int>(HcalEndcap)) && (etaR == 17) && (lay == 1))
0261 etaR = 18;
0262 if (det == static_cast<int>(HcalForward)) {
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) {
0424 dmax = maxHFDepth(zside * eta, phi);
0425 } else if (det == 4) {
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) {
0441 lmin = 1;
0442 } else if (det == 4) {
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)) {
0470 fioff = hpar->phioff[2];
0471 fibin = hpar->phitable[ieta - hpar->etaMin[2]];
0472 if (unitPhi(fibin) > 2) {
0473 fioff = hpar->phioff[4];
0474 }
0475 } else {
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
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
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
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
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 }
