Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-17 23:36:21

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