Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:49:33

0001 #include "Geometry/HcalCommonData/interface/HcalGeomParameters.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 
0006 #include "DataFormats/Math/interface/GeantUnits.h"
0007 #include "DetectorDescription/Core/interface/DDutils.h"
0008 #include "DetectorDescription/Core/interface/DDValue.h"
0009 #include "DetectorDescription/Core/interface/DDFilter.h"
0010 #include "DetectorDescription/Core/interface/DDSolid.h"
0011 #include "DetectorDescription/Core/interface/DDConstant.h"
0012 #include "DetectorDescription/RegressionTest/interface/DDErrorDetection.h"
0013 #include "CondFormats/GeometryObjects/interface/HcalParameters.h"
0014 
0015 //#define EDM_ML_DEBUG
0016 using namespace geant_units::operators;
0017 
0018 static const double tan10deg = std::tan(10._deg);
0019 
0020 void HcalGeomParameters::getConstRHO(std::vector<double>& rHO) const {
0021   rHO.emplace_back(rminHO_);
0022   for (double i : etaHO_)
0023     rHO.emplace_back(i);
0024 }
0025 
0026 std::vector<int> HcalGeomParameters::getModHalfHBHE(const int type) const {
0027   std::vector<int> modHalf;
0028   if (type == 0) {
0029     modHalf.emplace_back(nmodHB_);
0030     modHalf.emplace_back(nzHB_);
0031   } else {
0032     modHalf.emplace_back(nmodHE_);
0033     modHalf.emplace_back(nzHE_);
0034   }
0035   return modHalf;
0036 }
0037 
0038 void HcalGeomParameters::loadGeometry(const DDFilteredView& _fv, HcalParameters& php) {
0039   DDFilteredView fv = _fv;
0040   bool dodet = true;
0041   bool hf = false;
0042   clear(php);
0043 
0044   while (dodet) {
0045     DDTranslation t = fv.translation();
0046     std::vector<int> copy = fv.copyNumbers();
0047     const DDSolid& sol = fv.logicalPart().solid();
0048     int idet = 0, lay = -1;
0049     int nsiz = static_cast<int>(copy.size());
0050     if (nsiz > 0)
0051       lay = copy[nsiz - 1] / 10;
0052     if (nsiz > 6)
0053       idet = copy[nsiz - 2] / 1000;
0054     else if (nsiz > 4)
0055       idet = copy[nsiz - 3] / 1000;
0056     else if (nsiz > 1)
0057       idet = copy[nsiz - 2] / 1000;
0058 #ifdef EDM_ML_DEBUG
0059     std::ostringstream st1;
0060     for (unsigned int k = 0; k < copy.size(); ++k)
0061       st1 << " " << copy[k];
0062     edm::LogVerbatim("HCalGeom") << "Name " << fv.logicalPart().solid().name() << " Copy " << copy.size() << ":"
0063                                  << st1.str();
0064 #endif
0065     double dx = 0, dy = 0, dz = 0, dx1 = 0, dx2 = 0;
0066     double alp(0);
0067     if (sol.shape() == DDSolidShape::ddbox) {
0068       const DDBox& box = static_cast<DDBox>(sol);
0069       dx = HcalGeomParameters::k_ScaleFromDDDToG4 * box.halfX();
0070       dy = HcalGeomParameters::k_ScaleFromDDDToG4 * box.halfY();
0071       dz = HcalGeomParameters::k_ScaleFromDDDToG4 * box.halfZ();
0072     } else if (sol.shape() == DDSolidShape::ddtrap) {
0073       const DDTrap& trp = static_cast<DDTrap>(sol);
0074       dx1 = HcalGeomParameters::k_ScaleFromDDDToG4 * trp.x1();
0075       dx2 = HcalGeomParameters::k_ScaleFromDDDToG4 * trp.x2();
0076       dx = 0.25 * HcalGeomParameters::k_ScaleFromDDDToG4 * (trp.x1() + trp.x2() + trp.x3() + trp.x4());
0077       dy = 0.5 * HcalGeomParameters::k_ScaleFromDDDToG4 * (trp.y1() + trp.y2());
0078       dz = HcalGeomParameters::k_ScaleFromDDDToG4 * trp.halfZ();
0079       alp = 0.5 * (trp.alpha1() + trp.alpha2());
0080     } else if (sol.shape() == DDSolidShape::ddtubs) {
0081       const DDTubs& tub = static_cast<DDTubs>(sol);
0082       dx = HcalGeomParameters::k_ScaleFromDDDToG4 * tub.rIn();
0083       dy = HcalGeomParameters::k_ScaleFromDDDToG4 * tub.rOut();
0084       dz = HcalGeomParameters::k_ScaleFromDDDToG4 * tub.zhalf();
0085     }
0086     if (idet == 3) {
0087       // HB
0088 #ifdef EDM_ML_DEBUG
0089       edm::LogVerbatim("HCalGeom") << "HB " << sol.name() << " Shape " << sol.shape() << " Layer " << lay << " R "
0090                                    << t.Rho();
0091 #endif
0092       if (lay >= 0 && lay < maxLayer_) {
0093         ib_[lay]++;
0094         rb_[lay] += (HcalGeomParameters::k_ScaleFromDDDToG4 * t.Rho());
0095         if (thkb_[lay] <= 0) {
0096           if (lay < 17)
0097             thkb_[lay] = dx;
0098           else
0099             thkb_[lay] = std::min(dx, dy);
0100         }
0101         if (lay < 17) {
0102           bool found = false;
0103           for (double k : rxb_) {
0104             if (std::abs(k - (HcalGeomParameters::k_ScaleFromDDDToG4 * t.Rho())) < 0.01) {
0105               found = true;
0106               break;
0107             }
0108           }
0109           if (!found) {
0110             rxb_.emplace_back(HcalGeomParameters::k_ScaleFromDDDToG4 * t.Rho());
0111             php.rhoxHB.emplace_back(HcalGeomParameters::k_ScaleFromDDDToG4 * t.Rho() * std::cos(t.phi()));
0112             php.zxHB.emplace_back(HcalGeomParameters::k_ScaleFromDDDToG4 * std::abs(t.z()));
0113             php.dyHB.emplace_back(2. * dy);
0114             php.dxHB.emplace_back(2. * dz);
0115             php.layHB.emplace_back(lay);
0116           }
0117         }
0118       }
0119       if (lay == 2) {
0120         int iz = copy[nsiz - 5];
0121         int fi = copy[nsiz - 4];
0122         unsigned int it1 = find(iz, izb_);
0123         if (it1 == izb_.size())
0124           izb_.emplace_back(iz);
0125         unsigned int it2 = find(fi, phib_);
0126         if (it2 == phib_.size())
0127           phib_.emplace_back(fi);
0128       }
0129       if (lay == 18) {
0130         int ifi = -1, ich = -1;
0131         if (nsiz > 2)
0132           ifi = copy[nsiz - 3];
0133         if (nsiz > 3)
0134           ich = copy[nsiz - 4];
0135         double z1 = std::abs((HcalGeomParameters::k_ScaleFromDDDToG4 * t.z()) + dz);
0136         double z2 = std::abs((HcalGeomParameters::k_ScaleFromDDDToG4 * t.z()) - dz);
0137         if (std::abs(z1 - z2) < 0.01)
0138           z1 = 0;
0139         if (ifi == 1 && ich == 4) {
0140           if (z1 > z2) {
0141             double tmp = z1;
0142             z1 = z2;
0143             z2 = tmp;
0144           }
0145           bool sok = true;
0146           for (unsigned int kk = 0; kk < php.zHO.size(); kk++) {
0147             if (std::abs(z2 - php.zHO[kk]) < 0.01) {
0148               sok = false;
0149               break;
0150             } else if (z2 < php.zHO[kk]) {
0151               php.zHO.resize(php.zHO.size() + 2);
0152               for (unsigned int kz = php.zHO.size() - 1; kz > kk + 1; kz = kz - 2) {
0153                 php.zHO[kz] = php.zHO[kz - 2];
0154                 php.zHO[kz - 1] = php.zHO[kz - 3];
0155               }
0156               php.zHO[kk + 1] = z2;
0157               php.zHO[kk] = z1;
0158               sok = false;
0159               break;
0160             }
0161           }
0162           if (sok) {
0163             php.zHO.emplace_back(z1);
0164             php.zHO.emplace_back(z2);
0165           }
0166 #ifdef EDM_ML_DEBUG
0167           edm::LogVerbatim("HCalGeom") << "Detector " << idet << " Lay " << lay << " fi " << ifi << " " << ich << " z "
0168                                        << z1 << " " << z2;
0169 #endif
0170         }
0171       }
0172     } else if (idet == 4) {
0173       // HE
0174 #ifdef EDM_ML_DEBUG
0175       edm::LogVerbatim("HCalGeom") << "HE " << sol.name() << " Shape " << sol.shape() << " Layer " << lay << " Z "
0176                                    << t.z();
0177 #endif
0178       if (lay >= 0 && lay < maxLayer_) {
0179         ie_[lay]++;
0180         ze_[lay] += std::abs(HcalGeomParameters::k_ScaleFromDDDToG4 * t.z());
0181         if (thke_[lay] <= 0)
0182           thke_[lay] = dz;
0183         double rinHE = HcalGeomParameters::k_ScaleFromDDDToG4 * t.Rho() * cos(alp) - dy;
0184         double routHE = HcalGeomParameters::k_ScaleFromDDDToG4 * t.Rho() * cos(alp) + dy;
0185         rminHE_[lay] += rinHE;
0186         rmaxHE_[lay] += routHE;
0187         bool found = false;
0188         for (double k : php.zxHE) {
0189           if (std::abs(k - std::abs(HcalGeomParameters::k_ScaleFromDDDToG4 * t.z())) < 0.01) {
0190             found = true;
0191             break;
0192           }
0193         }
0194         if (!found) {
0195           php.zxHE.emplace_back(HcalGeomParameters::k_ScaleFromDDDToG4 * std::abs(t.z()));
0196           php.rhoxHE.emplace_back(HcalGeomParameters::k_ScaleFromDDDToG4 * t.Rho() * std::cos(t.phi()));
0197           php.dyHE.emplace_back(dy * std::cos(t.phi()));
0198           dx1 -= 0.5 * (HcalGeomParameters::k_ScaleFromDDDToG4 * t.rho() - dy) * std::cos(t.phi()) * tan10deg;
0199           dx2 -= 0.5 * (HcalGeomParameters::k_ScaleFromDDDToG4 * t.rho() + dy) * std::cos(t.phi()) * tan10deg;
0200           php.dx1HE.emplace_back(-dx1);
0201           php.dx2HE.emplace_back(-dx2);
0202           php.layHE.emplace_back(lay);
0203         }
0204       }
0205       if (copy[nsiz - 1] == kHELayer1_ || copy[nsiz - 1] == kHELayer2_) {
0206         int iz = copy[nsiz - 7];
0207         int fi = copy[nsiz - 5];
0208         unsigned int it1 = find(iz, ize_);
0209         if (it1 == ize_.size())
0210           ize_.emplace_back(iz);
0211         unsigned int it2 = find(fi, phie_);
0212         if (it2 == phie_.size())
0213           phie_.emplace_back(fi);
0214       }
0215     } else if (idet == 5) {
0216       // HF
0217       if (!hf) {
0218         const std::vector<double>& paras = sol.parameters();
0219 #ifdef EDM_ML_DEBUG
0220         edm::LogVerbatim("HCalGeom") << "HF " << sol.name() << " Shape " << sol.shape() << " Z " << t.z() << " with "
0221                                      << paras.size() << " Parameters";
0222         for (unsigned j = 0; j < paras.size(); j++)
0223           edm::LogVerbatim("HCalGeom") << "HF Parameter[" << j << "] = " << paras[j];
0224 #endif
0225         if (sol.shape() == DDSolidShape::ddpolycone_rrz) {
0226           int nz = (int)(paras.size()) - 3;
0227           dzVcal_ = 0.5 * HcalGeomParameters::k_ScaleFromDDDToG4 * (paras[nz] - paras[3]);
0228           hf = true;
0229         } else if (sol.shape() == DDSolidShape::ddtubs || sol.shape() == DDSolidShape::ddcons) {
0230           dzVcal_ = HcalGeomParameters::k_ScaleFromDDDToG4 * paras[0];
0231           hf = true;
0232         } else if (sol.shape() == DDSolidShape::ddtrap) {
0233           dzVcal_ = HcalGeomParameters::k_ScaleFromDDDToG4 * paras[0];
0234           hf = true;
0235         }
0236 #ifdef EDM_ML_DEBUG
0237         edm::LogVerbatim("HCalGeom") << "Sets for Detector " << idet << " for " << sol.name() << " Flag " << hf
0238                                      << " DZ " << dzVcal_;
0239 #endif
0240       }
0241 #ifdef EDM_ML_DEBUG
0242     } else {
0243       edm::LogVerbatim("HCalGeom") << "Unknown Detector " << idet << " for " << sol.name() << " Shape " << sol.shape()
0244                                    << " R " << t.Rho() << " Z " << t.z();
0245 #endif
0246     }
0247     dodet = fv.next();
0248   }
0249 
0250   loadfinal(php);
0251 }
0252 
0253 void HcalGeomParameters::loadGeometry(const cms::DDCompactView& cpv, HcalParameters& php) {
0254   const cms::DDFilter filter("OnlyForHcalSimNumbering", "HCAL");
0255   cms::DDFilteredView fv(cpv, filter);
0256   clear(php);
0257   bool hf(false);
0258   while (fv.firstChild()) {
0259     auto t = fv.translation();
0260     std::vector<double> paras = fv.parameters();
0261     std::vector<int> copy = fv.copyNos();
0262     int idet = 0, lay = -1;
0263     int nsiz = static_cast<int>(copy.size());
0264     if (nsiz > 0)
0265       lay = copy[0] / 10;
0266     if (nsiz > 6)
0267       idet = copy[1] / 1000;
0268     else if (nsiz > 4)
0269       idet = copy[2] / 1000;
0270     else if (nsiz > 1)
0271       idet = copy[1] / 1000;
0272 #ifdef EDM_ML_DEBUG
0273     std::ostringstream st1;
0274     for (unsigned int k = 0; k < copy.size(); ++k)
0275       st1 << " " << copy[k];
0276     edm::LogVerbatim("HCalGeom") << "Name " << fv.name() << " Copy " << copy.size() << ":" << st1.str();
0277     for (unsigned int n = 0; n < copy.size(); ++n)
0278       edm::LogVerbatim("HCalGeom") << "[" << n << "] " << copy[n];
0279     edm::LogVerbatim("HCalGeom") << "Detector " << idet << " Layer " << lay << " parameters: " << paras.size();
0280     for (unsigned int n = 0; n < paras.size(); ++n)
0281       edm::LogVerbatim("HCalGeom") << "[" << n << "] " << paras[n];
0282 #endif
0283     double dx = 0, dy = 0, dz = 0, dx1 = 0, dx2 = 0;
0284     double alp(0);
0285     if (dd4hep::isA<dd4hep::Box>(fv.solid())) {
0286       dx = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[0];
0287       dy = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[1];
0288       dz = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[2];
0289     } else if (dd4hep::isA<dd4hep::Trap>(fv.solid())) {
0290       dx1 = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[4];
0291       dx2 = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[5];
0292       dx = 0.25 * HcalGeomParameters::k_ScaleFromDD4hepToG4 * (paras[4] + paras[5] + paras[8] + paras[9]);
0293       dy = 0.5 * HcalGeomParameters::k_ScaleFromDD4hepToG4 * (paras[3] + paras[7]);
0294       dz = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[0];
0295       alp = 0.5 * (paras[6] + paras[10]);
0296     } else if (dd4hep::isA<dd4hep::Tube>(fv.solid())) {
0297       dx = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[0];
0298       dy = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[1];
0299       dz = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[2];
0300     }
0301     if (idet == 3) {
0302       // HB
0303 #ifdef EDM_ML_DEBUG
0304       edm::LogVerbatim("HCalGeom") << "HB " << fv.name() << " Shape " << cms::dd::name(cms::DDSolidShapeMap, fv.shape())
0305                                    << " Layer " << lay << " R " << t.Rho();
0306 #endif
0307       if (lay >= 0 && lay < maxLayer_) {
0308         ib_[lay]++;
0309         rb_[lay] += (HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho());
0310         if (thkb_[lay] <= 0) {
0311           if (lay < 17)
0312             thkb_[lay] = dx;
0313           else
0314             thkb_[lay] = std::min(dx, dy);
0315         }
0316         if (lay < 17) {
0317           bool found = false;
0318           for (double k : rxb_) {
0319             if (std::abs(k - (HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho())) < 0.01) {
0320               found = true;
0321               break;
0322             }
0323           }
0324           if (!found) {
0325             rxb_.emplace_back(HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho());
0326             php.rhoxHB.emplace_back(HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho() * std::cos(t.phi()));
0327             php.zxHB.emplace_back(std::abs(HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z()));
0328             php.dyHB.emplace_back(2. * dy);
0329             php.dxHB.emplace_back(2. * dz);
0330             php.layHB.emplace_back(lay);
0331           }
0332         }
0333       }
0334       if (lay == 2) {
0335         int iz = copy[4];
0336         int fi = copy[3];
0337         unsigned int it1 = find(iz, izb_);
0338         if (it1 == izb_.size())
0339           izb_.emplace_back(iz);
0340         unsigned int it2 = find(fi, phib_);
0341         if (it2 == phib_.size())
0342           phib_.emplace_back(fi);
0343       }
0344       if (lay == 18) {
0345         int ifi = -1, ich = -1;
0346         if (nsiz > 2)
0347           ifi = copy[2];
0348         if (nsiz > 3)
0349           ich = copy[3];
0350         double z1 = std::abs((HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z()) + dz);
0351         double z2 = std::abs((HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z()) - dz);
0352         if (std::abs(z1 - z2) < 0.01)
0353           z1 = 0;
0354         if (ifi == 1 && ich == 4) {
0355           if (z1 > z2) {
0356             double tmp = z1;
0357             z1 = z2;
0358             z2 = tmp;
0359           }
0360           bool sok = true;
0361           for (unsigned int kk = 0; kk < php.zHO.size(); kk++) {
0362             if (std::abs(z2 - php.zHO[kk]) < 0.01) {
0363               sok = false;
0364               break;
0365             } else if (z2 < php.zHO[kk]) {
0366               php.zHO.resize(php.zHO.size() + 2);
0367               for (unsigned int kz = php.zHO.size() - 1; kz > kk + 1; kz = kz - 2) {
0368                 php.zHO[kz] = php.zHO[kz - 2];
0369                 php.zHO[kz - 1] = php.zHO[kz - 3];
0370               }
0371               php.zHO[kk + 1] = z2;
0372               php.zHO[kk] = z1;
0373               sok = false;
0374               break;
0375             }
0376           }
0377           if (sok) {
0378             php.zHO.emplace_back(z1);
0379             php.zHO.emplace_back(z2);
0380           }
0381 #ifdef EDM_ML_DEBUG
0382           edm::LogVerbatim("HCalGeom") << "Detector " << idet << " Lay " << lay << " fi " << ifi << " " << ich << " z "
0383                                        << z1 << " " << z2;
0384 #endif
0385         }
0386       }
0387     } else if (idet == 4) {
0388       // HE
0389 #ifdef EDM_ML_DEBUG
0390       edm::LogVerbatim("HCalGeom") << "HE " << fv.name() << " Shape " << cms::dd::name(cms::DDSolidShapeMap, fv.shape())
0391                                    << " Layer " << lay << " Z " << t.z();
0392 #endif
0393       if (lay >= 0 && lay < maxLayer_) {
0394         ie_[lay]++;
0395         ze_[lay] += std::abs(HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z());
0396         if (thke_[lay] <= 0)
0397           thke_[lay] = dz;
0398         double rinHE = HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho() * cos(alp) - dy;
0399         double routHE = HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho() * cos(alp) + dy;
0400         rminHE_[lay] += rinHE;
0401         rmaxHE_[lay] += routHE;
0402         bool found = false;
0403         for (double k : php.zxHE) {
0404           if (std::abs(k - std::abs(HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z())) < 0.01) {
0405             found = true;
0406             break;
0407           }
0408         }
0409         if (!found) {
0410           php.zxHE.emplace_back(std::abs(HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z()));
0411           php.rhoxHE.emplace_back(HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho() * std::cos(t.phi()));
0412           php.dyHE.emplace_back(dy * std::cos(t.phi()));
0413           dx1 -= 0.5 * (HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho() - dy) * std::cos(t.phi()) * tan10deg;
0414           dx2 -= 0.5 * (HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho() + dy) * std::cos(t.phi()) * tan10deg;
0415           php.dx1HE.emplace_back(-dx1);
0416           php.dx2HE.emplace_back(-dx2);
0417           php.layHE.emplace_back(lay);
0418         }
0419       }
0420       if (copy[0] == kHELayer1_ || copy[0] == kHELayer2_) {
0421         int iz = copy[6];
0422         int fi = copy[4];
0423         unsigned int it1 = find(iz, ize_);
0424         if (it1 == ize_.size())
0425           ize_.emplace_back(iz);
0426         unsigned int it2 = find(fi, phie_);
0427         if (it2 == phie_.size())
0428           phie_.emplace_back(fi);
0429       }
0430     } else if (idet == 5) {
0431       // HF
0432       if (!hf) {
0433 #ifdef EDM_ML_DEBUG
0434         edm::LogVerbatim("HCalGeom") << "HF " << fv.name() << " Shape "
0435                                      << cms::dd::name(cms::DDSolidShapeMap, fv.shape()) << " Z "
0436                                      << HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z() << " with " << paras.size()
0437                                      << " Parameters";
0438         for (unsigned j = 0; j < paras.size(); j++)
0439           edm::LogVerbatim("HCalGeom") << "HF Parameter[" << j << "] = " << paras[j];
0440 #endif
0441         if (dd4hep::isA<dd4hep::Polycone>(fv.solid())) {
0442           int nz = (int)(paras.size()) - 3;
0443           dzVcal_ = 0.5 * HcalGeomParameters::k_ScaleFromDD4hepToG4 * (paras[nz] - paras[3]);
0444           hf = true;
0445         } else if (dd4hep::isA<dd4hep::Tube>(fv.solid()) || dd4hep::isA<dd4hep::ConeSegment>(fv.solid())) {
0446           dzVcal_ = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[2];
0447           hf = true;
0448         } else if (dd4hep::isA<dd4hep::Trap>(fv.solid())) {
0449           dzVcal_ = HcalGeomParameters::k_ScaleFromDD4hepToG4 * paras[0];
0450           hf = true;
0451         }
0452 #ifdef EDM_ML_DEBUG
0453         edm::LogVerbatim("HCalGeom") << "Sets for Detector " << idet << " for " << fv.name() << " Flag " << hf << " DZ "
0454                                      << dzVcal_;
0455 #endif
0456       }
0457 #ifdef EDM_ML_DEBUG
0458     } else {
0459       edm::LogVerbatim("HCalGeom") << "Unknown Detector " << idet << " for " << fv.name() << " Shape "
0460                                    << cms::dd::name(cms::DDSolidShapeMap, fv.shape()) << " R "
0461                                    << (HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.Rho()) << " Z "
0462                                    << (HcalGeomParameters::k_ScaleFromDD4hepToG4 * t.z());
0463 #endif
0464     }
0465   }
0466   loadfinal(php);
0467 }
0468 
0469 unsigned int HcalGeomParameters::find(int element, std::vector<int>& array) const {
0470   unsigned int id = array.size();
0471   for (unsigned int i = 0; i < array.size(); i++) {
0472     if (element == array[i]) {
0473       id = i;
0474       break;
0475     }
0476   }
0477   return id;
0478 }
0479 
0480 double HcalGeomParameters::getEta(double r, double z) const {
0481   double tmp = 0;
0482   if (z != 0)
0483     tmp = -log(tan(0.5 * atan(r / z)));
0484 #ifdef EDM_ML_DEBUG
0485   edm::LogVerbatim("HCalGeom") << "HcalGeomParameters::getEta " << r << " " << z << " ==> " << tmp;
0486 #endif
0487   return tmp;
0488 }
0489 
0490 void HcalGeomParameters::clear(HcalParameters& php) {
0491   php.rhoxHB.clear();
0492   php.zxHB.clear();
0493   php.dyHB.clear();
0494   php.dxHB.clear();
0495   php.layHB.clear();
0496   php.layHE.clear();
0497   php.zxHE.clear();
0498   php.rhoxHE.clear();
0499   php.dyHE.clear();
0500   php.dx1HE.clear();
0501   php.dx2HE.clear();
0502 
0503   // Initialize all variables
0504   nzHB_ = nmodHB_ = 0;
0505   nzHE_ = nmodHE_ = 0;
0506   for (int i = 0; i < 4; ++i)
0507     etaHO_[i] = 0;
0508   zVcal_ = dzVcal_ = dlShort_ = 0;
0509   rminHO_ = dzVcal_ = -1.;
0510   for (int i = 0; i < maxLayer_; ++i) {
0511     rb_.emplace_back(0.0);
0512     ze_.emplace_back(0.0);
0513     thkb_.emplace_back(-1.0);
0514     thke_.emplace_back(-1.0);
0515     ib_.emplace_back(0);
0516     ie_.emplace_back(0);
0517     rminHE_.emplace_back(0.0);
0518     rmaxHE_.emplace_back(0.0);
0519   }
0520 }
0521 
0522 void HcalGeomParameters::loadfinal(HcalParameters& php) {
0523   int ibmx = 0, iemx = 0;
0524   for (int i = 0; i < maxLayer_; i++) {
0525     if (ib_[i] > 0) {
0526       rb_[i] /= static_cast<double>(ib_[i]);
0527       ibmx = i + 1;
0528     }
0529     if (ie_[i] > 0) {
0530       ze_[i] /= static_cast<double>(ie_[i]);
0531       rminHE_[i] /= static_cast<double>(ie_[i]);
0532       rmaxHE_[i] /= static_cast<double>(ie_[i]);
0533       iemx = i + 1;
0534     }
0535 #ifdef EDM_ML_DEBUG
0536     edm::LogVerbatim("HCalGeom") << "Index " << i << " Barrel " << ib_[i] << " " << rb_[i] << " Endcap " << ie_[i]
0537                                  << " " << ze_[i] << ":" << rminHE_[i] << ":" << rmaxHE_[i];
0538 #endif
0539   }
0540   for (int i = 4; i >= 0; i--) {
0541     if (ib_[i] == 0) {
0542       rb_[i] = rb_[i + 1];
0543       thkb_[i] = thkb_[i + 1];
0544     }
0545     if (ie_[i] == 0) {
0546       ze_[i] = ze_[i + 1];
0547       thke_[i] = thke_[i + 1];
0548     }
0549 #ifdef EDM_ML_DEBUG
0550     if (ib_[i] == 0 || ie_[i] == 0)
0551       edm::LogVerbatim("HCalGeom") << "Index " << i << " Barrel " << ib_[i] << " " << rb_[i] << " Endcap " << ie_[i]
0552                                    << " " << ze_[i];
0553 #endif
0554   }
0555 
0556 #ifdef EDM_ML_DEBUG
0557   for (unsigned int k = 0; k < php.layHB.size(); ++k)
0558     edm::LogVerbatim("HCalGeom") << "HB: " << php.layHB[k] << " R " << rxb_[k] << " " << php.rhoxHB[k] << " Z "
0559                                  << php.zxHB[k] << " DY " << php.dyHB[k] << " DZ " << php.dxHB[k];
0560   for (unsigned int k = 0; k < php.layHE.size(); ++k)
0561     edm::LogVerbatim("HCalGeom") << "HE: " << php.layHE[k] << " R " << php.rhoxHE[k] << " Z " << php.zxHE[k]
0562                                  << " X1|X2 " << php.dx1HE[k] << "|" << php.dx2HE[k] << " DY " << php.dyHE[k];
0563   edm::LogVerbatim("HCalGeom") << "HcalGeomParameters: Maximum Layer for HB " << ibmx << " for HE " << iemx
0564                                << " extent " << dzVcal_;
0565 #endif
0566 
0567   if (ibmx > 0) {
0568     php.rHB.resize(ibmx);
0569     php.drHB.resize(ibmx);
0570     for (int i = 0; i < ibmx; i++) {
0571       php.rHB[i] = rb_[i];
0572       php.drHB[i] = thkb_[i];
0573 #ifdef EDM_ML_DEBUG
0574       edm::LogVerbatim("HCalGeom") << "HcalGeomParameters: php.rHB[" << i << "] = " << php.rHB[i] << " php.drHB[" << i
0575                                    << "] = " << php.drHB[i];
0576 #endif
0577     }
0578   }
0579   if (iemx > 0) {
0580     php.zHE.resize(iemx);
0581     php.dzHE.resize(iemx);
0582     for (int i = 0; i < iemx; i++) {
0583       php.zHE[i] = ze_[i];
0584       php.dzHE[i] = thke_[i];
0585 #ifdef EDM_ML_DEBUG
0586       edm::LogVerbatim("HCalGeom") << "HcalGeomParameters: php.zHE[" << i << "] = " << php.zHE[i] << " php.dzHE[" << i
0587                                    << "] = " << php.dzHE[i];
0588 #endif
0589     }
0590   }
0591 
0592   nzHB_ = static_cast<int>(izb_.size());
0593   nmodHB_ = static_cast<int>(phib_.size());
0594 #ifdef EDM_ML_DEBUG
0595   edm::LogVerbatim("HCalGeom") << "HcalGeomParameters::loadGeometry: " << nzHB_ << " barrel half-sectors";
0596   for (int i = 0; i < nzHB_; i++)
0597     edm::LogVerbatim("HCalGeom") << "Section " << i << " Copy number " << izb_[i];
0598   edm::LogVerbatim("HCalGeom") << "HcalGeomParameters::loadGeometry: " << nmodHB_ << " barrel modules";
0599   for (int i = 0; i < nmodHB_; i++)
0600     edm::LogVerbatim("HCalGeom") << "Module " << i << " Copy number " << phib_[i];
0601 #endif
0602 
0603   nzHE_ = static_cast<int>(ize_.size());
0604   nmodHE_ = static_cast<int>(phie_.size());
0605 #ifdef EDM_ML_DEBUG
0606   edm::LogVerbatim("HCalGeom") << "HcalGeomParameters::loadGeometry: " << nzHE_ << " endcap half-sectors";
0607   for (int i = 0; i < nzHE_; i++)
0608     edm::LogVerbatim("HCalGeom") << "Section " << i << " Copy number " << ize_[i];
0609   edm::LogVerbatim("HCalGeom") << "HcalGeomParameters::loadGeometry: " << nmodHE_ << " endcap modules";
0610   for (int i = 0; i < nmodHE_; i++)
0611     edm::LogVerbatim("HCalGeom") << "Module " << i << " Copy number " << phie_[i];
0612 #endif
0613 
0614 #ifdef EDM_ML_DEBUG
0615   edm::LogVerbatim("HCalGeom") << "HO has Z of size " << php.zHO.size();
0616   for (unsigned int kk = 0; kk < php.zHO.size(); kk++)
0617     edm::LogVerbatim("HCalGeom") << "ZHO[" << kk << "] = " << php.zHO[kk];
0618 #endif
0619   if (ibmx > 17 && php.zHO.size() > 4) {
0620     rminHO_ = php.rHB[17] - 100.0;
0621     etaHO_[0] = getEta(0.5 * (php.rHB[17] + php.rHB[18]), php.zHO[1]);
0622     etaHO_[1] = getEta(php.rHB[18] + php.drHB[18], php.zHO[2]);
0623     etaHO_[2] = getEta(php.rHB[18] - php.drHB[18], php.zHO[3]);
0624     etaHO_[3] = getEta(php.rHB[18] + php.drHB[18], php.zHO[4]);
0625   }
0626 #ifdef EDM_ML_DEBUG
0627   edm::LogVerbatim("HCalGeom") << "HO Eta boundaries " << etaHO_[0] << " " << etaHO_[1] << " " << etaHO_[2] << " "
0628                                << etaHO_[3];
0629   edm::LogVerbatim("HCalGeom") << "HO Parameters " << rminHO_ << " " << php.zHO.size();
0630   for (unsigned int i = 0; i < php.zHO.size(); ++i)
0631     edm::LogVerbatim("HCalGeom") << " zho[" << i << "] = " << php.zHO[i];
0632 #endif
0633 }