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