Back to home page

Project CMSSW displayed by LXR

 
 

    


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 //#define EDM_ML_DEBUG
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))) {  // Use Eta Table
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) {  // HFn
0363       lmin = 1;
0364     } else if (itype == 3) {  //HO
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   //Eta grouping
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   // First eta table
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   // Then Phi bins
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   //Now the depths
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   //Now the geometry constants
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   //Special RBX
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   //Map of special DetId's
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 }