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