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