Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:20:50 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!