File indexing completed on 2021-12-17 23:36:21
0001 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "DataFormats/Math/interface/GeantUnits.h"
0006 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0007 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0008 #include <algorithm>
0009 #include <cmath>
0010
0011
0012 using namespace geant_units::operators;
0013
0014 enum { kHOSizePreLS1 = 2160, kHFSizePreLS1 = 1728 };
0015
0016 HcalDDDRecConstants::HcalDDDRecConstants(const HcalParameters* hp, const HcalDDDSimConstants& hc)
0017 : hpar(hp), hcons(hc) {
0018 #ifdef EDM_ML_DEBUG
0019 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants::HcalDDDRecConstants (const HcalParameters* hp) constructor";
0020 #endif
0021 initialize();
0022 }
0023
0024 HcalDDDRecConstants::~HcalDDDRecConstants() {
0025 #ifdef EDM_ML_DEBUG
0026 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants::destructed!!!";
0027 #endif
0028 }
0029
0030 std::vector<int> HcalDDDRecConstants::getDepth(const unsigned int& eta, const bool& extra) const {
0031 if (!extra) {
0032 std::vector<HcalParameters::LayerItem>::const_iterator last = hpar->layerGroupEtaRec.begin();
0033 for (std::vector<HcalParameters::LayerItem>::const_iterator it = hpar->layerGroupEtaRec.begin();
0034 it != hpar->layerGroupEtaRec.end();
0035 ++it) {
0036 if (it->layer == eta + 1)
0037 return it->layerGroup;
0038 if (it->layer > eta + 1)
0039 return last->layerGroup;
0040 last = it;
0041 }
0042 return last->layerGroup;
0043 } else {
0044 std::map<int, int> layers;
0045 hcons.ldMap()->getLayerDepth(eta + 1, layers);
0046 std::vector<int> depths;
0047 for (unsigned int lay = 0; lay < layers.size(); ++lay)
0048 depths.emplace_back(layers[lay + 1]);
0049 return depths;
0050 }
0051 }
0052
0053 std::vector<int> HcalDDDRecConstants::getDepth(const int& det,
0054 const int& phi,
0055 const int& zside,
0056 const unsigned int& eta) const {
0057 std::map<int, int> layers;
0058 hcons.ldMap()->getLayerDepth(det, eta + 1, phi, zside, layers);
0059 if (layers.empty()) {
0060 return getDepth(eta, false);
0061 } else {
0062 std::vector<int> depths;
0063 for (unsigned int lay = 0; lay < layers.size(); ++lay)
0064 depths.emplace_back(layers[lay + 1]);
0065 return depths;
0066 }
0067 }
0068
0069 std::vector<HcalDDDRecConstants::HcalEtaBin> HcalDDDRecConstants::getEtaBins(const int& itype) const {
0070 std::vector<HcalDDDRecConstants::HcalEtaBin> bins;
0071 unsigned int type = (itype == 0) ? 0 : 1;
0072 HcalSubdetector subdet = HcalSubdetector(type + 1);
0073 std::vector<int> phiSp;
0074 HcalSubdetector subdetSp = HcalSubdetector(hcons.ldMap()->validDet(phiSp));
0075 std::map<int, int> layers;
0076 for (int iz = 0; iz < 2; ++iz) {
0077 int zside = 2 * iz - 1;
0078 for (int ieta = iEtaMin[type]; ieta <= iEtaMax[type]; ++ieta) {
0079 std::vector<std::pair<int, double>> phis = getPhis(subdet, ieta);
0080 std::vector<std::pair<int, double>> phiUse;
0081 getLayerDepth(ieta, layers);
0082 if (subdet == subdetSp) {
0083 for (auto& phi : phis) {
0084 if (std::find(phiSp.begin(), phiSp.end(), (zside * phi.first)) == phiSp.end()) {
0085 phiUse.emplace_back(phi);
0086 }
0087 }
0088 } else {
0089 phiUse.insert(phiUse.end(), phis.begin(), phis.end());
0090 }
0091 if (!phiUse.empty())
0092 getOneEtaBin(subdet, ieta, zside, phiUse, layers, false, bins);
0093 }
0094 }
0095 if (subdetSp == subdet) {
0096 for (int ieta = iEtaMin[type]; ieta <= iEtaMax[type]; ++ieta) {
0097 std::vector<std::pair<int, double>> phis = getPhis(subdet, ieta);
0098 for (int iz = 0; iz < 2; ++iz) {
0099 int zside = 2 * iz - 1;
0100 std::vector<std::pair<int, double>> phiUse;
0101 for (int i : phiSp) {
0102 for (auto& phi : phis) {
0103 if (i == zside * phi.first) {
0104 phiUse.emplace_back(phi);
0105 break;
0106 }
0107 }
0108 }
0109 if (!phiUse.empty()) {
0110 hcons.ldMap()->getLayerDepth(subdet, ieta, phiUse[0].first, zside, layers);
0111 getOneEtaBin(subdet, ieta, zside, phiUse, layers, true, bins);
0112 }
0113 }
0114 }
0115 }
0116 #ifdef EDM_ML_DEBUG
0117 edm::LogVerbatim("HCalGeom") << "Prepares " << bins.size() << " eta bins for type " << type;
0118 for (unsigned int i = 0; i < bins.size(); ++i) {
0119 edm::LogVerbatim("HCalGeom") << "Bin[" << i << "]: Eta = (" << bins[i].ieta << ":" << bins[i].etaMin << ":"
0120 << bins[i].etaMax << "), Zside = " << bins[i].zside << ", phis = ("
0121 << bins[i].phis.size() << ":" << bins[i].dphi << ") and " << bins[i].layer.size()
0122 << " depths (start) " << bins[i].depthStart;
0123 for (unsigned int k = 0; k < bins[i].layer.size(); ++k)
0124 edm::LogVerbatim("HCalGeom") << " [" << k << "] " << bins[i].layer[k].first << ":" << bins[i].layer[k].second;
0125 edm::LogVerbatim("HCalGeom") << "Phi sets";
0126 for (unsigned int k = 0; k < bins[i].phis.size(); ++k)
0127 edm::LogVerbatim("HCalGeom") << "[" << k << "] " << bins[i].phis[k].first << ":" << bins[i].phis[k].second;
0128 }
0129 #endif
0130 return bins;
0131 }
0132
0133 std::pair<double, double> HcalDDDRecConstants::getEtaPhi(const int& subdet, const int& ieta, const int& iphi) const {
0134 int ietaAbs = (ieta > 0) ? ieta : -ieta;
0135 double eta(0), phi(0);
0136 if ((subdet == static_cast<int>(HcalBarrel)) || (subdet == static_cast<int>(HcalEndcap)) ||
0137 (subdet == static_cast<int>(HcalOuter))) {
0138 int unit = hcons.unitPhi(phibin[ietaAbs - 1]);
0139 int kphi = (unit == 2) ? ((iphi - 1) / 2 + 1) : iphi;
0140 double foff = (ietaAbs <= iEtaMax[0]) ? hpar->phioff[0] : hpar->phioff[1];
0141 eta = 0.5 * (etaTable[ietaAbs - 1] + etaTable[ietaAbs]);
0142 phi = foff + (kphi - 0.5) * phibin[ietaAbs - 1];
0143 } else {
0144 ietaAbs -= iEtaMin[2];
0145 int unit = hcons.unitPhi(hpar->phitable[ietaAbs - 1]);
0146 int kphi = (unit == 4) ? ((iphi - 3) / 4 + 1) : ((iphi - 1) / 2 + 1);
0147 double foff = (unit > 2) ? hpar->phioff[4] : hpar->phioff[2];
0148 eta = 0.5 * (hpar->etaTableHF[ietaAbs - 1] + hpar->etaTableHF[ietaAbs]);
0149 phi = foff + (kphi - 0.5) * hpar->phitable[ietaAbs - 1];
0150 }
0151 if (ieta < 0)
0152 eta = -eta;
0153 if (phi > M_PI)
0154 phi -= (2 * M_PI);
0155 #ifdef EDM_ML_DEBUG
0156 edm::LogVerbatim("HCalGeom") << "getEtaPhi: subdet|ieta|iphi " << subdet << "|" << ieta << "|" << iphi << " eta|phi "
0157 << eta << "|" << phi;
0158 #endif
0159 return std::pair<double, double>(eta, phi);
0160 }
0161
0162 HcalDDDRecConstants::HcalID HcalDDDRecConstants::getHCID(int subdet, int keta, int iphi, int lay, int idepth) const {
0163 int ieta = (keta > 0) ? keta : -keta;
0164 int zside = (keta > 0) ? 1 : -1;
0165 int eta(ieta), phi(iphi), depth(idepth);
0166 if ((subdet == static_cast<int>(HcalOuter)) ||
0167 ((subdet == static_cast<int>(HcalBarrel)) && (lay > maxLayerHB_ + 1))) {
0168 subdet = static_cast<int>(HcalOuter);
0169 depth = 4;
0170 } else if (subdet == static_cast<int>(HcalBarrel) || subdet == static_cast<int>(HcalEndcap)) {
0171 eta = ietaMap[ieta - 1];
0172 int unit = phiUnitS[ieta - 1];
0173 int phi0 = (iphi - 1) / (hpar->phigroup[eta - 1]);
0174 if (unit == 2) {
0175 phi0 = (iphi + 1) / 2;
0176 phi0 = (phi0 - 1) / (hpar->phigroup[eta - 1]);
0177 } else if (unit == 4) {
0178 phi0 = (iphi + 1) / 4;
0179 phi0 = (phi0 - 1) / (hpar->phigroup[eta - 1]);
0180 }
0181 ++phi0;
0182 unit = hcons.unitPhi(phibin[eta - 1]);
0183 phi = hcons.phiNumber(phi0, unit);
0184 depth = hcons.findDepth(subdet, eta, phi, zside, lay - 1);
0185 if (depth <= 0)
0186 depth = layerGroup(eta - 1, lay - 1);
0187 if (eta == iEtaMin[1]) {
0188 if (subdet == static_cast<int>(HcalBarrel)) {
0189 if (depth > hcons.getDepthEta16(subdet, phi, zside))
0190 depth = hcons.getDepthEta16(subdet, phi, zside);
0191 } else {
0192 if (depth < hcons.getDepthEta16(subdet, phi, zside))
0193 depth = hcons.getDepthEta16(subdet, phi, zside);
0194 }
0195 } else if (eta == hpar->noff[0] && lay > 1) {
0196 int kphi = phi + int((hpar->phioff[3] + 0.1) / phibin[eta - 1]);
0197 kphi = (kphi - 1) % 4 + 1;
0198 if (kphi == 2 || kphi == 3)
0199 depth = layerGroup(eta - 1, lay - 2);
0200 } else if (eta == hpar->noff[1] && depth > hcons.getDepthEta29(phi, zside, 0)) {
0201 eta -= hcons.getDepthEta29(phi, zside, 1);
0202 }
0203 }
0204 #ifdef EDM_ML_DEBUG
0205 edm::LogVerbatim("HCalGeom") << "getHCID: input " << subdet << ":" << ieta << ":" << iphi << ":" << idepth << ":"
0206 << lay << " output " << eta << ":" << phi << ":" << depth;
0207 #endif
0208 return HcalDDDRecConstants::HcalID(subdet, eta, phi, depth);
0209 }
0210
0211 std::vector<HcalDDDRecConstants::HFCellParameters> HcalDDDRecConstants::getHFCellParameters() const {
0212 std::vector<HcalDDDRecConstants::HFCellParameters> cells;
0213 unsigned int nEta = hcons.getPhiTableHF().size();
0214 if (maxDepth[2] > 0) {
0215 for (unsigned int k = 0; k < nEta; ++k) {
0216 int ieta = iEtaMin[2] + k;
0217 int dphi = (int)(0.001 + hcons.getPhiTableHF()[k] / (5._deg));
0218 int iphi = (dphi == 4) ? 3 : 1;
0219 int nphi = 72 / dphi;
0220 double rMin = hcons.getRTableHF()[nEta - k - 1] / CLHEP::cm;
0221 double rMax = hcons.getRTableHF()[nEta - k] / CLHEP::cm;
0222 HcalDDDRecConstants::HFCellParameters cell1(ieta, 1, iphi, dphi, nphi, rMin, rMax);
0223 cells.emplace_back(cell1);
0224 HcalDDDRecConstants::HFCellParameters cell2(-ieta, 1, iphi, dphi, nphi, rMin, rMax);
0225 cells.emplace_back(cell2);
0226 }
0227 }
0228 if (maxDepth[2] > 2) {
0229 if (!hcons.getIdHF2QIE().empty()) {
0230 for (unsigned int k = 0; k < hcons.getIdHF2QIE().size(); ++k) {
0231 int ieta = hcons.getIdHF2QIE()[k].ieta();
0232 int ind = std::abs(ieta) - iEtaMin[2];
0233 int dphi = (int)(0.001 + hcons.getPhiTableHF()[ind] / (5._deg));
0234 int iphi = hcons.getIdHF2QIE()[k].iphi();
0235 double rMin = hcons.getRTableHF()[nEta - ind - 1] / CLHEP::cm;
0236 double rMax = hcons.getRTableHF()[nEta - ind] / CLHEP::cm;
0237 HcalDDDRecConstants::HFCellParameters cell1(ieta, 3, iphi, dphi, 1, rMin, rMax);
0238 cells.emplace_back(cell1);
0239 }
0240 } else {
0241 for (unsigned int k = 0; k < nEta; ++k) {
0242 int ieta = iEtaMin[2] + k;
0243 int dphi = (int)(0.001 + hcons.getPhiTableHF()[k] / (5._deg));
0244 int iphi = (dphi == 4) ? 3 : 1;
0245 int nphi = 72 / dphi;
0246 double rMin = hcons.getRTableHF()[nEta - k - 1] / CLHEP::cm;
0247 double rMax = hcons.getRTableHF()[nEta - k] / CLHEP::cm;
0248 HcalDDDRecConstants::HFCellParameters cell1(ieta, 3, iphi, dphi, nphi, rMin, rMax);
0249 cells.emplace_back(cell1);
0250 HcalDDDRecConstants::HFCellParameters cell2(-ieta, 3, iphi, dphi, nphi, rMin, rMax);
0251 cells.emplace_back(cell2);
0252 }
0253 }
0254 }
0255 #ifdef EDM_ML_DEBUG
0256 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants returns " << cells.size() << " HF cell parameters";
0257 for (unsigned int k = 0; k < cells.size(); ++k)
0258 edm::LogVerbatim("HCalGeom") << "Cell[" << k << "] : (" << cells[k].ieta << ", " << cells[k].depth << ", "
0259 << cells[k].firstPhi << ", " << cells[k].stepPhi << ", " << cells[k].nPhi << ", "
0260 << cells[k].rMin << ", " << cells[k].rMax << ")";
0261 #endif
0262 return cells;
0263 }
0264
0265 void HcalDDDRecConstants::getLayerDepth(const int& ieta, std::map<int, int>& layers) const {
0266 layers.clear();
0267 for (unsigned int l = 0; l < layerGroupSize(ieta - 1); ++l) {
0268 int lay = l + 1;
0269 layers[lay] = layerGroup(ieta - 1, l);
0270 }
0271 #ifdef EDM_ML_DEBUG
0272 edm::LogVerbatim("HCalGeom") << "getLayerDepth::Input " << ieta << " Output " << layers.size() << " entries";
0273 for (std::map<int, int>::iterator itr = layers.begin(); itr != layers.end(); ++itr)
0274 edm::LogVerbatim("HCalGeom") << " [" << itr->first << "] " << itr->second;
0275 #endif
0276 }
0277
0278 int HcalDDDRecConstants::getLayerBack(const int& idet, const int& ieta, const int& iphi, const int& depth) const {
0279 int subdet = (idet == 1) ? 1 : 2;
0280 int zside = (ieta > 0) ? 1 : -1;
0281 int eta = zside * ieta;
0282 int layBack = hcons.ldMap()->getLayerBack(subdet, eta, iphi, zside, depth);
0283 int laymax = hcons.getLastLayer(subdet, ieta);
0284 if (layBack < 0 && eta <= hpar->etaMax[1]) {
0285 for (unsigned int k = 0; k < layerGroupSize(eta - 1); ++k) {
0286 if (depth + 1 == (int)layerGroup(eta - 1, k)) {
0287 layBack = k - 1;
0288 break;
0289 }
0290 }
0291 }
0292 if (layBack < 0 || layBack > laymax)
0293 layBack = laymax;
0294 #ifdef EDM_ML_DEBUG
0295 edm::LogVerbatim("HCalGeom") << "getLayerBack::Input " << idet << ":" << ieta << ":" << iphi << ":" << depth
0296 << " Output " << layBack;
0297 #endif
0298 return layBack;
0299 }
0300
0301 int HcalDDDRecConstants::getLayerFront(const int& idet, const int& ieta, const int& iphi, const int& depth) const {
0302 int subdet = (idet == 1) ? 1 : 2;
0303 int zside = (ieta > 0) ? 1 : -1;
0304 int eta = zside * ieta;
0305 int layFront = hcons.ldMap()->getLayerFront(subdet, eta, iphi, zside, depth);
0306 int laymin = hcons.getFrontLayer(subdet, ieta);
0307 if ((layFront < 0) || ((subdet == static_cast<int>(HcalEndcap)) && (eta == 16))) {
0308 if ((subdet == static_cast<int>(HcalEndcap)) && (eta == 16)) {
0309 layFront = laymin;
0310 } else if (eta <= hpar->etaMax[1]) {
0311 for (unsigned int k = 0; k < layerGroupSize(eta - 1); ++k) {
0312 if (depth == (int)layerGroup(eta - 1, k)) {
0313 if ((int)(k) >= laymin) {
0314 layFront = k;
0315 break;
0316 }
0317 }
0318 }
0319 }
0320 } else {
0321 if (layFront < laymin)
0322 layFront = laymin;
0323 }
0324 #ifdef EDM_ML_DEBUG
0325 edm::LogVerbatim("HCalGeom") << "getLayerFront::Input " << idet << ":" << ieta << ":" << iphi << ":" << depth
0326 << " Output " << layFront;
0327 #endif
0328 return layFront;
0329 }
0330
0331 int HcalDDDRecConstants::getMaxDepth(const int& itype, const int& ieta, const int& iphi, const int& zside) const {
0332 unsigned int type = (itype == 0) ? 0 : 1;
0333 int lmax = hcons.getMaxDepth(type + 1, ieta, iphi, zside, true);
0334 if (lmax < 0) {
0335 unsigned int lymax = (type == 0) ? maxLayerHB_ + 1 : maxLayer_ + 1;
0336 lmax = 0;
0337 if (layerGroupSize(ieta - 1) > 0) {
0338 if (layerGroupSize(ieta - 1) < lymax)
0339 lymax = layerGroupSize(ieta - 1);
0340 lmax = (int)(layerGroup(ieta - 1, lymax - 1));
0341 if (type == 0 && ieta == iEtaMax[type])
0342 lmax = hcons.getDepthEta16M(1);
0343 if (type == 1 && ieta >= hpar->noff[1])
0344 lmax = hcons.getDepthEta29M(0, false);
0345 }
0346 }
0347 #ifdef EDM_ML_DEBUG
0348 edm::LogVerbatim("HCalGeom") << "getMaxDepth::Input " << itype << ":" << ieta << ":" << iphi << ":" << zside
0349 << " Output " << lmax;
0350 #endif
0351 return lmax;
0352 }
0353
0354 int HcalDDDRecConstants::getMinDepth(const int& itype, const int& ieta, const int& iphi, const int& zside) const {
0355 int lmin = hcons.getMinDepth(itype + 1, ieta, iphi, zside, true);
0356 if (lmin < 0) {
0357 if (itype == 2) {
0358 lmin = 1;
0359 } else if (itype == 3) {
0360 lmin = maxDepth[3];
0361 } else {
0362 unsigned int type = (itype == 0) ? 0 : 1;
0363 if (layerGroupSize(ieta - 1) > 0) {
0364 if (type == 1 && ieta == iEtaMin[type])
0365 lmin = hcons.getDepthEta16M(2);
0366 else
0367 lmin = (int)(layerGroup(ieta - 1, 0));
0368 }
0369 }
0370 }
0371 return lmin;
0372 }
0373
0374 std::vector<std::pair<int, double>> HcalDDDRecConstants::getPhis(const int& subdet, const int& ieta) const {
0375 std::vector<std::pair<int, double>> phis;
0376 int ietaAbs = (ieta > 0) ? ieta : -ieta;
0377 int keta = (subdet != HcalForward) ? etaSimValu[ietaAbs - 1].first : ietaAbs;
0378 std::pair<double, double> ficons = hcons.getPhiCons(subdet, keta);
0379 double fioff = ficons.first;
0380 double dphi = (subdet != HcalForward) ? phibin[ietaAbs - 1] : ficons.second;
0381 int nphi = int((2._pi + 0.1 * dphi) / dphi);
0382 int units = hcons.unitPhi(subdet, keta);
0383 for (int ifi = 0; ifi < nphi; ++ifi) {
0384 double phi = -fioff + (ifi + 0.5) * dphi;
0385 int iphi = hcons.phiNumber(ifi + 1, units);
0386 phis.emplace_back(std::pair<int, double>(iphi, phi));
0387 }
0388 #ifdef EDM_ML_DEBUG
0389 edm::LogVerbatim("HCalGeom") << "getEtaPhi: subdet|ieta|iphi " << subdet << "|" << ieta << " with " << phis.size()
0390 << " phi bins";
0391 for (unsigned int k = 0; k < phis.size(); ++k)
0392 edm::LogVerbatim("HCalGeom") << "[" << k << "] iphi " << phis[k].first << " phi "
0393 << convertRadToDeg(phis[k].second);
0394 #endif
0395 return phis;
0396 }
0397
0398 int HcalDDDRecConstants::getPhiZOne(std::vector<std::pair<int, int>>& phiz) const {
0399 phiz.clear();
0400 int subdet = hcons.ldMap()->getSubdet();
0401 if (subdet > 0) {
0402 std::vector<int> phis = hcons.ldMap()->getPhis();
0403 for (int k : phis) {
0404 int zside = (k > 0) ? 1 : -1;
0405 int phi = (k > 0) ? k : -k;
0406 phiz.emplace_back(std::pair<int, int>(phi, zside));
0407 }
0408 }
0409 #ifdef EDM_ML_DEBUG
0410 edm::LogVerbatim("HCalGeom") << "Special RBX for detector " << subdet << " with " << phiz.size() << " phi/z bins";
0411 for (unsigned int k = 0; k < phiz.size(); ++k)
0412 edm::LogVerbatim("HCalGeom") << " [" << k << "] " << phiz[k].first << ":" << phiz[k].second;
0413 #endif
0414 return subdet;
0415 }
0416
0417 double HcalDDDRecConstants::getRZ(const int& subdet, const int& ieta, const int& depth) const {
0418 return getRZ(subdet, ieta, 1, depth);
0419 }
0420
0421 double HcalDDDRecConstants::getRZ(const int& subdet, const int& ieta, const int& iphi, const int& depth) const {
0422 int layf = getLayerFront(subdet, ieta, iphi, depth);
0423 double rz =
0424 (layf < 0) ? 0.0 : ((subdet == static_cast<int>(HcalBarrel)) ? (gconsHB[layf].first) : (gconsHE[layf].first));
0425 #ifdef EDM_ML_DEBUG
0426 edm::LogVerbatim("HCalGeom") << "getRZ: subdet|ieta|ipho|depth " << subdet << "|" << ieta << "|" << iphi << "|"
0427 << depth << " lay|rz " << layf << "|" << rz;
0428 #endif
0429 return rz;
0430 }
0431
0432 double HcalDDDRecConstants::getRZ(const int& subdet, const int& layer) const {
0433 double rz(0);
0434 if (layer > 0 && layer <= (int)(layerGroupSize(0)))
0435 rz = ((subdet == static_cast<int>(HcalBarrel)) ? (gconsHB[layer - 1].first) : (gconsHE[layer - 1].first));
0436 #ifdef EDM_ML_DEBUG
0437 edm::LogVerbatim("HCalGeom") << "getRZ: subdet|layer " << subdet << "|" << layer << " rz " << rz;
0438 #endif
0439 return rz;
0440 }
0441
0442 std::pair<double, double> HcalDDDRecConstants::getRZ(const HcalDetId& id) const {
0443 int subdet = id.subdetId();
0444 int ieta = id.ieta();
0445 int iphi = id.iphi();
0446 int depth = id.depth();
0447 int zside = (subdet == static_cast<int>(HcalBarrel)) ? 1 : id.zside();
0448 int layf = getLayerFront(subdet, ieta, iphi, depth);
0449 double rzf = (layf < 0)
0450 ? 0.0
0451 : ((subdet == static_cast<int>(HcalBarrel)) ? zside * (gconsHB[layf].first - gconsHB[layf].second)
0452 : zside * (gconsHE[layf].first - gconsHE[layf].second));
0453 int layb = getLayerBack(subdet, ieta, iphi, depth);
0454 double rzb = (layb < 0)
0455 ? 0.0
0456 : ((subdet == static_cast<int>(HcalBarrel)) ? zside * (gconsHB[layb].first + gconsHB[layb].second)
0457 : zside * (gconsHE[layb].first + gconsHE[layb].second));
0458 #ifdef EDM_ML_DEBUG
0459 edm::LogVerbatim("HCalGeom") << "getRZ: subdet|ieta|ipho|depth " << subdet << "|" << ieta << "|" << iphi << "|"
0460 << depth << " lay|rz (front) " << layf << "|" << rzf << " lay|rz (back) " << layb << "|"
0461 << rzb;
0462 #endif
0463 return std::pair<double, double>(rzf, rzb);
0464 }
0465
0466 std::vector<HcalDDDRecConstants::HcalActiveLength> HcalDDDRecConstants::getThickActive(const int& type) const {
0467 std::vector<HcalDDDRecConstants::HcalActiveLength> actives;
0468 std::vector<HcalDDDRecConstants::HcalEtaBin> bins = getEtaBins(type);
0469 #ifdef EDM_ML_DEBUG
0470 unsigned int kount(0);
0471 #endif
0472 for (auto& bin : bins) {
0473 int ieta = bin.ieta;
0474 int zside = bin.zside;
0475 int stype = (bin.phis.size() > 4) ? 0 : 1;
0476 int layf = getLayerFront(type + 1, zside * ieta, bin.phis[0].first, bin.depthStart) + 1;
0477 int layl = hcons.getLastLayer(type + 1, zside * ieta) + 1;
0478 double eta = 0.5 * (bin.etaMin + bin.etaMax);
0479 double theta = 2 * atan(exp(-eta));
0480 double scale = 1.0 / ((type == 0) ? sin(theta) : cos(theta));
0481 int depth = bin.depthStart;
0482 #ifdef EDM_ML_DEBUG
0483 edm::LogVerbatim("HCalGeom") << "Eta " << ieta << " zside " << zside << " depth " << depth << " Layers " << layf
0484 << ":" << layl << ":" << bin.layer.size();
0485 for (auto ll : bin.layer)
0486 edm::LogVerbatim("HCalGeom") << "Layer " << ll.first << ":" << ll.second;
0487 for (auto phi : bin.phis)
0488 edm::LogVerbatim("HCalGeom") << "Phi " << phi.first << ":" << convertRadToDeg(phi.second);
0489 #endif
0490 for (unsigned int i = 0; i < bin.layer.size(); ++i) {
0491 double thick(0);
0492 int lmin = (type == 1 && ieta == iEtaMin[1]) ? layf : std::max(bin.layer[i].first, layf);
0493 int lmax = std::min(bin.layer[i].second, layl);
0494 for (int j = lmin; j <= lmax; ++j) {
0495 double t = ((type == 0) ? gconsHB[j - 1].second : gconsHE[j - 1].second);
0496 if ((type == 1) && (ieta <= 18))
0497 t = gconsHE[j].second;
0498 if (t > 0)
0499 thick += t;
0500 }
0501 #ifdef EDM_ML_DEBUG
0502 edm::LogVerbatim("HCalGeom") << "Type " << type << " L " << lmin << ":" << lmax << " T " << thick;
0503 #endif
0504 thick *= (2. * scale);
0505 HcalDDDRecConstants::HcalActiveLength active(ieta, depth, zside, stype, zside * eta, thick);
0506 for (auto phi : bin.phis)
0507 active.iphis.emplace_back(phi.first);
0508 actives.emplace_back(active);
0509 ++depth;
0510 #ifdef EDM_ML_DEBUG
0511 kount++;
0512 edm::LogVerbatim("HCalGeom") << "getThickActive: [" << kount << "] eta:" << active.ieta << ":" << active.eta
0513 << " zside " << active.zside << " depth " << active.depth << " type " << active.stype
0514 << " thick " << active.thick;
0515 #endif
0516 }
0517 }
0518 return actives;
0519 }
0520
0521 std::vector<HcalCellType> HcalDDDRecConstants::HcalCellTypes(HcalSubdetector subdet) const {
0522 if (subdet == HcalBarrel || subdet == HcalEndcap) {
0523 std::vector<HcalCellType> cells;
0524 int isub = (subdet == HcalBarrel) ? 0 : 1;
0525 std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = getEtaBins(isub);
0526 std::vector<int> missPhi;
0527 for (const auto& etabin : etabins) {
0528 std::vector<HcalCellType> temp;
0529 std::vector<int> count;
0530 std::vector<double> dmin, dmax;
0531 for (unsigned int il = 0; il < etabin.layer.size(); ++il) {
0532 HcalCellType cell(subdet, etabin.ieta, etabin.zside, 0, HcalCellType::HcalCell());
0533 temp.emplace_back(cell);
0534 count.emplace_back(0);
0535 dmin.emplace_back(0);
0536 dmax.emplace_back(0);
0537 }
0538 int ieta = etabin.ieta;
0539 for (int keta = etaSimValu[ieta - 1].first; keta <= etaSimValu[ieta - 1].second; ++keta) {
0540 std::vector<HcalCellType> cellsm = hcons.HcalCellTypes(subdet, keta, -1);
0541 for (unsigned int il = 0; il < etabin.layer.size(); ++il) {
0542 for (auto& ic : cellsm) {
0543 if (ic.depthSegment() >= etabin.layer[il].first && ic.depthSegment() <= etabin.layer[il].second &&
0544 ic.etaBin() == temp[il].etaBin() && ic.zside() == temp[il].zside()) {
0545 if (count[il] == 0) {
0546 temp[il] = ic;
0547 dmin[il] = ic.depthMin();
0548 dmax[il] = ic.depthMax();
0549 }
0550 ++count[il];
0551 if (ic.depthMin() < dmin[il])
0552 dmin[il] = ic.depthMin();
0553 if (ic.depthMax() > dmax[il])
0554 dmax[il] = ic.depthMax();
0555 }
0556 }
0557 }
0558 }
0559 for (unsigned int il = 0; il < etabin.layer.size(); ++il) {
0560 int depth = etabin.depthStart + (int)(il);
0561 temp[il].setEta(ieta, etabin.etaMin, etabin.etaMax);
0562 temp[il].setDepth(depth, dmin[il], dmax[il]);
0563 double foff = (etabin.ieta <= iEtaMax[0]) ? hpar->phioff[0] : hpar->phioff[1];
0564 int unit = hcons.unitPhi(etabin.dphi);
0565 temp[il].setPhi(etabin.phis, missPhi, foff, etabin.dphi, unit);
0566 cells.emplace_back(temp[il]);
0567 }
0568 }
0569 #ifdef EDM_ML_DEBUG
0570 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants: found " << cells.size() << " cells for sub-detector type "
0571 << isub;
0572 for (unsigned int ic = 0; ic < cells.size(); ++ic)
0573 edm::LogVerbatim("HCalGeom") << "Cell[" << ic << "] " << cells[ic];
0574 #endif
0575 return cells;
0576 } else {
0577 return hcons.HcalCellTypes(subdet, -1, -1);
0578 }
0579 }
0580
0581 bool HcalDDDRecConstants::mergedDepthList29(int ieta, int iphi, int depth) const {
0582 int eta = std::abs(ieta);
0583 int zside = (ieta > 0) ? 1 : -1;
0584 int etamin = iEtaMax[1] - hcons.getDepthEta29(iphi, zside, 1);
0585 if ((eta >= etamin) && (eta <= iEtaMax[1])) {
0586 int depthMax = getMaxDepth(1, etamin, iphi, zside);
0587 int depthMin = hcons.getDepthEta29(iphi, zside, 0) + 1;
0588 if (depth >= depthMin && depth <= depthMax)
0589 return true;
0590 }
0591 return false;
0592 }
0593
0594 std::vector<int> HcalDDDRecConstants::mergedDepthList29(int ieta, int iphi) const {
0595 std::vector<int> depths;
0596 int eta = std::abs(ieta);
0597 int zside = (ieta > 0) ? 1 : -1;
0598 int etamin = iEtaMax[1] - hcons.getDepthEta29(iphi, zside, 1);
0599 if ((eta >= etamin) && (eta <= iEtaMax[1])) {
0600 int depthMax = getMaxDepth(1, etamin, iphi, zside);
0601 int depthMin = hcons.getDepthEta29(iphi, zside, 0) + 1;
0602 depths.reserve(depthMax - depthMin + 1);
0603 for (int depth = depthMin; depth <= depthMax; ++depth)
0604 depths.emplace_back(depth);
0605 }
0606 return depths;
0607 }
0608
0609 unsigned int HcalDDDRecConstants::numberOfCells(HcalSubdetector subdet) const {
0610 if (subdet == HcalBarrel || subdet == HcalEndcap) {
0611 unsigned int num = 0;
0612 std::vector<HcalCellType> cellTypes = HcalCellTypes(subdet);
0613 for (auto& cellType : cellTypes) {
0614 num += (unsigned int)(cellType.nPhiBins());
0615 }
0616 #ifdef EDM_ML_DEBUG
0617 edm::LogInfo("HCalGeom") << "HcalDDDRecConstants:numberOfCells " << cellTypes.size() << " " << num
0618 << " for subdetector " << subdet;
0619 #endif
0620 return num;
0621 } else {
0622 return hcons.numberOfCells(subdet);
0623 }
0624 }
0625
0626 unsigned int HcalDDDRecConstants::nCells(HcalSubdetector subdet) const {
0627 if (subdet == HcalBarrel || subdet == HcalEndcap) {
0628 int isub = (subdet == HcalBarrel) ? 0 : 1;
0629 std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = getEtaBins(isub);
0630 unsigned int ncell(0);
0631 for (auto& etabin : etabins) {
0632 ncell += ((etabin.phis.size()) * (etabin.layer.size()));
0633 }
0634 return ncell;
0635 } else if (subdet == HcalOuter) {
0636 return kHOSizePreLS1;
0637 } else if (subdet == HcalForward) {
0638 return (unsigned int)(hcons.numberOfCells(subdet));
0639 } else {
0640 return 0;
0641 }
0642 }
0643
0644 unsigned int HcalDDDRecConstants::nCells() const {
0645 return (nCells(HcalBarrel) + nCells(HcalEndcap) + nCells(HcalOuter) + nCells(HcalForward));
0646 }
0647
0648 HcalDetId HcalDDDRecConstants::mergedDepthDetId(const HcalDetId& id) const {
0649 std::map<HcalDetId, HcalDetId>::const_iterator itr = detIdSp_.find(id);
0650 if (itr == detIdSp_.end())
0651 return id;
0652 else
0653 return itr->second;
0654 }
0655
0656 HcalDetId HcalDDDRecConstants::idFront(const HcalDetId& id) const {
0657 HcalDetId hid(id);
0658 std::map<HcalDetId, std::vector<HcalDetId>>::const_iterator itr = detIdSpR_.find(id);
0659 if (itr != detIdSpR_.end())
0660 hid = HcalDetId(id.subdet(), id.ieta(), id.iphi(), (itr->second)[0].depth());
0661 return hid;
0662 }
0663
0664 HcalDetId HcalDDDRecConstants::idBack(const HcalDetId& id) const {
0665 HcalDetId hid(id);
0666 std::map<HcalDetId, std::vector<HcalDetId>>::const_iterator itr = detIdSpR_.find(id);
0667 if (itr != detIdSpR_.end())
0668 hid = HcalDetId(id.subdet(), id.ieta(), id.iphi(), (itr->second).back().depth());
0669 return hid;
0670 }
0671
0672 void HcalDDDRecConstants::unmergeDepthDetId(const HcalDetId& id, std::vector<HcalDetId>& ids) const {
0673 ids.clear();
0674 std::map<HcalDetId, std::vector<HcalDetId>>::const_iterator itr = detIdSpR_.find(id);
0675 if (itr == detIdSpR_.end()) {
0676 ids.emplace_back(id);
0677 } else {
0678 for (auto k : itr->second) {
0679 HcalDetId hid(id.subdet(), id.ieta(), id.iphi(), k.depth());
0680 ids.emplace_back(hid);
0681 }
0682 }
0683 }
0684
0685 void HcalDDDRecConstants::specialRBXHBHE(const std::vector<HcalDetId>& idsOld, std::vector<HcalDetId>& idsNew) const {
0686 for (auto k : idsOld) {
0687 std::map<HcalDetId, HcalDetId>::const_iterator itr = detIdSp_.find(k);
0688 if (itr == detIdSp_.end())
0689 idsNew.emplace_back(k);
0690 else
0691 idsNew.emplace_back(itr->second);
0692 }
0693 }
0694
0695 bool HcalDDDRecConstants::specialRBXHBHE(bool tobemerged, std::vector<HcalDetId>& ids) const {
0696 if (tobemerged) {
0697 std::map<HcalDetId, HcalDetId>::const_iterator itr;
0698 for (itr = detIdSp_.begin(); itr != detIdSp_.end(); ++itr)
0699 ids.emplace_back(itr->first);
0700 } else {
0701 std::map<HcalDetId, std::vector<HcalDetId>>::const_iterator itr;
0702 for (itr = detIdSpR_.begin(); itr != detIdSpR_.end(); ++itr)
0703 ids.emplace_back(itr->first);
0704 }
0705 return (!ids.empty());
0706 }
0707
0708 void HcalDDDRecConstants::getOneEtaBin(HcalSubdetector subdet,
0709 int ieta,
0710 int zside,
0711 std::vector<std::pair<int, double>>& phis,
0712 std::map<int, int>& layers,
0713 bool planOne,
0714 std::vector<HcalDDDRecConstants::HcalEtaBin>& bins) const {
0715 unsigned int lymax = (subdet == HcalBarrel) ? maxLayerHB_ + 1 : maxLayer_ + 1;
0716 int type = (subdet == HcalBarrel) ? 0 : 1;
0717 double dphi = phibin[ieta - 1];
0718 HcalDDDRecConstants::HcalEtaBin etabin =
0719 HcalDDDRecConstants::HcalEtaBin(ieta, zside, dphi, etaTable[ieta - 1], etaTable[ieta]);
0720 etabin.phis.insert(etabin.phis.end(), phis.begin(), phis.end());
0721 int n = (ieta == iEtaMax[type]) ? 0 : 1;
0722 HcalDDDRecConstants::HcalEtaBin etabin0 =
0723 HcalDDDRecConstants::HcalEtaBin(ieta, zside, dphi, etaTable[ieta - 1], etaTable[ieta + n]);
0724 etabin0.depthStart = hcons.getDepthEta29(phis[0].first, zside, 0) + 1;
0725 int dstart = -1;
0726 int lmin(0), lmax(0);
0727
0728 std::map<int, int>::iterator itr = layers.begin();
0729 if (!layers.empty()) {
0730 int dep = itr->second;
0731 if (subdet == HcalEndcap && ieta == iEtaMin[type])
0732 dep = hcons.getDepthEta16(subdet, phis[0].first, zside);
0733 unsigned lymx0 = (layers.size() > lymax) ? lymax : layers.size();
0734 #ifdef EDM_ML_DEBUG
0735 edm::LogVerbatim("HCalGeom") << "Eta " << ieta << ":" << hpar->noff[1] << " zside " << zside << " lymax " << lymx0
0736 << ":" << lymax << " Depth " << dep << ":" << itr->second;
0737 unsigned int l(0);
0738 for (itr = layers.begin(); itr != layers.end(); ++itr, ++l)
0739 edm::LogVerbatim("HCalGeom") << "Layer [" << l << "] " << itr->first << ":" << itr->second;
0740 edm::LogVerbatim("HCalGeom") << "With " << phis.size() << " phis";
0741 for (unsigned int l = 0; l < phis.size(); ++l)
0742 edm::LogVerbatim("HCalGeom") << "[" << l << "] " << phis[l].first << ":" << convertRadToDeg(phis[l].second);
0743 #endif
0744 for (itr = layers.begin(); itr != layers.end(); ++itr) {
0745 if (itr->first <= (int)(lymx0)) {
0746 if (itr->second == dep) {
0747 if (lmin == 0)
0748 lmin = itr->first;
0749 lmax = itr->first;
0750 } else if (itr->second > dep) {
0751 if (dstart < 0)
0752 dstart = dep;
0753 int lmax0 = (lmax >= lmin) ? lmax : lmin;
0754 if (subdet == HcalEndcap && ieta + 1 == hpar->noff[1] && dep > hcons.getDepthEta29(phis[0].first, zside, 0)) {
0755 etabin0.layer.emplace_back(std::pair<int, int>(lmin, lmax0));
0756 } else {
0757 etabin.layer.emplace_back(std::pair<int, int>(lmin, lmax0));
0758 }
0759 lmin = itr->first;
0760 lmax = lmin - 1;
0761 dep = itr->second;
0762 }
0763 if (subdet == HcalBarrel && ieta == iEtaMax[type] && dep > hcons.getDepthEta16M(1))
0764 break;
0765 if (subdet == HcalEndcap && ieta == hpar->noff[1] && dep > hcons.getDepthEta29M(0, planOne)) {
0766 lmax = lymx0;
0767 break;
0768 }
0769 if (itr->first == (int)(lymx0))
0770 lmax = lymx0;
0771 }
0772 }
0773 if (lmax >= lmin) {
0774 if (ieta + 1 == hpar->noff[1]) {
0775 etabin0.layer.emplace_back(std::pair<int, int>(lmin, lmax));
0776 etabin0.phis.insert(etabin0.phis.end(), phis.begin(), phis.end());
0777 bins.emplace_back(etabin0);
0778 #ifdef EDM_ML_DEBUG
0779 edm::LogVerbatim("HCalGeom") << "etabin0: dStatrt " << etabin0.depthStart << " layers " << etabin0.layer.size()
0780 << ":" << lmin << ":" << lmax << " phis " << phis.size();
0781 for (unsigned int k = 0; k < etabin0.layer.size(); ++k)
0782 edm::LogVerbatim("HCalGeom") << " [" << k << "] " << etabin0.layer[k].first << ":" << etabin0.layer[k].second;
0783 #endif
0784 } else if (ieta == hpar->noff[1]) {
0785 } else {
0786 etabin.layer.emplace_back(std::pair<int, int>(lmin, lmax));
0787 if (dstart < 0)
0788 dstart = dep;
0789 }
0790 }
0791 }
0792 etabin.depthStart = dstart;
0793 bins.emplace_back(etabin);
0794 #ifdef EDM_ML_DEBUG
0795 edm::LogVerbatim("HCalGeom") << "etabin: dStatrt " << etabin.depthStart << " layers " << etabin.layer.size() << ":"
0796 << lmin << ":" << lmax << " phis " << etabin.phis.size();
0797 for (unsigned int k = 0; k < etabin.layer.size(); ++k)
0798 edm::LogVerbatim("HCalGeom") << "[" << k << "] " << etabin.layer[k].first << ":" << etabin.layer[k].second;
0799 #endif
0800 }
0801
0802 void HcalDDDRecConstants::initialize(void) {
0803
0804 int nEta = (int)(hpar->etagroup.size());
0805 if (nEta != (int)(hpar->phigroup.size())) {
0806 edm::LogError("HCalGeom") << "HcalDDDRecConstants: sizes of the vectors "
0807 << " etaGroup (" << nEta << ") and phiGroup (" << hpar->phigroup.size()
0808 << ") do not match";
0809 throw cms::Exception("DDException") << "HcalDDDRecConstants: inconsistent array sizes" << nEta << ":"
0810 << hpar->phigroup.size();
0811 }
0812
0813
0814 iEtaMin = hpar->etaMin;
0815 iEtaMax = hpar->etaMax;
0816 etaTable.clear();
0817 ietaMap.clear();
0818 etaSimValu.clear();
0819 int ieta(0), ietaHB(0), ietaHE(0), ietaHEM(0);
0820 etaTable.emplace_back(hpar->etaTable[ieta]);
0821 for (int i = 0; i < nEta; ++i) {
0822 int ef = ieta + 1;
0823 ieta += (hpar->etagroup[i]);
0824 if (ieta >= (int)(hpar->etaTable.size())) {
0825 edm::LogError("HCalGeom") << "HcalDDDRecConstants: Going beyond the array boundary " << hpar->etaTable.size()
0826 << " at index " << i << " of etaTable from SimConstant";
0827 throw cms::Exception("DDException")
0828 << "HcalDDDRecConstants: Going beyond the array boundary " << hpar->etaTable.size() << " at index " << i
0829 << " of etaTable from SimConstant";
0830 } else {
0831 etaTable.emplace_back(hpar->etaTable[ieta]);
0832 etaSimValu.emplace_back(std::pair<int, int>(ef, ieta));
0833 }
0834 for (int k = 0; k < (hpar->etagroup[i]); ++k)
0835 ietaMap.emplace_back(i + 1);
0836 if (ieta <= hpar->etaMax[0])
0837 ietaHB = i + 1;
0838 if (ieta <= hpar->etaMin[1])
0839 ietaHE = i + 1;
0840 if (ieta <= hpar->etaMax[1])
0841 ietaHEM = i + 1;
0842 }
0843 iEtaMin[1] = ietaHE;
0844 iEtaMax[0] = ietaHB;
0845 iEtaMax[1] = ietaHEM;
0846
0847
0848 nPhiBins.clear();
0849 for (unsigned int k = 0; k < 4; ++k)
0850 nPhiBins.emplace_back(0);
0851 ieta = 0;
0852 phibin.clear();
0853 phiUnitS.clear();
0854 for (int i = 0; i < nEta; ++i) {
0855 double dphi = (hpar->phigroup[i]) * (hpar->phibin[ieta]);
0856 phibin.emplace_back(dphi);
0857 int nphi = (int)((2._pi + 0.001) / dphi);
0858 if (ieta <= iEtaMax[0]) {
0859 if (nphi > nPhiBins[0])
0860 nPhiBins[3] = nPhiBins[0] = nphi;
0861 }
0862 if (ieta >= iEtaMin[1]) {
0863 if (nphi > nPhiBins[1])
0864 nPhiBins[1] = nphi;
0865 }
0866 ieta += (hpar->etagroup[i]);
0867 }
0868 for (unsigned int i = 1; i < hpar->etaTable.size(); ++i) {
0869 int unit = hcons.unitPhi(hpar->phibin[i - 1]);
0870 phiUnitS.emplace_back(unit);
0871 }
0872 for (double i : hpar->phitable) {
0873 int nphi = (int)((2._pi + 0.001) / i);
0874 if (nphi > nPhiBins[2])
0875 nPhiBins[2] = nphi;
0876 }
0877 #ifdef EDM_ML_DEBUG
0878 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants: Modified eta/deltaphi table for " << nEta << " bins";
0879 for (int i = 0; i < nEta; ++i)
0880 edm::LogVerbatim("HCalGeom") << "Eta[" << i << "] = " << etaTable[i] << ":" << etaTable[i + 1] << ":"
0881 << etaSimValu[i].first << ":" << etaSimValu[i].second << " PhiBin[" << i
0882 << "] = " << convertRadToDeg(phibin[i]);
0883 for (unsigned int i = 0; i < phiUnitS.size(); ++i)
0884 edm::LogVerbatim("HCalGeom") << " PhiUnitS[" << i << "] = " << phiUnitS[i];
0885 for (unsigned int i = 0; i < nPhiBins.size(); ++i)
0886 edm::LogVerbatim("HCalGeom") << " nPhiBins[" << i << "] = " << nPhiBins[i];
0887 for (unsigned int i = 0; i < hpar->etaTableHF.size(); ++i)
0888 edm::LogVerbatim("HCalGeom") << " EtaTableHF[" << i << "] = " << hpar->etaTableHF[i];
0889 for (unsigned int i = 0; i < hpar->phitable.size(); ++i)
0890 edm::LogVerbatim("HCalGeom") << " PhiBinHF[" << i << "] = " << hpar->phitable[i];
0891 #endif
0892
0893
0894 maxDepth = hpar->maxDepth;
0895 maxDepth[0] = maxDepth[1] = 0;
0896 for (int i = 0; i < nEta; ++i) {
0897 unsigned int imx = layerGroupSize(i);
0898 int laymax = (imx > 0) ? layerGroup(i, imx - 1) : 0;
0899 if (i < iEtaMax[0]) {
0900 int laymax0 = (imx > 16) ? layerGroup(i, 16) : laymax;
0901 if (i + 1 == iEtaMax[0])
0902 laymax0 = hcons.getDepthEta16M(1);
0903 #ifdef EDM_ML_DEBUG
0904 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:HB " << i << " " << imx << " " << laymax << " " << laymax0;
0905 #endif
0906 if (maxDepth[0] < laymax0)
0907 maxDepth[0] = laymax0;
0908 }
0909 if (i >= iEtaMin[1] - 1 && i < iEtaMax[1]) {
0910 #ifdef EDM_ML_DEBUG
0911 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:HE " << i << " " << imx << " " << laymax;
0912 #endif
0913 if (maxDepth[1] < laymax)
0914 maxDepth[1] = laymax;
0915 }
0916 }
0917 #ifdef EDM_ML_DEBUG
0918 for (int i = 0; i < 4; ++i)
0919 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:Detector Type[" << i << "] iEta " << iEtaMin[i] << ":"
0920 << iEtaMax[i] << " MaxDepth " << maxDepth[i];
0921 #endif
0922
0923
0924 nModule[0] = hpar->modHB[0];
0925 nHalves[0] = hpar->modHB[1];
0926 for (unsigned int i = 0; i < hpar->rHB.size(); ++i) {
0927 gconsHB.emplace_back(std::pair<double, double>(hpar->rHB[i] / CLHEP::cm, hpar->drHB[i] / CLHEP::cm));
0928 }
0929 #ifdef EDM_ML_DEBUG
0930 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:HB with " << nModule[0] << " modules and " << nHalves[0]
0931 << " halves and " << gconsHB.size() << " layers";
0932 for (unsigned int i = 0; i < gconsHB.size(); ++i)
0933 edm::LogVerbatim("HCalGeom") << "rHB[" << i << "] = " << gconsHB[i].first << " +- " << gconsHB[i].second;
0934 #endif
0935 nModule[1] = hpar->modHE[0];
0936 nHalves[1] = hpar->modHE[1];
0937 for (unsigned int i = 0; i < hpar->zHE.size(); ++i) {
0938 gconsHE.emplace_back(std::pair<double, double>(hpar->zHE[i] / CLHEP::cm, hpar->dzHE[i] / CLHEP::cm));
0939 }
0940 #ifdef EDM_ML_DEBUG
0941 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:HE with " << nModule[1] << " modules and " << nHalves[1]
0942 << " halves and " << gconsHE.size() << " layers";
0943 for (unsigned int i = 0; i < gconsHE.size(); ++i)
0944 edm::LogVerbatim("HCalGeom") << "zHE[" << i << "] = " << gconsHE[i].first << " +- " << gconsHE[i].second;
0945 #endif
0946
0947
0948 depthMaxSp_ = hcons.getMaxDepthDet(0);
0949 if (depthMaxSp_.first == 0) {
0950 depthMaxSp_ = depthMaxDf_ = std::pair<int, int>(2, maxDepth[1]);
0951 } else if (depthMaxSp_.first == 1) {
0952 depthMaxDf_ = std::pair<int, int>(1, maxDepth[0]);
0953 if (depthMaxSp_.second > maxDepth[0])
0954 maxDepth[0] = depthMaxSp_.second;
0955 } else {
0956 depthMaxDf_ = std::pair<int, int>(2, maxDepth[1]);
0957 if (depthMaxSp_.second > maxDepth[1])
0958 maxDepth[1] = depthMaxSp_.second;
0959 }
0960 #ifdef EDM_ML_DEBUG
0961 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:Detector type and maximum depth for all RBX "
0962 << depthMaxDf_.first << ":" << depthMaxDf_.second << " and for special RBX "
0963 << depthMaxSp_.first << ":" << depthMaxSp_.second;
0964 #endif
0965
0966
0967 std::vector<int> phis;
0968 HcalSubdetector subdet = HcalSubdetector(hcons.ldMap()->validDet(phis));
0969 detIdSp_.clear();
0970 detIdSpR_.clear();
0971 if ((subdet == HcalBarrel) || (subdet == HcalEndcap)) {
0972 int phi = (phis[0] > 0) ? phis[0] : -phis[0];
0973 int zside = (phis[0] > 0) ? 1 : -1;
0974 int lymax = (subdet == HcalBarrel) ? maxLayerHB_ + 1 : maxLayer_ + 1;
0975 std::pair<int, int> etas = hcons.ldMap()->validEta();
0976 for (int eta = etas.first; eta <= etas.second; ++eta) {
0977 std::map<int, std::pair<int, int>> oldDep;
0978 int depth(0);
0979 int lmin = layerGroup(eta - 1, 0);
0980 for (int lay = 0; lay < lymax; ++lay) {
0981 int depc = layerGroup(eta - 1, lay);
0982 if (depth != depc) {
0983 if (depth != 0)
0984 oldDep[depth] = std::pair<int, int>(lmin, lay - 1);
0985 depth = depc;
0986 lmin = lay;
0987 }
0988 }
0989 if (depth != 0)
0990 oldDep[depth] = std::pair<int, int>(lmin, lymax - 1);
0991 #ifdef EDM_ML_DEBUG
0992 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:Eta|Phi|Zside " << eta << ":" << phi << ":" << zside
0993 << " with " << oldDep.size() << " old Depths";
0994 unsigned int kk(0);
0995 for (std::map<int, std::pair<int, int>>::const_iterator itr = oldDep.begin(); itr != oldDep.end(); ++itr, ++kk)
0996 edm::LogVerbatim("HCalGeom") << "[" << kk << "] " << itr->first << " --> " << itr->second.first << ":"
0997 << itr->second.second;
0998 #endif
0999 std::pair<int, int> depths = hcons.ldMap()->getDepths(eta);
1000 for (int ndepth = depths.first; ndepth <= depths.second; ++ndepth) {
1001 bool flag = ((subdet == HcalBarrel && eta == iEtaMax[0] && ndepth > hcons.getDepthEta16(subdet, phi, zside)) ||
1002 (subdet == HcalEndcap && eta == iEtaMin[1] && ndepth < hcons.getDepthEta16(subdet, phi, zside)));
1003 if (!flag) {
1004 std::vector<int> count(oldDep.size(), 0);
1005 int layFront = hcons.ldMap()->getLayerFront(subdet, eta, phi, zside, ndepth);
1006 int layBack = hcons.ldMap()->getLayerBack(subdet, eta, phi, zside, ndepth);
1007 for (int lay = layFront; lay <= layBack; ++lay) {
1008 unsigned int l(0);
1009 for (std::map<int, std::pair<int, int>>::iterator itr = oldDep.begin(); itr != oldDep.end(); ++itr, ++l) {
1010 if (lay >= (itr->second).first && lay <= (itr->second).second) {
1011 ++count[l];
1012 break;
1013 }
1014 }
1015 }
1016 int odepth(0), maxlay(0);
1017 unsigned int l(0);
1018 for (std::map<int, std::pair<int, int>>::iterator itr = oldDep.begin(); itr != oldDep.end(); ++itr, ++l) {
1019 if (count[l] > maxlay) {
1020 odepth = itr->first;
1021 maxlay = count[l];
1022 }
1023 }
1024 #ifdef EDM_ML_DEBUG
1025 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:New Depth " << ndepth << " old Depth " << odepth
1026 << " max " << maxlay;
1027 #endif
1028 for (int k : phis) {
1029 zside = (k > 0) ? 1 : -1;
1030 phi = (k > 0) ? k : -k;
1031 if (subdet == HcalEndcap && eta == hpar->noff[1] && ndepth > hcons.getDepthEta29M(0, true))
1032 break;
1033 HcalDetId newId(subdet, zside * eta, phi, ndepth);
1034 HcalDetId oldId(subdet, zside * eta, phi, odepth);
1035 detIdSp_[newId] = oldId;
1036 std::vector<HcalDetId> ids;
1037 std::map<HcalDetId, std::vector<HcalDetId>>::iterator itr = detIdSpR_.find(oldId);
1038 if (itr != detIdSpR_.end())
1039 ids = itr->second;
1040 ids.emplace_back(newId);
1041 detIdSpR_[oldId] = ids;
1042 }
1043 }
1044 }
1045 }
1046 #ifdef EDM_ML_DEBUG
1047 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:Map for merging new channels to old channel"
1048 << " IDs with " << detIdSp_.size() << " entries";
1049 int l(0);
1050 for (auto itr : detIdSp_) {
1051 edm::LogVerbatim("HCalGeom") << "[" << l << "] Special " << itr.first << " Standard " << itr.second;
1052 ++l;
1053 }
1054 edm::LogVerbatim("HCalGeom") << "HcalDDDRecConstants:Reverse Map for mapping old to new IDs with "
1055 << detIdSpR_.size() << " entries";
1056 l = 0;
1057 for (auto itr : detIdSpR_) {
1058 edm::LogVerbatim("HCalGeom") << "[" << l << "] Standard " << itr.first << " Special";
1059 for (auto itr1 : itr.second)
1060 edm::LogVerbatim("HCalGeom") << "ID " << (itr1);
1061 ++l;
1062 }
1063 #endif
1064 }
1065 }
1066
1067 unsigned int HcalDDDRecConstants::layerGroupSize(int eta) const {
1068 unsigned int k = 0;
1069 for (auto const& it : hpar->layerGroupEtaRec) {
1070 if (it.layer == (unsigned int)(eta + 1)) {
1071 return it.layerGroup.size();
1072 }
1073 if (it.layer > (unsigned int)(eta + 1))
1074 break;
1075 k = it.layerGroup.size();
1076 }
1077 return k;
1078 }
1079
1080 unsigned int HcalDDDRecConstants::layerGroup(int eta, int i) const {
1081 unsigned int k = 0;
1082 for (auto const& it : hpar->layerGroupEtaRec) {
1083 if (it.layer == (unsigned int)(eta + 1)) {
1084 return it.layerGroup.at(i);
1085 }
1086 if (it.layer > (unsigned int)(eta + 1))
1087 break;
1088 k = it.layerGroup.at(i);
1089 }
1090 return k;
1091 }
