Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-06-22 02:23:46

0001 #include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.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 "DataFormats/Math/interface/Point3D.h"
0008 #include "DetectorDescription/Core/interface/DDConstant.h"
0009 #include "DetectorDescription/Core/interface/DDFilter.h"
0010 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0011 #include "DetectorDescription/Core/interface/DDSolid.h"
0012 #include "DetectorDescription/Core/interface/DDValue.h"
0013 #include "DetectorDescription/Core/interface/DDutils.h"
0014 #include "DetectorDescription/RegressionTest/interface/DDErrorDetection.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalProperty.h"
0016 #include "Geometry/HGCalCommonData/interface/HGCalTileIndex.h"
0017 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0018 #include "Geometry/HGCalCommonData/interface/HGCalWaferIndex.h"
0019 #include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
0020 #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h"
0021 
0022 #include <algorithm>
0023 #include <sstream>
0024 #include <unordered_set>
0025 
0026 //#define EDM_ML_DEBUG
0027 using namespace geant_units::operators;
0028 
0029 const double tolerance = 0.001;
0030 const double tolmin = 1.e-20;
0031 
0032 HGCalGeomParameters::HGCalGeomParameters() : sqrt3_(std::sqrt(3.0)) {
0033 #ifdef EDM_ML_DEBUG
0034   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters::HGCalGeomParameters "
0035                                 << "constructor";
0036 #endif
0037 }
0038 
0039 void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv,
0040                                               HGCalParameters& php,
0041                                               const std::string& sdTag1,
0042                                               const DDCompactView* cpv,
0043                                               const std::string& sdTag2,
0044                                               const std::string& sdTag3,
0045                                               HGCalGeometryMode::WaferMode mode) {
0046   DDFilteredView fv = _fv;
0047   bool dodet(true);
0048   std::map<int, HGCalGeomParameters::layerParameters> layers;
0049   std::vector<HGCalParameters::hgtrform> trforms;
0050   std::vector<bool> trformUse;
0051 
0052   while (dodet) {
0053     const DDSolid& sol = fv.logicalPart().solid();
0054     // Layers first
0055     std::vector<int> copy = fv.copyNumbers();
0056     int nsiz = static_cast<int>(copy.size());
0057     int lay = (nsiz > 0) ? copy[nsiz - 1] : 0;
0058     int zp = (nsiz > 2) ? copy[nsiz - 3] : -1;
0059     if (zp != 1)
0060       zp = -1;
0061     if (lay == 0) {
0062       throw cms::Exception("DDException") << "Funny layer # " << lay << " zp " << zp << " in " << nsiz << " components";
0063     } else {
0064       if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0065         php.layer_.emplace_back(lay);
0066       auto itr = layers.find(lay);
0067       if (itr == layers.end()) {
0068         double rin(0), rout(0);
0069         double zz = HGCalParameters::k_ScaleFromDDD * fv.translation().Z();
0070         if ((sol.shape() == DDSolidShape::ddpolyhedra_rz) || (sol.shape() == DDSolidShape::ddpolyhedra_rrz)) {
0071           const DDPolyhedra& polyhedra = static_cast<DDPolyhedra>(sol);
0072           const std::vector<double>& rmin = polyhedra.rMinVec();
0073           const std::vector<double>& rmax = polyhedra.rMaxVec();
0074           rin = 0.5 * HGCalParameters::k_ScaleFromDDD * (rmin[0] + rmin[1]);
0075           rout = 0.5 * HGCalParameters::k_ScaleFromDDD * (rmax[0] + rmax[1]);
0076         } else if (sol.shape() == DDSolidShape::ddtubs) {
0077           const DDTubs& tube = static_cast<DDTubs>(sol);
0078           rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
0079           rout = HGCalParameters::k_ScaleFromDDD * tube.rOut();
0080         }
0081         HGCalGeomParameters::layerParameters laypar(rin, rout, zz);
0082         layers[lay] = laypar;
0083       }
0084       DD3Vector x, y, z;
0085       fv.rotation().GetComponents(x, y, z);
0086       const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0087       const CLHEP::HepRotation hr(rotation);
0088       double xx = HGCalParameters::k_ScaleFromDDD * fv.translation().X();
0089       if (std::abs(xx) < tolerance)
0090         xx = 0;
0091       double yy = HGCalParameters::k_ScaleFromDDD * fv.translation().Y();
0092       if (std::abs(yy) < tolerance)
0093         yy = 0;
0094       double zz = HGCalParameters::k_ScaleFromDDD * fv.translation().Z();
0095       const CLHEP::Hep3Vector h3v(xx, yy, zz);
0096       HGCalParameters::hgtrform mytrf;
0097       mytrf.zp = zp;
0098       mytrf.lay = lay;
0099       mytrf.sec = 0;
0100       mytrf.subsec = 0;
0101       mytrf.h3v = h3v;
0102       mytrf.hr = hr;
0103       trforms.emplace_back(mytrf);
0104       trformUse.emplace_back(false);
0105     }
0106     dodet = fv.next();
0107   }
0108 
0109   // Then wafers
0110   // This assumes layers are build starting from 1 (which on 25 Jan 2016, they
0111   // were) to ensure that new copy numbers are always added to the end of the
0112   // list.
0113   std::unordered_map<int32_t, int32_t> copies;
0114   HGCalParameters::layer_map copiesInLayers(layers.size() + 1);
0115   std::vector<int32_t> wafer2copy;
0116   std::vector<HGCalGeomParameters::cellParameters> wafers;
0117   std::string attribute = "Volume";
0118   DDValue val1(attribute, sdTag2, 0.0);
0119   DDSpecificsMatchesValueFilter filter1{val1};
0120   DDFilteredView fv1(*cpv, filter1);
0121   bool ok = fv1.firstChild();
0122   if (!ok) {
0123     throw cms::Exception("DDException") << "Attribute " << val1 << " not found but needed.";
0124   } else {
0125     dodet = true;
0126     std::unordered_set<std::string> names;
0127     while (dodet) {
0128       const DDSolid& sol = fv1.logicalPart().solid();
0129       const std::string& name = fv1.logicalPart().name().name();
0130       std::vector<int> copy = fv1.copyNumbers();
0131       int nsiz = static_cast<int>(copy.size());
0132       int wafer = (nsiz > 0) ? copy[nsiz - 1] : 0;
0133       int layer = (nsiz > 1) ? copy[nsiz - 2] : 0;
0134       if (nsiz < 2) {
0135         throw cms::Exception("DDException") << "Funny wafer # " << wafer << " in " << nsiz << " components";
0136       } else if (layer > static_cast<int>(layers.size())) {
0137         edm::LogWarning("HGCalGeom") << "Funny wafer # " << wafer << " Layer " << layer << ":" << layers.size()
0138                                      << " among " << nsiz << " components";
0139       } else {
0140         auto itr = copies.find(wafer);
0141         auto cpy = copiesInLayers[layer].find(wafer);
0142         if (itr != copies.end() && cpy == copiesInLayers[layer].end()) {
0143           copiesInLayers[layer][wafer] = itr->second;
0144         }
0145         if (itr == copies.end()) {
0146           copies[wafer] = wafer2copy.size();
0147           copiesInLayers[layer][wafer] = wafer2copy.size();
0148           double xx = HGCalParameters::k_ScaleFromDDD * fv1.translation().X();
0149           if (std::abs(xx) < tolerance)
0150             xx = 0;
0151           double yy = HGCalParameters::k_ScaleFromDDD * fv1.translation().Y();
0152           if (std::abs(yy) < tolerance)
0153             yy = 0;
0154           wafer2copy.emplace_back(wafer);
0155           GlobalPoint p(xx, yy, HGCalParameters::k_ScaleFromDDD * fv1.translation().Z());
0156           HGCalGeomParameters::cellParameters cell(false, wafer, p);
0157           wafers.emplace_back(cell);
0158           if (names.count(name) == 0) {
0159             std::vector<double> zv, rv;
0160             if (mode == HGCalGeometryMode::Polyhedra) {
0161               const DDPolyhedra& polyhedra = static_cast<DDPolyhedra>(sol);
0162               zv = polyhedra.zVec();
0163               rv = polyhedra.rMaxVec();
0164             } else {
0165               const DDExtrudedPolygon& polygon = static_cast<DDExtrudedPolygon>(sol);
0166               zv = polygon.zVec();
0167               rv = polygon.xVec();
0168             }
0169             php.waferR_ = 2.0 * HGCalParameters::k_ScaleFromDDDToG4 * rv[0] * tan30deg_;
0170             php.waferSize_ = HGCalParameters::k_ScaleFromDDD * rv[0];
0171             double dz = 0.5 * HGCalParameters::k_ScaleFromDDDToG4 * (zv[1] - zv[0]);
0172 #ifdef EDM_ML_DEBUG
0173             edm::LogVerbatim("HGCalGeom")
0174                 << "Mode " << mode << " R " << php.waferSize_ << ":" << php.waferR_ << " z " << dz;
0175 #endif
0176             HGCalParameters::hgtrap mytr;
0177             mytr.lay = 1;
0178             mytr.bl = php.waferR_;
0179             mytr.tl = php.waferR_;
0180             mytr.h = php.waferR_;
0181             mytr.dz = dz;
0182             mytr.alpha = 0.0;
0183             mytr.cellSize = waferSize_;
0184             php.fillModule(mytr, false);
0185             names.insert(name);
0186           }
0187         }
0188       }
0189       dodet = fv1.next();
0190     }
0191   }
0192 
0193   // Finally the cells
0194   std::map<int, int> wafertype;
0195   std::map<int, HGCalGeomParameters::cellParameters> cellsf, cellsc;
0196   DDValue val2(attribute, sdTag3, 0.0);
0197   DDSpecificsMatchesValueFilter filter2{val2};
0198   DDFilteredView fv2(*cpv, filter2);
0199   ok = fv2.firstChild();
0200   if (!ok) {
0201     throw cms::Exception("DDException") << "Attribute " << val2 << " not found but needed.";
0202   } else {
0203     dodet = true;
0204     while (dodet) {
0205       const DDSolid& sol = fv2.logicalPart().solid();
0206       const std::string& name = sol.name().name();
0207       std::vector<int> copy = fv2.copyNumbers();
0208       int nsiz = static_cast<int>(copy.size());
0209       int cellx = (nsiz > 0) ? copy[nsiz - 1] : 0;
0210       int wafer = (nsiz > 1) ? copy[nsiz - 2] : 0;
0211       int cell = HGCalTypes::getUnpackedCell6(cellx);
0212       int type = HGCalTypes::getUnpackedCellType6(cellx);
0213       if (type != 1 && type != 2) {
0214         throw cms::Exception("DDException")
0215             << "Funny cell # " << cell << " type " << type << " in " << nsiz << " components";
0216       } else {
0217         auto ktr = wafertype.find(wafer);
0218         if (ktr == wafertype.end())
0219           wafertype[wafer] = type;
0220         bool newc(false);
0221         std::map<int, HGCalGeomParameters::cellParameters>::iterator itr;
0222         double cellsize = php.cellSize_[0];
0223         if (type == 1) {
0224           itr = cellsf.find(cell);
0225           newc = (itr == cellsf.end());
0226         } else {
0227           itr = cellsc.find(cell);
0228           newc = (itr == cellsc.end());
0229           cellsize = php.cellSize_[1];
0230         }
0231         if (newc) {
0232           bool half = (name.find("Half") != std::string::npos);
0233           double xx = HGCalParameters::k_ScaleFromDDD * fv2.translation().X();
0234           double yy = HGCalParameters::k_ScaleFromDDD * fv2.translation().Y();
0235           if (half) {
0236             math::XYZPointD p1(-2.0 * cellsize / 9.0, 0, 0);
0237             math::XYZPointD p2 = fv2.rotation()(p1);
0238             xx += (HGCalParameters::k_ScaleFromDDD * (p2.X()));
0239             yy += (HGCalParameters::k_ScaleFromDDD * (p2.Y()));
0240 #ifdef EDM_ML_DEBUG
0241             if (std::abs(p2.X()) < HGCalParameters::tol)
0242               p2.SetX(0.0);
0243             if (std::abs(p2.Z()) < HGCalParameters::tol)
0244               p2.SetZ(0.0);
0245             edm::LogVerbatim("HGCalGeom") << "Wafer " << wafer << " Type " << type << " Cell " << cellx << " local "
0246                                           << xx << ":" << yy << " new " << p1 << ":" << p2;
0247 #endif
0248           }
0249           HGCalGeomParameters::cellParameters cp(half, wafer, GlobalPoint(xx, yy, 0));
0250           if (type == 1) {
0251             cellsf[cell] = cp;
0252           } else {
0253             cellsc[cell] = cp;
0254           }
0255         }
0256       }
0257       dodet = fv2.next();
0258     }
0259   }
0260 
0261   loadGeometryHexagon(
0262       layers, trforms, trformUse, copies, copiesInLayers, wafer2copy, wafers, wafertype, cellsf, cellsc, php);
0263 }
0264 
0265 void HGCalGeomParameters::loadGeometryHexagon(const cms::DDCompactView* cpv,
0266                                               HGCalParameters& php,
0267                                               const std::string& sdTag1,
0268                                               const std::string& sdTag2,
0269                                               const std::string& sdTag3,
0270                                               HGCalGeometryMode::WaferMode mode) {
0271   const cms::DDFilter filter("Volume", sdTag1);
0272   cms::DDFilteredView fv((*cpv), filter);
0273   std::map<int, HGCalGeomParameters::layerParameters> layers;
0274   std::vector<HGCalParameters::hgtrform> trforms;
0275   std::vector<bool> trformUse;
0276   std::vector<std::pair<int, int> > trused;
0277 
0278   while (fv.firstChild()) {
0279     const std::vector<double>& pars = fv.parameters();
0280     // Layers first
0281     std::vector<int> copy = fv.copyNos();
0282     int nsiz = static_cast<int>(copy.size());
0283     int lay = (nsiz > 0) ? copy[0] : 0;
0284     int zp = (nsiz > 2) ? copy[2] : -1;
0285     if (zp != 1)
0286       zp = -1;
0287     if (lay == 0) {
0288       throw cms::Exception("DDException") << "Funny layer # " << lay << " zp " << zp << " in " << nsiz << " components";
0289     } else {
0290       if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0291         php.layer_.emplace_back(lay);
0292       auto itr = layers.find(lay);
0293       double zz = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z();
0294       if (itr == layers.end()) {
0295         double rin(0), rout(0);
0296         if (dd4hep::isA<dd4hep::Polyhedra>(fv.solid())) {
0297           rin = 0.5 * HGCalParameters::k_ScaleFromDD4hep * (pars[5] + pars[8]);
0298           rout = 0.5 * HGCalParameters::k_ScaleFromDD4hep * (pars[6] + pars[9]);
0299         } else if (dd4hep::isA<dd4hep::Tube>(fv.solid())) {
0300           dd4hep::Tube tubeSeg(fv.solid());
0301           rin = HGCalParameters::k_ScaleFromDD4hep * tubeSeg.rMin();
0302           rout = HGCalParameters::k_ScaleFromDD4hep * tubeSeg.rMax();
0303         }
0304         HGCalGeomParameters::layerParameters laypar(rin, rout, zz);
0305         layers[lay] = laypar;
0306       }
0307       std::pair<int, int> layz(lay, zp);
0308       if (std::find(trused.begin(), trused.end(), layz) == trused.end()) {
0309         trused.emplace_back(layz);
0310         DD3Vector x, y, z;
0311         fv.rotation().GetComponents(x, y, z);
0312         const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0313         const CLHEP::HepRotation hr(rotation);
0314         double xx = HGCalParameters::k_ScaleFromDD4hep * fv.translation().X();
0315         if (std::abs(xx) < tolerance)
0316           xx = 0;
0317         double yy = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Y();
0318         if (std::abs(yy) < tolerance)
0319           yy = 0;
0320         double zz = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z();
0321         const CLHEP::Hep3Vector h3v(xx, yy, zz);
0322         HGCalParameters::hgtrform mytrf;
0323         mytrf.zp = zp;
0324         mytrf.lay = lay;
0325         mytrf.sec = 0;
0326         mytrf.subsec = 0;
0327         mytrf.h3v = h3v;
0328         mytrf.hr = hr;
0329         trforms.emplace_back(mytrf);
0330         trformUse.emplace_back(false);
0331       }
0332     }
0333   }
0334 
0335   // Then wafers
0336   // This assumes layers are build starting from 1 (which on 25 Jan 2016, they
0337   // were) to ensure that new copy numbers are always added to the end of the
0338   // list.
0339   std::unordered_map<int32_t, int32_t> copies;
0340   HGCalParameters::layer_map copiesInLayers(layers.size() + 1);
0341   std::vector<int32_t> wafer2copy;
0342   std::vector<HGCalGeomParameters::cellParameters> wafers;
0343   const cms::DDFilter filter1("Volume", sdTag2);
0344   cms::DDFilteredView fv1((*cpv), filter1);
0345   bool ok = fv1.firstChild();
0346   if (!ok) {
0347     throw cms::Exception("DDException") << "Attribute " << sdTag2 << " not found but needed.";
0348   } else {
0349     bool dodet = true;
0350     std::unordered_set<std::string> names;
0351     while (dodet) {
0352       const std::string name = static_cast<std::string>(fv1.name());
0353       std::vector<int> copy = fv1.copyNos();
0354       int nsiz = static_cast<int>(copy.size());
0355       int wafer = (nsiz > 0) ? copy[0] : 0;
0356       int layer = (nsiz > 1) ? copy[1] : 0;
0357       if (nsiz < 2) {
0358         throw cms::Exception("DDException") << "Funny wafer # " << wafer << " in " << nsiz << " components";
0359       } else if (layer > static_cast<int>(layers.size())) {
0360         edm::LogWarning("HGCalGeom") << "Funny wafer # " << wafer << " Layer " << layer << ":" << layers.size()
0361                                      << " among " << nsiz << " components";
0362       } else {
0363         auto itr = copies.find(wafer);
0364         auto cpy = copiesInLayers[layer].find(wafer);
0365         if (itr != copies.end() && cpy == copiesInLayers[layer].end()) {
0366           copiesInLayers[layer][wafer] = itr->second;
0367         }
0368         if (itr == copies.end()) {
0369           copies[wafer] = wafer2copy.size();
0370           copiesInLayers[layer][wafer] = wafer2copy.size();
0371           double xx = HGCalParameters::k_ScaleFromDD4hep * fv1.translation().X();
0372           if (std::abs(xx) < tolerance)
0373             xx = 0;
0374           double yy = HGCalParameters::k_ScaleFromDD4hep * fv1.translation().Y();
0375           if (std::abs(yy) < tolerance)
0376             yy = 0;
0377           wafer2copy.emplace_back(wafer);
0378           GlobalPoint p(xx, yy, HGCalParameters::k_ScaleFromDD4hep * fv1.translation().Z());
0379           HGCalGeomParameters::cellParameters cell(false, wafer, p);
0380           wafers.emplace_back(cell);
0381           if (names.count(name) == 0) {
0382             double zv[2], rv;
0383             const std::vector<double>& pars = fv1.parameters();
0384             if (mode == HGCalGeometryMode::Polyhedra) {
0385               zv[0] = pars[4];
0386               zv[1] = pars[7];
0387               rv = pars[6];
0388             } else {
0389               zv[0] = pars[3];
0390               zv[1] = pars[9];
0391               rv = pars[4];
0392             }
0393             php.waferR_ = 2.0 * HGCalParameters::k_ScaleFromDD4hepToG4 * rv * tan30deg_;
0394             php.waferSize_ = HGCalParameters::k_ScaleFromDD4hep * rv;
0395             double dz = 0.5 * HGCalParameters::k_ScaleFromDD4hepToG4 * (zv[1] - zv[0]);
0396 #ifdef EDM_ML_DEBUG
0397             edm::LogVerbatim("HGCalGeom")
0398                 << "Mode " << mode << " R " << php.waferSize_ << ":" << php.waferR_ << " z " << dz;
0399 #endif
0400             HGCalParameters::hgtrap mytr;
0401             mytr.lay = 1;
0402             mytr.bl = php.waferR_;
0403             mytr.tl = php.waferR_;
0404             mytr.h = php.waferR_;
0405             mytr.dz = dz;
0406             mytr.alpha = 0.0;
0407             mytr.cellSize = waferSize_;
0408             php.fillModule(mytr, false);
0409             names.insert(name);
0410           }
0411         }
0412       }
0413       dodet = fv1.firstChild();
0414     }
0415   }
0416 
0417   // Finally the cells
0418   std::map<int, int> wafertype;
0419   std::map<int, HGCalGeomParameters::cellParameters> cellsf, cellsc;
0420   const cms::DDFilter filter2("Volume", sdTag3);
0421   cms::DDFilteredView fv2((*cpv), filter2);
0422   ok = fv2.firstChild();
0423   if (!ok) {
0424     throw cms::Exception("DDException") << "Attribute " << sdTag3 << " not found but needed.";
0425   } else {
0426     bool dodet = true;
0427     while (dodet) {
0428       const std::string name = static_cast<std::string>(fv2.name());
0429       std::vector<int> copy = fv2.copyNos();
0430       int nsiz = static_cast<int>(copy.size());
0431       int cellx = (nsiz > 0) ? copy[0] : 0;
0432       int wafer = (nsiz > 1) ? copy[1] : 0;
0433       int cell = HGCalTypes::getUnpackedCell6(cellx);
0434       int type = HGCalTypes::getUnpackedCellType6(cellx);
0435       if (type != 1 && type != 2) {
0436         throw cms::Exception("DDException")
0437             << "Funny cell # " << cell << " type " << type << " in " << nsiz << " components";
0438       } else {
0439         auto ktr = wafertype.find(wafer);
0440         if (ktr == wafertype.end())
0441           wafertype[wafer] = type;
0442         bool newc(false);
0443         std::map<int, HGCalGeomParameters::cellParameters>::iterator itr;
0444         double cellsize = php.cellSize_[0];
0445         if (type == 1) {
0446           itr = cellsf.find(cell);
0447           newc = (itr == cellsf.end());
0448         } else {
0449           itr = cellsc.find(cell);
0450           newc = (itr == cellsc.end());
0451           cellsize = php.cellSize_[1];
0452         }
0453         if (newc) {
0454           bool half = (name.find("Half") != std::string::npos);
0455           double xx = HGCalParameters::k_ScaleFromDD4hep * fv2.translation().X();
0456           double yy = HGCalParameters::k_ScaleFromDD4hep * fv2.translation().Y();
0457           if (half) {
0458             math::XYZPointD p1(-2.0 * cellsize / 9.0, 0, 0);
0459             math::XYZPointD p2 = fv2.rotation()(p1);
0460             xx += (HGCalParameters::k_ScaleFromDDD * (p2.X()));
0461             yy += (HGCalParameters::k_ScaleFromDDD * (p2.Y()));
0462 #ifdef EDM_ML_DEBUG
0463             if (std::abs(p2.X()) < HGCalParameters::tol)
0464               p2.SetX(0.0);
0465             if (std::abs(p2.Z()) < HGCalParameters::tol)
0466               p2.SetZ(0.0);
0467             edm::LogVerbatim("HGCalGeom") << "Wafer " << wafer << " Type " << type << " Cell " << cellx << " local "
0468                                           << xx << ":" << yy << " new " << p1 << ":" << p2;
0469 #endif
0470           }
0471           HGCalGeomParameters::cellParameters cp(half, wafer, GlobalPoint(xx, yy, 0));
0472           if (type == 1) {
0473             cellsf[cell] = cp;
0474           } else {
0475             cellsc[cell] = cp;
0476           }
0477         }
0478       }
0479       dodet = fv2.firstChild();
0480     }
0481   }
0482 
0483   loadGeometryHexagon(
0484       layers, trforms, trformUse, copies, copiesInLayers, wafer2copy, wafers, wafertype, cellsf, cellsc, php);
0485 }
0486 
0487 void HGCalGeomParameters::loadGeometryHexagon(const std::map<int, HGCalGeomParameters::layerParameters>& layers,
0488                                               std::vector<HGCalParameters::hgtrform>& trforms,
0489                                               std::vector<bool>& trformUse,
0490                                               const std::unordered_map<int32_t, int32_t>& copies,
0491                                               const HGCalParameters::layer_map& copiesInLayers,
0492                                               const std::vector<int32_t>& wafer2copy,
0493                                               const std::vector<HGCalGeomParameters::cellParameters>& wafers,
0494                                               const std::map<int, int>& wafertype,
0495                                               const std::map<int, HGCalGeomParameters::cellParameters>& cellsf,
0496                                               const std::map<int, HGCalGeomParameters::cellParameters>& cellsc,
0497                                               HGCalParameters& php) {
0498   if (((cellsf.size() + cellsc.size()) == 0) || (wafers.empty()) || (layers.empty())) {
0499     throw cms::Exception("DDException") << "HGCalGeomParameters: mismatch between geometry and specpar: cells "
0500                                         << cellsf.size() << ":" << cellsc.size() << " wafers " << wafers.size()
0501                                         << " layers " << layers.size();
0502   }
0503 
0504   for (unsigned int i = 0; i < layers.size(); ++i) {
0505     for (auto& layer : layers) {
0506       if (layer.first == static_cast<int>(i + php.firstLayer_)) {
0507         php.layerIndex_.emplace_back(i);
0508         php.rMinLayHex_.emplace_back(layer.second.rmin);
0509         php.rMaxLayHex_.emplace_back(layer.second.rmax);
0510         php.zLayerHex_.emplace_back(layer.second.zpos);
0511         break;
0512       }
0513     }
0514   }
0515 
0516   for (unsigned int i = 0; i < php.layer_.size(); ++i) {
0517     for (unsigned int i1 = 0; i1 < trforms.size(); ++i1) {
0518       if (!trformUse[i1] && php.layerGroup_[trforms[i1].lay - 1] == static_cast<int>(i + 1)) {
0519         trforms[i1].h3v *= static_cast<double>(HGCalParameters::k_ScaleFromDDD);
0520         trforms[i1].lay = (i + 1);
0521         trformUse[i1] = true;
0522         php.fillTrForm(trforms[i1]);
0523         int nz(1);
0524         for (unsigned int i2 = i1 + 1; i2 < trforms.size(); ++i2) {
0525           if (!trformUse[i2] && trforms[i2].zp == trforms[i1].zp &&
0526               php.layerGroup_[trforms[i2].lay - 1] == static_cast<int>(i + 1)) {
0527             php.addTrForm(trforms[i2].h3v);
0528             nz++;
0529             trformUse[i2] = true;
0530           }
0531         }
0532         if (nz > 0) {
0533           php.scaleTrForm(double(1.0 / nz));
0534         }
0535       }
0536     }
0537   }
0538 
0539   double rmin = HGCalParameters::k_ScaleFromDDD * php.waferR_;
0540   for (unsigned i = 0; i < wafer2copy.size(); ++i) {
0541     php.waferCopy_.emplace_back(wafer2copy[i]);
0542     php.waferPosX_.emplace_back(wafers[i].xyz.x());
0543     php.waferPosY_.emplace_back(wafers[i].xyz.y());
0544     auto ktr = wafertype.find(wafer2copy[i]);
0545     int typet = (ktr == wafertype.end()) ? 0 : (ktr->second);
0546     php.waferTypeT_.emplace_back(typet);
0547     double r = wafers[i].xyz.perp();
0548     int type(3);
0549     for (int k = 1; k < 4; ++k) {
0550       if ((r + rmin) <= php.boundR_[k]) {
0551         type = k;
0552         break;
0553       }
0554     }
0555     php.waferTypeL_.emplace_back(type);
0556   }
0557   php.copiesInLayers_ = copiesInLayers;
0558   php.nSectors_ = static_cast<int>(php.waferCopy_.size());
0559 
0560   std::vector<HGCalGeomParameters::cellParameters>::const_iterator itrf = wafers.end();
0561   for (unsigned int i = 0; i < cellsf.size(); ++i) {
0562     auto itr = cellsf.find(i);
0563     if (itr == cellsf.end()) {
0564       throw cms::Exception("DDException") << "HGCalGeomParameters: missing info for fine cell number " << i;
0565     } else {
0566       double xx = (itr->second).xyz.x();
0567       double yy = (itr->second).xyz.y();
0568       int waf = (itr->second).wafer;
0569       std::pair<double, double> xy = cellPosition(wafers, itrf, waf, xx, yy);
0570       php.cellFineX_.emplace_back(xy.first);
0571       php.cellFineY_.emplace_back(xy.second);
0572       php.cellFineHalf_.emplace_back((itr->second).half);
0573     }
0574   }
0575   itrf = wafers.end();
0576   for (unsigned int i = 0; i < cellsc.size(); ++i) {
0577     auto itr = cellsc.find(i);
0578     if (itr == cellsc.end()) {
0579       throw cms::Exception("DDException") << "HGCalGeomParameters: missing info for coarse cell number " << i;
0580     } else {
0581       double xx = (itr->second).xyz.x();
0582       double yy = (itr->second).xyz.y();
0583       int waf = (itr->second).wafer;
0584       std::pair<double, double> xy = cellPosition(wafers, itrf, waf, xx, yy);
0585       php.cellCoarseX_.emplace_back(xy.first);
0586       php.cellCoarseY_.emplace_back(xy.second);
0587       php.cellCoarseHalf_.emplace_back((itr->second).half);
0588     }
0589   }
0590   int depth(0);
0591   for (unsigned int i = 0; i < php.layerGroup_.size(); ++i) {
0592     bool first(true);
0593     for (unsigned int k = 0; k < php.layerGroup_.size(); ++k) {
0594       if (php.layerGroup_[k] == static_cast<int>(i + 1)) {
0595         if (first) {
0596           php.depth_.emplace_back(i + 1);
0597           php.depthIndex_.emplace_back(depth);
0598           php.depthLayerF_.emplace_back(k);
0599           ++depth;
0600           first = false;
0601         }
0602       }
0603     }
0604   }
0605   HGCalParameters::hgtrap mytr = php.getModule(0, false);
0606   mytr.bl *= HGCalParameters::k_ScaleFromDDD;
0607   mytr.tl *= HGCalParameters::k_ScaleFromDDD;
0608   mytr.h *= HGCalParameters::k_ScaleFromDDD;
0609   mytr.dz *= HGCalParameters::k_ScaleFromDDD;
0610   mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
0611   double dz = mytr.dz;
0612   php.fillModule(mytr, true);
0613   mytr.dz = 2 * dz;
0614   php.fillModule(mytr, true);
0615   mytr.dz = 3 * dz;
0616   php.fillModule(mytr, true);
0617 #ifdef EDM_ML_DEBUG
0618   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.zLayerHex_.size() << " layers";
0619   for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
0620     int k = php.layerIndex_[i];
0621     edm::LogVerbatim("HGCalGeom") << "Layer[" << i << ":" << k << ":" << php.layer_[k]
0622                                   << "] with r = " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
0623                                   << " at z = " << php.zLayerHex_[i];
0624   }
0625   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters has " << php.depthIndex_.size() << " depths";
0626   for (unsigned int i = 0; i < php.depthIndex_.size(); ++i) {
0627     int k = php.depthIndex_[i];
0628     edm::LogVerbatim("HGCalGeom") << "Reco Layer[" << i << ":" << k << "]  First Layer " << php.depthLayerF_[i]
0629                                   << " Depth " << php.depth_[k];
0630   }
0631   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.nSectors_ << " wafers";
0632   for (unsigned int i = 0; i < php.waferCopy_.size(); ++i)
0633     edm::LogVerbatim("HGCalGeom") << "Wafer[" << i << ": " << php.waferCopy_[i] << "] type " << php.waferTypeL_[i]
0634                                   << ":" << php.waferTypeT_[i] << " at (" << php.waferPosX_[i] << ","
0635                                   << php.waferPosY_[i] << ",0)";
0636   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer radius " << php.waferR_ << " and dimensions of the "
0637                                 << "wafers:";
0638   edm::LogVerbatim("HGCalGeom") << "Sim[0] " << php.moduleLayS_[0] << " dx " << php.moduleBlS_[0] << ":"
0639                                 << php.moduleTlS_[0] << " dy " << php.moduleHS_[0] << " dz " << php.moduleDzS_[0]
0640                                 << " alpha " << php.moduleAlphaS_[0];
0641   for (unsigned int k = 0; k < php.moduleLayR_.size(); ++k)
0642     edm::LogVerbatim("HGCalGeom") << "Rec[" << k << "] " << php.moduleLayR_[k] << " dx " << php.moduleBlR_[k] << ":"
0643                                   << php.moduleTlR_[k] << " dy " << php.moduleHR_[k] << " dz " << php.moduleDzR_[k]
0644                                   << " alpha " << php.moduleAlphaR_[k];
0645   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.cellFineX_.size() << " fine cells in a  wafer";
0646   for (unsigned int i = 0; i < php.cellFineX_.size(); ++i)
0647     edm::LogVerbatim("HGCalGeom") << "Fine Cell[" << i << "] at (" << php.cellFineX_[i] << "," << php.cellFineY_[i]
0648                                   << ",0)";
0649   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.cellCoarseX_.size()
0650                                 << " coarse cells in a wafer";
0651   for (unsigned int i = 0; i < php.cellCoarseX_.size(); ++i)
0652     edm::LogVerbatim("HGCalGeom") << "Coarse Cell[" << i << "] at (" << php.cellCoarseX_[i] << ","
0653                                   << php.cellCoarseY_[i] << ",0)";
0654   edm::LogVerbatim("HGCalGeom") << "Obtained " << php.trformIndex_.size() << " transformation matrices";
0655   for (unsigned int k = 0; k < php.trformIndex_.size(); ++k) {
0656     edm::LogVerbatim("HGCalGeom") << "Matrix[" << k << "] (" << std::hex << php.trformIndex_[k] << std::dec
0657                                   << ") Translation (" << php.trformTranX_[k] << ", " << php.trformTranY_[k] << ", "
0658                                   << php.trformTranZ_[k] << " Rotation (" << php.trformRotXX_[k] << ", "
0659                                   << php.trformRotYX_[k] << ", " << php.trformRotZX_[k] << ", " << php.trformRotXY_[k]
0660                                   << ", " << php.trformRotYY_[k] << ", " << php.trformRotZY_[k] << ", "
0661                                   << php.trformRotXZ_[k] << ", " << php.trformRotYZ_[k] << ", " << php.trformRotZZ_[k]
0662                                   << ")";
0663   }
0664   edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
0665   for (unsigned int k = 0; k < php.copiesInLayers_.size(); ++k) {
0666     const auto& theModules = php.copiesInLayers_[k];
0667     edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
0668     int k2(0);
0669     for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
0670       edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
0671     }
0672   }
0673 #endif
0674 }
0675 
0676 void HGCalGeomParameters::loadGeometryHexagon8(const DDFilteredView& _fv, HGCalParameters& php, int firstLayer) {
0677   DDFilteredView fv = _fv;
0678   bool dodet(true);
0679   std::map<int, HGCalGeomParameters::layerParameters> layers;
0680   std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0681   int levelTop = 3 + std::max(php.levelT_[0], php.levelT_[1]);
0682 #ifdef EDM_ML_DEBUG
0683   int ntot(0);
0684 #endif
0685   while (dodet) {
0686 #ifdef EDM_ML_DEBUG
0687     ++ntot;
0688 #endif
0689     std::vector<int> copy = fv.copyNumbers();
0690     int nsiz = static_cast<int>(copy.size());
0691     if (nsiz < levelTop) {
0692       int lay = copy[nsiz - 1];
0693       int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
0694       if (zside != 1)
0695         zside = -1;
0696       const DDSolid& sol = fv.logicalPart().solid();
0697 #ifdef EDM_ML_DEBUG
0698       edm::LogVerbatim("HGCalGeom") << sol.name() << " shape " << sol.shape() << " size " << nsiz << ":" << levelTop
0699                                     << " lay " << lay << " z " << zside;
0700 #endif
0701       if (lay == 0) {
0702         throw cms::Exception("DDException")
0703             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0704       } else if (sol.shape() == DDSolidShape::ddtubs) {
0705         if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0706           php.layer_.emplace_back(lay);
0707         const DDTubs& tube = static_cast<DDTubs>(sol);
0708         double rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
0709         double rout = HGCalParameters::k_ScaleFromDDD * tube.rOut();
0710         auto itr = layers.find(lay);
0711         if (itr == layers.end()) {
0712           double zp = HGCalParameters::k_ScaleFromDDD * fv.translation().Z();
0713           HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
0714           layers[lay] = laypar;
0715         } else {
0716           (itr->second).rmin = std::min(rin, (itr->second).rmin);
0717           (itr->second).rmax = std::max(rout, (itr->second).rmax);
0718         }
0719         if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
0720           DD3Vector x, y, z;
0721           fv.rotation().GetComponents(x, y, z);
0722           const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0723           const CLHEP::HepRotation hr(rotation);
0724           double xx =
0725               ((std::abs(fv.translation().X()) < tolerance) ? 0
0726                                                             : HGCalParameters::k_ScaleFromDDD * fv.translation().X());
0727           double yy =
0728               ((std::abs(fv.translation().Y()) < tolerance) ? 0
0729                                                             : HGCalParameters::k_ScaleFromDDD * fv.translation().Y());
0730           const CLHEP::Hep3Vector h3v(xx, yy, HGCalParameters::k_ScaleFromDDD * fv.translation().Z());
0731           HGCalParameters::hgtrform mytrf;
0732           mytrf.zp = zside;
0733           mytrf.lay = lay;
0734           mytrf.sec = 0;
0735           mytrf.subsec = 0;
0736           mytrf.h3v = h3v;
0737           mytrf.hr = hr;
0738           trforms[std::make_pair(lay, zside)] = mytrf;
0739         }
0740       }
0741     }
0742     dodet = fv.next();
0743   }
0744 #ifdef EDM_ML_DEBUG
0745   edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot;
0746 #endif
0747   loadGeometryHexagon8(layers, trforms, firstLayer, php);
0748 }
0749 
0750 void HGCalGeomParameters::loadGeometryHexagon8(const cms::DDCompactView* cpv,
0751                                                HGCalParameters& php,
0752                                                const std::string& sdTag1,
0753                                                int firstLayer) {
0754   const cms::DDFilter filter("Volume", sdTag1);
0755   cms::DDFilteredView fv((*cpv), filter);
0756   std::map<int, HGCalGeomParameters::layerParameters> layers;
0757   std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0758   int levelTop = 3 + std::max(php.levelT_[0], php.levelT_[1]);
0759 #ifdef EDM_ML_DEBUG
0760   int ntot(0);
0761 #endif
0762   while (fv.firstChild()) {
0763 #ifdef EDM_ML_DEBUG
0764     ++ntot;
0765 #endif
0766     // Layers first
0767     int nsiz = static_cast<int>(fv.level());
0768     if (nsiz < levelTop) {
0769       std::vector<int> copy = fv.copyNos();
0770       int lay = copy[0];
0771       int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
0772       if (zside != 1)
0773         zside = -1;
0774 #ifdef EDM_ML_DEBUG
0775       edm::LogVerbatim("HGCalGeom") << fv.name() << " shape " << cms::dd::name(cms::DDSolidShapeMap, fv.shape())
0776                                     << " size " << nsiz << ":" << levelTop << " lay " << lay << " z " << zside << ":"
0777                                     << php.levelZSide_;
0778 #endif
0779       if (lay == 0) {
0780         throw cms::Exception("DDException")
0781             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0782       } else if (fv.shape() == cms::DDSolidShape::ddtubs) {
0783         if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0784           php.layer_.emplace_back(lay);
0785         const std::vector<double>& pars = fv.parameters();
0786         double rin = HGCalParameters::k_ScaleFromDD4hep * pars[0];
0787         double rout = HGCalParameters::k_ScaleFromDD4hep * pars[1];
0788         auto itr = layers.find(lay);
0789         if (itr == layers.end()) {
0790           double zp = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z();
0791           HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
0792           layers[lay] = laypar;
0793         } else {
0794           (itr->second).rmin = std::min(rin, (itr->second).rmin);
0795           (itr->second).rmax = std::max(rout, (itr->second).rmax);
0796         }
0797         if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
0798           DD3Vector x, y, z;
0799           fv.rotation().GetComponents(x, y, z);
0800           const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0801           const CLHEP::HepRotation hr(rotation);
0802           double xx = ((std::abs(fv.translation().X()) < tolerance)
0803                            ? 0
0804                            : HGCalParameters::k_ScaleFromDD4hep * fv.translation().X());
0805           double yy = ((std::abs(fv.translation().Y()) < tolerance)
0806                            ? 0
0807                            : HGCalParameters::k_ScaleFromDD4hep * fv.translation().Y());
0808           const CLHEP::Hep3Vector h3v(xx, yy, HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z());
0809           HGCalParameters::hgtrform mytrf;
0810           mytrf.zp = zside;
0811           mytrf.lay = lay;
0812           mytrf.sec = 0;
0813           mytrf.subsec = 0;
0814           mytrf.h3v = h3v;
0815           mytrf.hr = hr;
0816           trforms[std::make_pair(lay, zside)] = mytrf;
0817         }
0818       }
0819     }
0820   }
0821 #ifdef EDM_ML_DEBUG
0822   edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot;
0823 #endif
0824   loadGeometryHexagon8(layers, trforms, firstLayer, php);
0825 }
0826 
0827 void HGCalGeomParameters::loadGeometryHexagonModule(const DDCompactView* cpv,
0828                                                     HGCalParameters& php,
0829                                                     const std::string& sdTag1,
0830                                                     const std::string& sdTag2,
0831                                                     int firstLayer) {
0832 #ifdef EDM_ML_DEBUG
0833   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters (DDD)::loadGeometryHexagonModule called with tags " << sdTag1
0834                                 << ":" << sdTag2 << " firstLayer " << firstLayer << ":" << php.firstMixedLayer_;
0835   int ntot1(0), ntot2(0);
0836 #endif
0837   std::map<int, HGCalGeomParameters::layerParameters> layers;
0838   std::map<std::pair<int, int>, double> zvals;
0839   std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0840   int levelTop = php.levelT_[0];
0841 
0842   std::string attribute = "Volume";
0843   DDValue val1(attribute, sdTag2, 0.0);
0844   DDSpecificsMatchesValueFilter filter1{val1};
0845   DDFilteredView fv1(*cpv, filter1);
0846   bool dodet = fv1.firstChild();
0847   while (dodet) {
0848 #ifdef EDM_ML_DEBUG
0849     ++ntot1;
0850 #endif
0851     std::vector<int> copy = fv1.copyNumbers();
0852     int nsiz = static_cast<int>(copy.size());
0853     if (levelTop < nsiz) {
0854       int lay = copy[levelTop];
0855       int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
0856       if (zside != 1)
0857         zside = -1;
0858       if (lay == 0) {
0859         throw cms::Exception("DDException")
0860             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0861       } else {
0862         if (zvals.find(std::make_pair(lay, zside)) == zvals.end()) {
0863           zvals[std::make_pair(lay, zside)] = HGCalParameters::k_ScaleFromDDD * fv1.translation().Z();
0864 #ifdef EDM_ML_DEBUG
0865           std::ostringstream st1;
0866           st1 << "Name0 " << fv1.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
0867               << nsiz;
0868           for (const auto& c : copy)
0869             st1 << ":" << c;
0870           st1 << " Z " << zvals[std::make_pair(lay, zside)];
0871           edm::LogVerbatim("HGCalGeom") << st1.str();
0872 #endif
0873         }
0874       }
0875     }
0876     dodet = fv1.next();
0877   }
0878 
0879   DDValue val2(attribute, sdTag1, 0.0);
0880   DDSpecificsMatchesValueFilter filter2{val2};
0881   DDFilteredView fv2(*cpv, filter2);
0882   dodet = fv2.firstChild();
0883   while (dodet) {
0884     std::vector<int> copy = fv2.copyNumbers();
0885     int nsiz = static_cast<int>(copy.size());
0886 #ifdef EDM_ML_DEBUG
0887     ++ntot2;
0888     edm::LogVerbatim("HGCalGeom") << "loadGeometryHexagonModule:: nsiz " << nsiz << " Ltop " << levelTop;
0889 #endif
0890     if (levelTop < nsiz) {
0891       int lay = copy[levelTop];
0892       int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
0893       if (zside != 1)
0894         zside = -1;
0895       const DDSolid& sol = fv2.logicalPart().solid();
0896 #ifdef EDM_ML_DEBUG
0897       std::ostringstream st2;
0898       st2 << "Name1 " << sol.name() << " shape " << sol.shape() << " LTop " << levelTop << ":" << lay << " ZSide "
0899           << zside << ":" << php.levelZSide_ << " # of levels " << nsiz;
0900       for (const auto& c : copy)
0901         st2 << ":" << c;
0902       edm::LogVerbatim("HGCalGeom") << st2.str();
0903 #endif
0904       if (lay == 0) {
0905         throw cms::Exception("DDException")
0906             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0907       } else if ((sol.shape() == DDSolidShape::ddtubs) || (sol.shape() == DDSolidShape::ddbox)) {
0908         if (zvals.find(std::make_pair(lay, zside)) != zvals.end()) {
0909           if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0910             php.layer_.emplace_back(lay);
0911           auto itr = layers.find(lay);
0912           if (itr == layers.end()) {
0913             double rin(0), rout(0);
0914             if (sol.shape() == DDSolidShape::ddtubs) {
0915               const DDTubs& tube = static_cast<DDTubs>(sol);
0916               rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
0917               rout = (php.firstMixedLayer_ > 0 && lay >= php.firstMixedLayer_)
0918                          ? php.radiusMixBoundary_[lay - php.firstMixedLayer_]
0919                          : HGCalParameters::k_ScaleFromDDD * tube.rOut();
0920             } else {
0921               const DDBox& box = static_cast<DDBox>(sol);
0922               rout = HGCalParameters::k_ScaleFromDDD * box.halfX();
0923             }
0924             double zp = zvals[std::make_pair(lay, 1)];
0925             HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
0926             layers[lay] = laypar;
0927 #ifdef EDM_ML_DEBUG
0928             std::ostringstream st3;
0929             st3 << "Name1 " << fv2.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
0930                 << nsiz;
0931             for (const auto& c : copy)
0932               st3 << ":" << c;
0933             st3 << " R " << rin << ":" << rout;
0934             edm::LogVerbatim("HGCalGeom") << st3.str();
0935 #endif
0936           }
0937 
0938           if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
0939             DD3Vector x, y, z;
0940             fv2.rotation().GetComponents(x, y, z);
0941             const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0942             const CLHEP::HepRotation hr(rotation);
0943             double xx = ((std::abs(fv2.translation().X()) < tolerance)
0944                              ? 0
0945                              : HGCalParameters::k_ScaleFromDDD * fv2.translation().X());
0946             double yy = ((std::abs(fv2.translation().Y()) < tolerance)
0947                              ? 0
0948                              : HGCalParameters::k_ScaleFromDDD * fv2.translation().Y());
0949             const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
0950             HGCalParameters::hgtrform mytrf;
0951             mytrf.zp = zside;
0952             mytrf.lay = lay;
0953             mytrf.sec = 0;
0954             mytrf.subsec = 0;
0955             mytrf.h3v = h3v;
0956             mytrf.hr = hr;
0957             trforms[std::make_pair(lay, zside)] = mytrf;
0958 #ifdef EDM_ML_DEBUG
0959             edm::LogVerbatim("HGCalGeom") << "Translation " << h3v;
0960 #endif
0961           }
0962         }
0963       }
0964     }
0965     dodet = fv2.next();
0966   }
0967 #ifdef EDM_ML_DEBUG
0968   edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot1 << ":" << ntot2;
0969 #endif
0970   loadGeometryHexagon8(layers, trforms, firstLayer, php);
0971 }
0972 
0973 void HGCalGeomParameters::loadGeometryHexagonModule(const cms::DDCompactView* cpv,
0974                                                     HGCalParameters& php,
0975                                                     const std::string& sdTag1,
0976                                                     const std::string& sdTag2,
0977                                                     int firstLayer) {
0978 #ifdef EDM_ML_DEBUG
0979   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters (DD4hep)::loadGeometryHexagonModule called with tags " << sdTag1
0980                                 << ":" << sdTag2 << " firstLayer " << firstLayer;
0981   int ntot1(0), ntot2(0);
0982 #endif
0983   std::map<int, HGCalGeomParameters::layerParameters> layers;
0984   std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0985   std::map<std::pair<int, int>, double> zvals;
0986   int levelTop = php.levelT_[0];
0987 
0988   const cms::DDFilter filter1("Volume", sdTag2);
0989   cms::DDFilteredView fv1((*cpv), filter1);
0990   while (fv1.firstChild()) {
0991 #ifdef EDM_ML_DEBUG
0992     ++ntot1;
0993 #endif
0994     int nsiz = static_cast<int>(fv1.level());
0995     if (nsiz > levelTop) {
0996       std::vector<int> copy = fv1.copyNos();
0997       int lay = copy[nsiz - levelTop - 1];
0998       int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
0999       if (zside != 1)
1000         zside = -1;
1001       if (lay == 0) {
1002         throw cms::Exception("DDException")
1003             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
1004       } else {
1005         if (zvals.find(std::make_pair(lay, zside)) == zvals.end()) {
1006           zvals[std::make_pair(lay, zside)] = HGCalParameters::k_ScaleFromDD4hep * fv1.translation().Z();
1007 #ifdef EDM_ML_DEBUG
1008           std::ostringstream st1;
1009           st1 << "Name0 " << fv1.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
1010               << nsiz;
1011           for (const auto& c : copy)
1012             st1 << ":" << c;
1013           st1 << " Z " << zvals[std::make_pair(lay, zside)];
1014           edm::LogVerbatim("HGCalGeom") << st1.str();
1015 #endif
1016         }
1017       }
1018     }
1019   }
1020 
1021   const cms::DDFilter filter2("Volume", sdTag1);
1022   cms::DDFilteredView fv2((*cpv), filter2);
1023   while (fv2.firstChild()) {
1024     // Layers first
1025     int nsiz = static_cast<int>(fv2.level());
1026 #ifdef EDM_ML_DEBUG
1027     ++ntot2;
1028     edm::LogVerbatim("HGCalGeom") << "loadGeometryHexagonModule:: nsiz " << nsiz << " Ltop " << levelTop;
1029 #endif
1030     if (nsiz > levelTop) {
1031       std::vector<int> copy = fv2.copyNos();
1032       int lay = copy[nsiz - levelTop - 1];
1033       int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
1034       if (zside != 1)
1035         zside = -1;
1036 #ifdef EDM_ML_DEBUG
1037       std::ostringstream st2;
1038       st2 << "Name1 " << fv2.name() << "Shape " << cms::dd::name(cms::DDSolidShapeMap, fv2.shape()) << " LTop "
1039           << levelTop << ":" << lay << " ZSide " << zside << ":" << php.levelZSide_ << " # of levels " << nsiz;
1040       for (const auto& c : copy)
1041         st2 << ":" << c;
1042       edm::LogVerbatim("HGCalGeom") << st2.str();
1043 #endif
1044       if (lay == 0) {
1045         throw cms::Exception("DDException")
1046             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
1047       } else {
1048         if (zvals.find(std::make_pair(lay, zside)) != zvals.end()) {
1049           if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
1050             php.layer_.emplace_back(lay);
1051           auto itr = layers.find(lay);
1052           if (itr == layers.end()) {
1053             const std::vector<double>& pars = fv2.parameters();
1054             double rin(0), rout(0);
1055             if (dd4hep::isA<dd4hep::Box>(fv2.solid())) {
1056               rout = HGCalParameters::k_ScaleFromDD4hep * pars[0];
1057             } else {
1058               rin = HGCalParameters::k_ScaleFromDD4hep * pars[0];
1059               rout = (php.firstMixedLayer_ > 0 && lay >= php.firstMixedLayer_)
1060                          ? php.radiusMixBoundary_[lay - php.firstMixedLayer_]
1061                          : HGCalParameters::k_ScaleFromDD4hep * pars[1];
1062             }
1063             double zp = zvals[std::make_pair(lay, 1)];
1064             HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
1065             layers[lay] = laypar;
1066 #ifdef EDM_ML_DEBUG
1067             std::ostringstream st3;
1068             st3 << "Name2 " << fv2.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
1069                 << nsiz;
1070             for (const auto& c : copy)
1071               st3 << ":" << c;
1072             st3 << " R " << rin << ":" << rout;
1073             edm::LogVerbatim("HGCalGeom") << st3.str();
1074 #endif
1075           }
1076 
1077           if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
1078             DD3Vector x, y, z;
1079             fv2.rotation().GetComponents(x, y, z);
1080             const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
1081             const CLHEP::HepRotation hr(rotation);
1082             double xx = ((std::abs(fv2.translation().X()) < tolerance)
1083                              ? 0
1084                              : HGCalParameters::k_ScaleFromDD4hep * fv2.translation().X());
1085             double yy = ((std::abs(fv2.translation().Y()) < tolerance)
1086                              ? 0
1087                              : HGCalParameters::k_ScaleFromDD4hep * fv2.translation().Y());
1088             const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
1089             HGCalParameters::hgtrform mytrf;
1090             mytrf.zp = zside;
1091             mytrf.lay = lay;
1092             mytrf.sec = 0;
1093             mytrf.subsec = 0;
1094             mytrf.h3v = h3v;
1095             mytrf.hr = hr;
1096             trforms[std::make_pair(lay, zside)] = mytrf;
1097 #ifdef EDM_ML_DEBUG
1098             edm::LogVerbatim("HGCalGeom") << "Translation " << h3v;
1099 #endif
1100           }
1101         }
1102       }
1103     }
1104   }
1105 #ifdef EDM_ML_DEBUG
1106   edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot1 << ":" << ntot2;
1107 #endif
1108   loadGeometryHexagon8(layers, trforms, firstLayer, php);
1109 }
1110 
1111 void HGCalGeomParameters::loadGeometryHexagon8(const std::map<int, HGCalGeomParameters::layerParameters>& layers,
1112                                                std::map<std::pair<int, int>, HGCalParameters::hgtrform>& trforms,
1113                                                const int& firstLayer,
1114                                                HGCalParameters& php) {
1115   double rmin(0), rmax(0);
1116   for (unsigned int i = 0; i < layers.size(); ++i) {
1117     for (auto& layer : layers) {
1118       if (layer.first == static_cast<int>(i + firstLayer)) {
1119         php.layerIndex_.emplace_back(i);
1120         php.rMinLayHex_.emplace_back(layer.second.rmin);
1121         php.rMaxLayHex_.emplace_back(layer.second.rmax);
1122         php.zLayerHex_.emplace_back(layer.second.zpos);
1123         if (i == 0) {
1124           rmin = layer.second.rmin;
1125           rmax = layer.second.rmax;
1126         } else {
1127           if (rmin > layer.second.rmin)
1128             rmin = layer.second.rmin;
1129           if (rmax < layer.second.rmax)
1130             rmax = layer.second.rmax;
1131         }
1132         break;
1133       }
1134     }
1135   }
1136   php.rLimit_.emplace_back(rmin);
1137   php.rLimit_.emplace_back(rmax);
1138   php.depth_ = php.layer_;
1139   php.depthIndex_ = php.layerIndex_;
1140   php.depthLayerF_ = php.layerIndex_;
1141 
1142   for (unsigned int i = 0; i < php.layer_.size(); ++i) {
1143     for (auto& trform : trforms) {
1144       if (trform.first.first == static_cast<int>(i + firstLayer)) {
1145         php.fillTrForm(trform.second);
1146       }
1147     }
1148   }
1149 
1150 #ifdef EDM_ML_DEBUG
1151   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Minimum/maximum R " << php.rLimit_[0] << ":" << php.rLimit_[1];
1152   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.zLayerHex_.size() << " layers";
1153   for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
1154     int k = php.layerIndex_[i];
1155     edm::LogVerbatim("HGCalGeom") << "Layer[" << i << ":" << k << ":" << php.layer_[k]
1156                                   << "] with r = " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
1157                                   << " at z = " << php.zLayerHex_[i];
1158   }
1159   edm::LogVerbatim("HGCalGeom") << "Obtained " << php.trformIndex_.size() << " transformation matrices";
1160   for (unsigned int k = 0; k < php.trformIndex_.size(); ++k) {
1161     edm::LogVerbatim("HGCalGeom") << "Matrix[" << k << "] (" << std::hex << php.trformIndex_[k] << std::dec
1162                                   << ") Translation (" << php.trformTranX_[k] << ", " << php.trformTranY_[k] << ", "
1163                                   << php.trformTranZ_[k] << " Rotation (" << php.trformRotXX_[k] << ", "
1164                                   << php.trformRotYX_[k] << ", " << php.trformRotZX_[k] << ", " << php.trformRotXY_[k]
1165                                   << ", " << php.trformRotYY_[k] << ", " << php.trformRotZY_[k] << ", "
1166                                   << php.trformRotXZ_[k] << ", " << php.trformRotYZ_[k] << ", " << php.trformRotZZ_[k]
1167                                   << ")";
1168   }
1169 #endif
1170 }
1171 
1172 void HGCalGeomParameters::loadSpecParsHexagon(const DDFilteredView& fv,
1173                                               HGCalParameters& php,
1174                                               const DDCompactView* cpv,
1175                                               const std::string& sdTag1,
1176                                               const std::string& sdTag2) {
1177   DDsvalues_type sv(fv.mergedSpecifics());
1178   php.boundR_ = getDDDArray("RadiusBound", sv, 4);
1179   rescale(php.boundR_, HGCalParameters::k_ScaleFromDDD);
1180   php.rLimit_ = getDDDArray("RadiusLimits", sv, 2);
1181   rescale(php.rLimit_, HGCalParameters::k_ScaleFromDDD);
1182   php.levelT_ = dbl_to_int(getDDDArray("LevelTop", sv, 0));
1183 
1184   // Grouping of layers
1185   php.layerGroup_ = dbl_to_int(getDDDArray("GroupingZFine", sv, 0));
1186   php.layerGroupM_ = dbl_to_int(getDDDArray("GroupingZMid", sv, 0));
1187   php.layerGroupO_ = dbl_to_int(getDDDArray("GroupingZOut", sv, 0));
1188   php.slopeMin_ = getDDDArray("Slope", sv, 1);
1189   const auto& dummy2 = getDDDArray("LayerOffset", sv, 0);
1190   if (!dummy2.empty())
1191     php.layerOffset_ = dummy2[0];
1192   else
1193     php.layerOffset_ = 0;
1194 
1195   // Wafer size
1196   std::string attribute = "Volume";
1197   DDSpecificsMatchesValueFilter filter1{DDValue(attribute, sdTag1, 0.0)};
1198   DDFilteredView fv1(*cpv, filter1);
1199   if (fv1.firstChild()) {
1200     DDsvalues_type sv(fv1.mergedSpecifics());
1201     const auto& dummy = getDDDArray("WaferSize", sv, 0);
1202     waferSize_ = dummy[0];
1203   }
1204 
1205   // Cell size
1206   DDSpecificsMatchesValueFilter filter2{DDValue(attribute, sdTag2, 0.0)};
1207   DDFilteredView fv2(*cpv, filter2);
1208   if (fv2.firstChild()) {
1209     DDsvalues_type sv(fv2.mergedSpecifics());
1210     php.cellSize_ = getDDDArray("CellSize", sv, 0);
1211   }
1212 
1213   loadSpecParsHexagon(php);
1214 }
1215 
1216 void HGCalGeomParameters::loadSpecParsHexagon(const cms::DDFilteredView& fv,
1217                                               HGCalParameters& php,
1218                                               const std::string& sdTag1,
1219                                               const std::string& sdTag2,
1220                                               const std::string& sdTag3,
1221                                               const std::string& sdTag4) {
1222   php.boundR_ = fv.get<std::vector<double> >(sdTag4, "RadiusBound");
1223   rescale(php.boundR_, HGCalParameters::k_ScaleFromDD4hep);
1224   php.rLimit_ = fv.get<std::vector<double> >(sdTag4, "RadiusLimits");
1225   rescale(php.rLimit_, HGCalParameters::k_ScaleFromDD4hep);
1226   php.levelT_ = dbl_to_int(fv.get<std::vector<double> >(sdTag4, "LevelTop"));
1227 
1228   // Grouping of layers
1229   php.layerGroup_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZFine"));
1230   php.layerGroupM_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZMid"));
1231   php.layerGroupO_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZOut"));
1232   php.slopeMin_ = fv.get<std::vector<double> >(sdTag4, "Slope");
1233   if (php.slopeMin_.empty())
1234     php.slopeMin_.emplace_back(0);
1235 
1236   // Wafer size
1237   const auto& dummy = fv.get<std::vector<double> >(sdTag2, "WaferSize");
1238   waferSize_ = dummy[0] * HGCalParameters::k_ScaleFromDD4hepToG4;
1239 
1240   // Cell size
1241   php.cellSize_ = fv.get<std::vector<double> >(sdTag3, "CellSize");
1242   rescale(php.cellSize_, HGCalParameters::k_ScaleFromDD4hepToG4);
1243 
1244   // Layer Offset
1245   const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1246   if (!dummy2.empty()) {
1247     php.layerOffset_ = dummy2[0];
1248   } else {
1249     php.layerOffset_ = 0;
1250   }
1251 
1252   loadSpecParsHexagon(php);
1253 }
1254 
1255 void HGCalGeomParameters::loadSpecParsHexagon(const HGCalParameters& php) {
1256 #ifdef EDM_ML_DEBUG
1257   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer radius ranges"
1258                                 << " for cell grouping " << php.boundR_[0] << ":" << php.boundR_[1] << ":"
1259                                 << php.boundR_[2] << ":" << php.boundR_[3];
1260   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Minimum/maximum R " << php.rLimit_[0] << ":" << php.rLimit_[1];
1261   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: LevelTop " << php.levelT_[0];
1262   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: minimum slope " << php.slopeMin_[0] << " and layer groupings "
1263                                 << "for the 3 ranges:";
1264   for (unsigned int k = 0; k < php.layerGroup_.size(); ++k)
1265     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerGroup_[k] << ":" << php.layerGroupM_[k] << ":"
1266                                   << php.layerGroupO_[k];
1267   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Wafer Size: " << waferSize_;
1268   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: " << php.cellSize_.size() << " cells of sizes:";
1269   for (unsigned int k = 0; k < php.cellSize_.size(); ++k)
1270     edm::LogVerbatim("HGCalGeom") << " [" << k << "] " << php.cellSize_[k];
1271   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: First Layer " << php.firstLayer_ << " and layer offset "
1272                                 << php.layerOffset_;
1273 #endif
1274 }
1275 
1276 void HGCalGeomParameters::loadSpecParsHexagon8(const DDFilteredView& fv, HGCalParameters& php) {
1277   DDsvalues_type sv(fv.mergedSpecifics());
1278   php.cellThickness_ = getDDDArray("CellThickness", sv, 3);
1279   rescale(php.cellThickness_, HGCalParameters::k_ScaleFromDDD);
1280   if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1281       (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1282     php.waferThickness_ = getDDDArray("WaferThickness", sv, 3);
1283     rescale(php.waferThickness_, HGCalParameters::k_ScaleFromDDD);
1284   } else {
1285     for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1286       php.waferThickness_.emplace_back(php.waferThick_);
1287   }
1288 
1289   php.radius100to200_ = getDDDArray("Radius100to200", sv, 5);
1290   php.radius200to300_ = getDDDArray("Radius200to300", sv, 5);
1291 
1292   const auto& dummy = getDDDArray("RadiusCuts", sv, 4);
1293   php.choiceType_ = static_cast<int>(dummy[0]);
1294   php.nCornerCut_ = static_cast<int>(dummy[1]);
1295   php.fracAreaMin_ = dummy[2];
1296   php.zMinForRad_ = HGCalParameters::k_ScaleFromDDD * dummy[3];
1297 
1298   php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1299   rescale(php.radiusMixBoundary_, HGCalParameters::k_ScaleFromDDD);
1300 
1301   php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1302   php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1303   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDDD);
1304   php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1305   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDDD);
1306 
1307   php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1308   php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1309   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDDD);
1310   php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1311   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDDD);
1312 
1313   php.zRanges_ = fv.vector("ZRanges");
1314   rescale(php.zRanges_, HGCalParameters::k_ScaleFromDDD);
1315 
1316   const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1317   php.layerOffset_ = dummy2[0];
1318   php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1319 
1320   if (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell) {
1321     const auto& dummy3 = fv.vector("CalibCellRadius");
1322     php.calibCellRHD_ = HGCalParameters::k_ScaleFromDDD * dummy3[0];
1323     php.calibCellFullHD_ = dbl_to_int(fv.vector("CalibCellFullHD"));
1324     php.calibCellPartHD_ = dbl_to_int(fv.vector("CalibCellPartHD"));
1325     php.calibCellRLD_ = HGCalParameters::k_ScaleFromDDD * dummy3[1];
1326     php.calibCellFullLD_ = dbl_to_int(fv.vector("CalibCellFullLD"));
1327     php.calibCellPartLD_ = dbl_to_int(fv.vector("CalibCellPartLD"));
1328   } else {
1329     php.calibCellRHD_ = php.calibCellRLD_ = 0;
1330     php.calibCellFullHD_.clear();
1331     php.calibCellPartHD_.clear();
1332     php.calibCellFullLD_.clear();
1333     php.calibCellPartLD_.clear();
1334   }
1335 
1336   loadSpecParsHexagon8(php);
1337 
1338   // Read in parameters from Philip's file
1339   if (php.waferMaskMode_ > 1) {
1340     std::vector<int> layerType, waferIndex, waferProperties;
1341     std::vector<double> cassetteShift;
1342     if (php.waferMaskMode_ == siliconFileEE) {
1343       waferIndex = dbl_to_int(fv.vector("WaferIndexEE"));
1344       waferProperties = dbl_to_int(fv.vector("WaferPropertiesEE"));
1345     } else if (php.waferMaskMode_ == siliconCassetteEE) {
1346       waferIndex = dbl_to_int(fv.vector("WaferIndexEE"));
1347       waferProperties = dbl_to_int(fv.vector("WaferPropertiesEE"));
1348       cassetteShift = fv.vector("CassetteShiftEE");
1349     } else if (php.waferMaskMode_ == siliconFileHE) {
1350       waferIndex = dbl_to_int(fv.vector("WaferIndexHE"));
1351       waferProperties = dbl_to_int(fv.vector("WaferPropertiesHE"));
1352     } else if (php.waferMaskMode_ == siliconCassetteHE) {
1353       waferIndex = dbl_to_int(fv.vector("WaferIndexHE"));
1354       waferProperties = dbl_to_int(fv.vector("WaferPropertiesHE"));
1355       cassetteShift = fv.vector("CassetteShiftHE");
1356     }
1357     if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1358         (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1359       if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1360         layerType = dbl_to_int(fv.vector("LayerTypesEE"));
1361       } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1362         layerType = dbl_to_int(fv.vector("LayerTypesHE"));
1363       }
1364     }
1365 
1366     php.cassetteShift_ = cassetteShift;
1367     rescale(php.cassetteShift_, HGCalParameters::k_ScaleFromDDD);
1368     loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1369   }
1370 }
1371 
1372 void HGCalGeomParameters::loadSpecParsHexagon8(const cms::DDFilteredView& fv,
1373                                                const cms::DDVectorsMap& vmap,
1374                                                HGCalParameters& php,
1375                                                const std::string& sdTag1) {
1376   php.cellThickness_ = fv.get<std::vector<double> >(sdTag1, "CellThickness");
1377   rescale(php.cellThickness_, HGCalParameters::k_ScaleFromDD4hep);
1378   if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1379       (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1380     php.waferThickness_ = fv.get<std::vector<double> >(sdTag1, "WaferThickness");
1381     rescale(php.waferThickness_, HGCalParameters::k_ScaleFromDD4hep);
1382   } else {
1383     for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1384       php.waferThickness_.emplace_back(php.waferThick_);
1385   }
1386 
1387   php.radius100to200_ = fv.get<std::vector<double> >(sdTag1, "Radius100to200");
1388   php.radius200to300_ = fv.get<std::vector<double> >(sdTag1, "Radius200to300");
1389 
1390   const auto& dummy = fv.get<std::vector<double> >(sdTag1, "RadiusCuts");
1391   if (dummy.size() > 3) {
1392     php.choiceType_ = static_cast<int>(dummy[0]);
1393     php.nCornerCut_ = static_cast<int>(dummy[1]);
1394     php.fracAreaMin_ = dummy[2];
1395     php.zMinForRad_ = HGCalParameters::k_ScaleFromDD4hep * dummy[3];
1396   } else {
1397     php.choiceType_ = php.nCornerCut_ = php.fracAreaMin_ = php.zMinForRad_ = 0;
1398   }
1399 
1400   php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1401   php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1402   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDD4hep);
1403   php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1404   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDD4hep);
1405 
1406   php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1407   php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1408   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDD4hep);
1409   php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1410   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDD4hep);
1411   unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1412   for (unsigned int k = 0; k < kmax; ++k)
1413     php.slopeTop_.emplace_back(0);
1414 
1415   const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1416   if (!dummy2.empty()) {
1417     php.layerOffset_ = dummy2[0];
1418   } else {
1419     php.layerOffset_ = 0;
1420   }
1421 
1422   php.calibCellRHD_ = php.calibCellRLD_ = 0;
1423   php.calibCellFullHD_.clear();
1424   php.calibCellPartHD_.clear();
1425   php.calibCellFullLD_.clear();
1426   php.calibCellPartLD_.clear();
1427   for (auto const& it : vmap) {
1428     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1429       for (const auto& i : it.second)
1430         php.radiusMixBoundary_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1431     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1432       for (const auto& i : it.second)
1433         php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1434     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1435       for (const auto& i : it.second)
1436         php.layerCenter_.emplace_back(std::round(i));
1437     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellRadius")) {
1438       php.calibCellRHD_ = HGCalParameters::k_ScaleFromDD4hep * it.second[0];
1439       php.calibCellRLD_ = HGCalParameters::k_ScaleFromDDD * it.second[1];
1440     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellFullHD")) {
1441       for (const auto& i : it.second)
1442         php.calibCellFullHD_.emplace_back(std::round(i));
1443     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellPartHD")) {
1444       for (const auto& i : it.second)
1445         php.calibCellPartHD_.emplace_back(std::round(i));
1446     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellFullLD")) {
1447       for (const auto& i : it.second)
1448         php.calibCellFullLD_.emplace_back(std::round(i));
1449     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellPartLD")) {
1450       for (const auto& i : it.second)
1451         php.calibCellPartLD_.emplace_back(std::round(i));
1452     }
1453   }
1454 
1455   loadSpecParsHexagon8(php);
1456 
1457   // Read in parameters from Philip's file
1458   if (php.waferMaskMode_ > 1) {
1459     std::vector<int> layerType, waferIndex, waferProperties;
1460     std::vector<double> cassetteShift;
1461     if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1462       for (auto const& it : vmap) {
1463         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferIndexEE")) {
1464           for (const auto& i : it.second)
1465             waferIndex.emplace_back(std::round(i));
1466         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferPropertiesEE")) {
1467           for (const auto& i : it.second)
1468             waferProperties.emplace_back(std::round(i));
1469         }
1470       }
1471       if (php.waferMaskMode_ == siliconCassetteEE) {
1472         for (auto const& it : vmap) {
1473           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftEE")) {
1474             for (const auto& i : it.second)
1475               cassetteShift.emplace_back(i);
1476           }
1477         }
1478       }
1479     } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1480       for (auto const& it : vmap) {
1481         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferIndexHE")) {
1482           for (const auto& i : it.second)
1483             waferIndex.emplace_back(std::round(i));
1484         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferPropertiesHE")) {
1485           for (const auto& i : it.second)
1486             waferProperties.emplace_back(std::round(i));
1487         }
1488       }
1489       if (php.waferMaskMode_ == siliconCassetteHE) {
1490         for (auto const& it : vmap) {
1491           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1492             for (const auto& i : it.second)
1493               cassetteShift.emplace_back(i);
1494           }
1495         }
1496       }
1497     }
1498     if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1499         (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1500       if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1501         for (auto const& it : vmap) {
1502           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerTypesEE")) {
1503             for (const auto& i : it.second)
1504               layerType.emplace_back(std::round(i));
1505           }
1506         }
1507       } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1508         for (auto const& it : vmap) {
1509           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerTypesHE")) {
1510             for (const auto& i : it.second)
1511               layerType.emplace_back(std::round(i));
1512           }
1513         }
1514       }
1515     }
1516 
1517     php.cassetteShift_ = cassetteShift;
1518     rescale(php.cassetteShift_, HGCalParameters::k_ScaleFromDD4hep);
1519     loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1520   }
1521 }
1522 
1523 void HGCalGeomParameters::loadSpecParsHexagon8(HGCalParameters& php) {
1524 #ifdef EDM_ML_DEBUG
1525   for (unsigned int k = 0; k < php.waferThickness_.size(); ++k)
1526     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer[" << k << "] Thickness " << php.waferThickness_[k];
1527   for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1528     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: cell[" << k << "] Thickness " << php.cellThickness_[k];
1529   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Polynomial "
1530                                 << "parameters for 120 to 200 micron "
1531                                 << "transition with " << php.radius100to200_.size() << " elements";
1532   for (unsigned int k = 0; k < php.radius100to200_.size(); ++k)
1533     edm::LogVerbatim("HGCalGeom") << "Element [" << k << "] " << php.radius100to200_[k];
1534   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Polynomial "
1535                                 << "parameters for 200 to 300 micron "
1536                                 << "transition with " << php.radius200to300_.size() << " elements";
1537   for (unsigned int k = 0; k < php.radius200to300_.size(); ++k)
1538     edm::LogVerbatim("HGCalGeom") << "Element [" << k << "] " << php.radius200to300_[k];
1539   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Parameters for the"
1540                                 << " transition " << php.choiceType_ << ":" << php.nCornerCut_ << ":"
1541                                 << php.fracAreaMin_ << ":" << php.zMinForRad_;
1542   for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k)
1543     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k];
1544   for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1545     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1546                                   << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1547   for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1548     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1549                                   << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1550   std::ostringstream st1;
1551   for (unsigned int k = 0; k < php.zRanges_.size(); ++k)
1552     st1 << ":" << php.zRanges_[k];
1553   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary[" << php.zRanges_.size() << "] " << st1.str();
1554   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1555                                 << php.layerCenter_.size();
1556   for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1557     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1558   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Calibration cells in HD having radius " << php.calibCellRHD_
1559                                 << " in arrays of size " << php.calibCellFullHD_.size() << ":"
1560                                 << php.calibCellPartHD_.size();
1561   for (unsigned int k = 0; k < php.calibCellFullHD_.size(); ++k)
1562     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.calibCellFullHD_[k] << ":" << php.calibCellPartHD_[k];
1563   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Calibration cells in LD having radius " << php.calibCellRLD_
1564                                 << " in arrays of size " << php.calibCellFullLD_.size() << ":"
1565                                 << php.calibCellPartLD_.size();
1566   for (unsigned int k = 0; k < php.calibCellFullLD_.size(); ++k)
1567     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.calibCellFullLD_[k] << ":" << php.calibCellPartLD_[k];
1568 #endif
1569 }
1570 
1571 void HGCalGeomParameters::loadSpecParsHexagon8(HGCalParameters& php,
1572                                                const std::vector<int>& layerType,
1573                                                const std::vector<int>& waferIndex,
1574                                                const std::vector<int>& waferProperties) {
1575   // Store parameters from Philip's file
1576   for (unsigned int k = 0; k < layerType.size(); ++k) {
1577     php.layerType_.emplace_back(HGCalTypes::layerType(layerType[k]));
1578 #ifdef EDM_ML_DEBUG
1579     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Type " << layerType[k] << ":" << php.layerType_.back();
1580 #endif
1581   }
1582   for (unsigned int k = 0; k < php.layerType_.size(); ++k) {
1583     double cth = (php.layerType_[k] == HGCalTypes::WaferCenterR) ? cos(php.layerRotation_) : 1.0;
1584     double sth = (php.layerType_[k] == HGCalTypes::WaferCenterR) ? sin(php.layerRotation_) : 0.0;
1585     php.layerRotV_.emplace_back(std::make_pair(cth, sth));
1586 #ifdef EDM_ML_DEBUG
1587     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Type " << php.layerType_[k] << " cos|sin(Theta) "
1588                                   << php.layerRotV_.back().first << ":" << php.layerRotV_.back().second;
1589 #endif
1590   }
1591   for (unsigned int k = 0; k < waferIndex.size(); ++k) {
1592     int partial = HGCalProperty::waferPartial(waferProperties[k]);
1593     int orient =
1594         ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) || (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
1595             ? HGCalProperty::waferOrient(waferProperties[k])
1596             : HGCalWaferMask::getRotation(php.waferZSide_, partial, HGCalProperty::waferOrient(waferProperties[k]));
1597     php.waferInfoMap_[waferIndex[k]] = HGCalParameters::waferInfo(HGCalProperty::waferThick(waferProperties[k]),
1598                                                                   partial,
1599                                                                   orient,
1600                                                                   HGCalProperty::waferCassette(waferProperties[k]));
1601 #ifdef EDM_ML_DEBUG
1602     edm::LogVerbatim("HGCalGeom") << "[" << k << ":" << waferIndex[k] << ":"
1603                                   << HGCalWaferIndex::waferLayer(waferIndex[k]) << ":"
1604                                   << HGCalWaferIndex::waferU(waferIndex[k]) << ":"
1605                                   << HGCalWaferIndex::waferV(waferIndex[k]) << "]  Thickness type "
1606                                   << HGCalProperty::waferThick(waferProperties[k]) << " Partial type " << partial
1607                                   << " Orientation " << HGCalProperty::waferOrient(waferProperties[k]) << ":" << orient;
1608 #endif
1609   }
1610 }
1611 
1612 void HGCalGeomParameters::loadSpecParsTrapezoid(const DDFilteredView& fv, HGCalParameters& php) {
1613   DDsvalues_type sv(fv.mergedSpecifics());
1614   php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1615   rescale(php.radiusMixBoundary_, HGCalParameters::k_ScaleFromDDD);
1616 
1617   php.nPhiBinBH_ = dbl_to_int(getDDDArray("NPhiBinBH", sv, 0));
1618   php.layerFrontBH_ = dbl_to_int(getDDDArray("LayerFrontBH", sv, 0));
1619   php.rMinLayerBH_ = getDDDArray("RMinLayerBH", sv, 0);
1620   rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDDD);
1621   php.nCellsFine_ = php.nPhiBinBH_[0];
1622   php.nCellsCoarse_ = php.nPhiBinBH_[1];
1623   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1624   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1625 
1626   php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1627   php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1628   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDDD);
1629   php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1630   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDDD);
1631 
1632   php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1633   php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1634   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDDD);
1635   php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1636   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDDD);
1637 
1638   php.zRanges_ = fv.vector("ZRanges");
1639   rescale(php.zRanges_, HGCalParameters::k_ScaleFromDDD);
1640 
1641   // Offsets
1642   const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1643   php.layerOffset_ = dummy2[0];
1644   php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1645 
1646   loadSpecParsTrapezoid(php);
1647 
1648   // tile parameters from Katja's file
1649   if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette)) {
1650     std::vector<int> tileIndx, tileProperty;
1651     std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4;
1652     std::vector<double> tileRMin, tileRMax;
1653     std::vector<int> tileRingMin, tileRingMax;
1654     std::vector<double> cassetteShift;
1655     tileIndx = dbl_to_int(fv.vector("TileIndex"));
1656     tileProperty = dbl_to_int(fv.vector("TileProperty"));
1657     tileHEX1 = dbl_to_int(fv.vector("TileHEX1"));
1658     tileHEX2 = dbl_to_int(fv.vector("TileHEX2"));
1659     tileHEX3 = dbl_to_int(fv.vector("TileHEX3"));
1660     tileHEX4 = dbl_to_int(fv.vector("TileHEX4"));
1661     tileRMin = fv.vector("TileRMin");
1662     tileRMax = fv.vector("TileRMax");
1663     rescale(tileRMin, HGCalParameters::k_ScaleFromDDD);
1664     rescale(tileRMax, HGCalParameters::k_ScaleFromDDD);
1665     tileRingMin = dbl_to_int(fv.vector("TileRingMin"));
1666     tileRingMax = dbl_to_int(fv.vector("TileRingMax"));
1667     if (php.waferMaskMode_ == scintillatorCassette) {
1668       if (php.cassettes_ > 0)
1669         php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1670       cassetteShift = fv.vector("CassetteShiftHE");
1671       rescale(cassetteShift, HGCalParameters::k_ScaleFromDDD);
1672     }
1673 
1674     php.cassetteShift_ = cassetteShift;
1675     loadSpecParsTrapezoid(php,
1676                           tileIndx,
1677                           tileProperty,
1678                           tileHEX1,
1679                           tileHEX2,
1680                           tileHEX3,
1681                           tileHEX4,
1682                           tileRMin,
1683                           tileRMax,
1684                           tileRingMin,
1685                           tileRingMax);
1686   }
1687 }
1688 
1689 void HGCalGeomParameters::loadSpecParsTrapezoid(const cms::DDFilteredView& fv,
1690                                                 const cms::DDVectorsMap& vmap,
1691                                                 HGCalParameters& php,
1692                                                 const std::string& sdTag1) {
1693   for (auto const& it : vmap) {
1694     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1695       for (const auto& i : it.second)
1696         php.radiusMixBoundary_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1697     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1698       for (const auto& i : it.second)
1699         php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1700     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1701       for (const auto& i : it.second)
1702         php.layerCenter_.emplace_back(std::round(i));
1703     }
1704   }
1705 
1706   php.nPhiBinBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "NPhiBinBH"));
1707   php.layerFrontBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "LayerFrontBH"));
1708   php.rMinLayerBH_ = fv.get<std::vector<double> >(sdTag1, "RMinLayerBH");
1709   rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDD4hep);
1710   php.nCellsFine_ = php.nPhiBinBH_[0];
1711   php.nCellsCoarse_ = php.nPhiBinBH_[1];
1712   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1713   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1714 
1715   php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1716   php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1717   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDD4hep);
1718   php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1719   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDD4hep);
1720 
1721   php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1722   php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1723   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDD4hep);
1724   php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1725   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDD4hep);
1726   unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1727   for (unsigned int k = 0; k < kmax; ++k)
1728     php.slopeTop_.emplace_back(0);
1729 
1730   const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1731   php.layerOffset_ = dummy2[0];
1732 
1733   loadSpecParsTrapezoid(php);
1734 
1735   // tile parameters from Katja's file
1736   if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette)) {
1737     std::vector<int> tileIndx, tileProperty;
1738     std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4;
1739     std::vector<double> tileRMin, tileRMax;
1740     std::vector<int> tileRingMin, tileRingMax;
1741     std::vector<double> cassetteShift;
1742     for (auto const& it : vmap) {
1743       if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileIndex")) {
1744         for (const auto& i : it.second)
1745           tileIndx.emplace_back(std::round(i));
1746       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileProperty")) {
1747         for (const auto& i : it.second)
1748           tileProperty.emplace_back(std::round(i));
1749       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX1")) {
1750         for (const auto& i : it.second)
1751           tileHEX1.emplace_back(std::round(i));
1752       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX2")) {
1753         for (const auto& i : it.second)
1754           tileHEX2.emplace_back(std::round(i));
1755       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX3")) {
1756         for (const auto& i : it.second)
1757           tileHEX3.emplace_back(std::round(i));
1758       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX4")) {
1759         for (const auto& i : it.second)
1760           tileHEX4.emplace_back(std::round(i));
1761       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin")) {
1762         for (const auto& i : it.second)
1763           tileRMin.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1764       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax")) {
1765         for (const auto& i : it.second)
1766           tileRMax.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1767       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin")) {
1768         for (const auto& i : it.second)
1769           tileRingMin.emplace_back(std::round(i));
1770       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax")) {
1771         for (const auto& i : it.second)
1772           tileRingMax.emplace_back(std::round(i));
1773       }
1774     }
1775     if (php.waferMaskMode_ == scintillatorCassette) {
1776       for (auto const& it : vmap) {
1777         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1778           for (const auto& i : it.second)
1779             cassetteShift.emplace_back(i);
1780         }
1781       }
1782     }
1783 
1784     rescale(cassetteShift, HGCalParameters::k_ScaleFromDD4hep);
1785     php.cassetteShift_ = cassetteShift;
1786     loadSpecParsTrapezoid(php,
1787                           tileIndx,
1788                           tileProperty,
1789                           tileHEX1,
1790                           tileHEX2,
1791                           tileHEX3,
1792                           tileHEX4,
1793                           tileRMin,
1794                           tileRMax,
1795                           tileRingMin,
1796                           tileRingMax);
1797   }
1798 }
1799 
1800 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php) {
1801 #ifdef EDM_ML_DEBUG
1802   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters:nCells " << php.nCellsFine_ << ":" << php.nCellsCoarse_
1803                                 << " cellSize: " << php.cellSize_[0] << ":" << php.cellSize_[1];
1804   for (unsigned int k = 0; k < php.layerFrontBH_.size(); ++k)
1805     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Type[" << k << "] Front Layer = " << php.layerFrontBH_[k]
1806                                   << " rMin = " << php.rMinLayerBH_[k];
1807   for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k) {
1808     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k]
1809                                   << " Nphi = " << php.scintCells(k + php.firstLayer_)
1810                                   << " dPhi = " << php.scintCellSize(k + php.firstLayer_);
1811   }
1812 
1813   for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1814     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1815                                   << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1816 
1817   for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1818     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1819                                   << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1820 
1821   std::ostringstream st1;
1822   for (unsigned int k = 0; k < php.zRanges_.size(); ++k)
1823     st1 << ":" << php.zRanges_[k];
1824   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary[" << php.zRanges_.size() << "] " << st1.str();
1825 
1826   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1827                                 << php.layerCenter_.size();
1828   for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1829     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1830 #endif
1831 }
1832 
1833 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php,
1834                                                 const std::vector<int>& tileIndx,
1835                                                 const std::vector<int>& tileProperty,
1836                                                 const std::vector<int>& tileHEX1,
1837                                                 const std::vector<int>& tileHEX2,
1838                                                 const std::vector<int>& tileHEX3,
1839                                                 const std::vector<int>& tileHEX4,
1840                                                 const std::vector<double>& tileRMin,
1841                                                 const std::vector<double>& tileRMax,
1842                                                 const std::vector<int>& tileRingMin,
1843                                                 const std::vector<int>& tileRingMax) {
1844   // tile parameters from Katja's file
1845   for (unsigned int k = 0; k < tileIndx.size(); ++k) {
1846     php.tileInfoMap_[tileIndx[k]] = HGCalParameters::tileInfo(HGCalProperty::tileType(tileProperty[k]),
1847                                                               HGCalProperty::tileSiPM(tileProperty[k]),
1848                                                               tileHEX1[k],
1849                                                               tileHEX2[k],
1850                                                               tileHEX3[k],
1851                                                               tileHEX4[k]);
1852 #ifdef EDM_ML_DEBUG
1853     edm::LogVerbatim("HGCalGeom") << "Tile[" << k << ":" << tileIndx[k] << "] "
1854                                   << " Type " << HGCalProperty::tileType(tileProperty[k]) << " SiPM "
1855                                   << HGCalProperty::tileSiPM(tileProperty[k]) << " HEX " << std::hex << tileHEX1[k]
1856                                   << ":" << tileHEX2[k] << ":" << tileHEX3[k] << ":" << tileHEX4[k] << std::dec;
1857 #endif
1858   }
1859 
1860   for (unsigned int k = 0; k < tileRMin.size(); ++k) {
1861     php.tileRingR_.emplace_back(tileRMin[k], tileRMax[k]);
1862 #ifdef EDM_ML_DEBUG
1863     edm::LogVerbatim("HGCalGeom") << "TileRingR[" << k << "] " << tileRMin[k] << ":" << tileRMax[k];
1864 #endif
1865   }
1866 
1867   for (unsigned k = 0; k < tileRingMin.size(); ++k) {
1868     php.tileRingRange_.emplace_back(tileRingMin[k], tileRingMax[k]);
1869 #ifdef EDM_ML_DEBUG
1870     edm::LogVerbatim("HGCalGeom") << "TileRingRange[" << k << "] " << tileRingMin[k] << ":" << tileRingMax[k];
1871 #endif
1872   }
1873 }
1874 
1875 void HGCalGeomParameters::loadWaferHexagon(HGCalParameters& php) {
1876   double waferW(HGCalParameters::k_ScaleFromDDD * waferSize_), rmin(HGCalParameters::k_ScaleFromDDD * php.waferR_);
1877   double rin(php.rLimit_[0]), rout(php.rLimit_[1]), rMaxFine(php.boundR_[1]);
1878 #ifdef EDM_ML_DEBUG
1879   edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << rmin << " R Limits: " << rin << ":" << rout
1880                                 << " Fine " << rMaxFine;
1881 #endif
1882   // Clear the vectors
1883   php.waferCopy_.clear();
1884   php.waferTypeL_.clear();
1885   php.waferTypeT_.clear();
1886   php.waferPosX_.clear();
1887   php.waferPosY_.clear();
1888   double dx = 0.5 * waferW;
1889   double dy = 3.0 * dx * tan(30._deg);
1890   double rr = 2.0 * dx * tan(30._deg);
1891   int ncol = static_cast<int>(2.0 * rout / waferW) + 1;
1892   int nrow = static_cast<int>(rout / (waferW * tan(30._deg))) + 1;
1893   int ns2 = (2 * ncol + 1) * (2 * nrow + 1) * php.layer_.size();
1894   int incm(0), inrm(0);
1895   HGCalParameters::layer_map copiesInLayers(php.layer_.size() + 1);
1896   HGCalParameters::waferT_map waferTypes(ns2 + 1);
1897 #ifdef EDM_ML_DEBUG
1898   int kount(0), ntot(0);
1899   edm::LogVerbatim("HGCalGeom") << "Row " << nrow << " Column " << ncol;
1900 #endif
1901   for (int nr = -nrow; nr <= nrow; ++nr) {
1902     int inr = (nr >= 0) ? nr : -nr;
1903     for (int nc = -ncol; nc <= ncol; ++nc) {
1904       int inc = (nc >= 0) ? nc : -nc;
1905       if (inr % 2 == inc % 2) {
1906         double xpos = nc * dx;
1907         double ypos = nr * dy;
1908         std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, rin, rout, true);
1909         double rpos = std::sqrt(xpos * xpos + ypos * ypos);
1910         int typet = (rpos < rMaxFine) ? 1 : 2;
1911         int typel(3);
1912         for (int k = 1; k < 4; ++k) {
1913           if ((rpos + rmin) <= php.boundR_[k]) {
1914             typel = k;
1915             break;
1916           }
1917         }
1918 #ifdef EDM_ML_DEBUG
1919         ++ntot;
1920 #endif
1921         if (corner.first > 0) {
1922           int copy = HGCalTypes::packTypeUV(typel, nc, nr);
1923           if (inc > incm)
1924             incm = inc;
1925           if (inr > inrm)
1926             inrm = inr;
1927 #ifdef EDM_ML_DEBUG
1928           kount++;
1929           edm::LogVerbatim("HGCalGeom") << kount << ":" << ntot << " Copy " << copy << " Type " << typel << ":" << typet
1930                                         << " Location " << corner.first << " Position " << xpos << ":" << ypos
1931                                         << " Layers " << php.layer_.size();
1932 #endif
1933           php.waferCopy_.emplace_back(copy);
1934           php.waferTypeL_.emplace_back(typel);
1935           php.waferTypeT_.emplace_back(typet);
1936           php.waferPosX_.emplace_back(xpos);
1937           php.waferPosY_.emplace_back(ypos);
1938           for (unsigned int il = 0; il < php.layer_.size(); ++il) {
1939             std::pair<int, int> corner =
1940                 HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, php.rMinLayHex_[il], php.rMaxLayHex_[il], true);
1941             if (corner.first > 0) {
1942               auto cpy = copiesInLayers[php.layer_[il]].find(copy);
1943               if (cpy == copiesInLayers[php.layer_[il]].end())
1944                 copiesInLayers[php.layer_[il]][copy] =
1945                     ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ? php.waferCopy_.size() : -1);
1946             }
1947             if ((corner.first > 0) && (corner.first < static_cast<int>(HGCalParameters::k_CornerSize))) {
1948               int wl = HGCalWaferIndex::waferIndex(php.layer_[il], copy, 0, true);
1949               waferTypes[wl] = corner;
1950             }
1951           }
1952         }
1953       }
1954     }
1955   }
1956   php.copiesInLayers_ = copiesInLayers;
1957   php.waferTypes_ = waferTypes;
1958   php.nSectors_ = static_cast<int>(php.waferCopy_.size());
1959   php.waferUVMax_ = 0;
1960 #ifdef EDM_ML_DEBUG
1961   edm::LogVerbatim("HGCalGeom") << "HGCalWaferHexagon: # of columns " << incm << " # of rows " << inrm << " and "
1962                                 << kount << ":" << ntot << " wafers; R " << rin << ":" << rout;
1963   edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
1964   for (unsigned int k = 0; k < copiesInLayers.size(); ++k) {
1965     const auto& theModules = copiesInLayers[k];
1966     edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
1967     int k2(0);
1968     for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
1969       edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
1970     }
1971   }
1972 #endif
1973 }
1974 
1975 void HGCalGeomParameters::loadWaferHexagon8(HGCalParameters& php) {
1976   double waferW(php.waferSize_);
1977   double waferS(php.sensorSeparation_);
1978   auto wType = std::make_unique<HGCalWaferType>(php.radius100to200_,
1979                                                 php.radius200to300_,
1980                                                 HGCalParameters::k_ScaleToDDD * (waferW + waferS),
1981                                                 HGCalParameters::k_ScaleToDDD * php.zMinForRad_,
1982                                                 php.choiceType_,
1983                                                 php.nCornerCut_,
1984                                                 php.fracAreaMin_);
1985 
1986   double rout(php.rLimit_[1]);
1987 #ifdef EDM_ML_DEBUG
1988   edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << waferS << " R Max: " << rout;
1989 #endif
1990   // Clear the vectors
1991   php.waferCopy_.clear();
1992   php.waferTypeL_.clear();
1993   php.waferTypeT_.clear();
1994   php.waferPosX_.clear();
1995   php.waferPosY_.clear();
1996   double r = 0.5 * (waferW + waferS);
1997   double R = 2.0 * r / sqrt3_;
1998   double dy = 0.75 * R;
1999   double r1 = 0.5 * waferW;
2000   double R1 = 2.0 * r1 / sqrt3_;
2001   int N = (r == 0) ? 3 : (static_cast<int>(0.5 * rout / r) + 3);
2002   int ns1 = (2 * N + 1) * (2 * N + 1);
2003   int ns2 = ns1 * php.zLayerHex_.size();
2004 #ifdef EDM_ML_DEBUG
2005   edm::LogVerbatim("HGCalGeom") << "wafer " << waferW << ":" << waferS << " r " << r << " dy " << dy << " N " << N
2006                                 << " sizes " << ns1 << ":" << ns2;
2007   std::vector<int> indtypes(ns1 + 1);
2008   indtypes.clear();
2009 #endif
2010   HGCalParameters::wafer_map wafersInLayers(ns1 + 1);
2011   HGCalParameters::wafer_map typesInLayers(ns2 + 1);
2012   HGCalParameters::waferT_map waferTypes(ns2 + 1);
2013   int ipos(0), lpos(0), uvmax(0), nwarn(0);
2014   std::vector<int> uvmx(php.zLayerHex_.size(), 0);
2015   for (int v = -N; v <= N; ++v) {
2016     for (int u = -N; u <= N; ++u) {
2017       int nr = 2 * v;
2018       int nc = -2 * u + v;
2019       double xpos = nc * r;
2020       double ypos = nr * dy;
2021       int indx = HGCalWaferIndex::waferIndex(0, u, v);
2022       php.waferCopy_.emplace_back(indx);
2023       php.waferPosX_.emplace_back(xpos);
2024       php.waferPosY_.emplace_back(ypos);
2025       wafersInLayers[indx] = ipos;
2026       ++ipos;
2027       std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, r1, R1, 0, rout, false);
2028       if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2029           ((corner.first > 0) && php.defineFull_)) {
2030         uvmax = std::max(uvmax, std::max(std::abs(u), std::abs(v)));
2031       }
2032       for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
2033         int copy = i + php.layerOffset_;
2034         std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], (waferW + waferS));
2035         int lay = php.layer_[php.layerIndex_[i]];
2036         double xpos0 = xpos + xyoff.first;
2037         double ypos0 = ypos + xyoff.second;
2038         double zpos = php.zLayerHex_[i];
2039         int kndx = HGCalWaferIndex::waferIndex(lay, u, v);
2040         int type(-1);
2041         if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2042             (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) || (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
2043           type = wType->getType(kndx, php.waferInfoMap_);
2044         if (type < 0)
2045           type = wType->getType(HGCalParameters::k_ScaleToDDD * xpos0,
2046                                 HGCalParameters::k_ScaleToDDD * ypos0,
2047                                 HGCalParameters::k_ScaleToDDD * zpos);
2048         php.waferTypeL_.emplace_back(type);
2049         typesInLayers[kndx] = lpos;
2050         ++lpos;
2051 #ifdef EDM_ML_DEBUG
2052         indtypes.emplace_back(kndx);
2053 #endif
2054         std::pair<int, int> corner =
2055             HGCalGeomTools::waferCorner(xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], false);
2056 #ifdef EDM_ML_DEBUG
2057         if (((corner.first == 0) && std::abs(u) < 5 && std::abs(v) < 5) || (std::abs(u) < 2 && std::abs(v) < 2)) {
2058           edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " R " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2059                                         << " u " << u << " v " << v << " with " << corner.first << " corners";
2060         }
2061 #endif
2062         if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2063             ((corner.first > 0) && php.defineFull_)) {
2064           uvmx[i] = std::max(uvmx[i], std::max(std::abs(u), std::abs(v)));
2065         }
2066         if ((corner.first < static_cast<int>(HGCalParameters::k_CornerSize)) && (corner.first > 0)) {
2067 #ifdef EDM_ML_DEBUG
2068           edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2069                                         << corner.first << ":" << corner.second;
2070 #endif
2071           int wl = HGCalWaferIndex::waferIndex(lay, u, v);
2072           if (php.waferMaskMode_ > 0) {
2073             bool v17OrLess = (php.mode_ < HGCalGeometryMode::Hexagon8CalibCell);
2074             std::pair<int, int> corner0 = HGCalWaferMask::getTypeMode(
2075                 xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], type, php.waferMaskMode_, v17OrLess);
2076             if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2077                 (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2078                 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
2079               auto itr = php.waferInfoMap_.find(wl);
2080               if (itr != php.waferInfoMap_.end()) {
2081                 int part = (itr->second).part;
2082                 int orient = (itr->second).orient;
2083                 bool ok = ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2084                            (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
2085                               ? true
2086                               : HGCalWaferMask::goodTypeMode(
2087                                     xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], part, orient, false);
2088                 if (ok)
2089                   corner0 = std::make_pair(part, (HGCalTypes::k_OffsetRotation + orient));
2090 #ifdef EDM_ML_DEBUG
2091                 edm::LogVerbatim("HGCalGeom")
2092                     << "Layer:u:v " << i << ":" << lay << ":" << u << ":" << v << " Part " << corner0.first << ":"
2093                     << part << " Orient " << corner0.second << ":" << orient << " Position " << xpos0 << ":" << ypos0
2094                     << " delta " << r1 << ":" << R1 << " Limit " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2095                     << " Compatibiliety Flag " << ok;
2096 #endif
2097                 if (!ok)
2098                   ++nwarn;
2099               }
2100             }
2101             waferTypes[wl] = corner0;
2102 #ifdef EDM_ML_DEBUG
2103             edm::LogVerbatim("HGCalGeom")
2104                 << "Layer " << lay << " u|v " << u << ":" << v << " Index " << std::hex << wl << std::dec << " pos "
2105                 << xpos0 << ":" << ypos0 << " R " << r1 << ":" << R1 << " Range " << php.rMinLayHex_[i] << ":"
2106                 << php.rMaxLayHex_[i] << type << ":" << php.waferMaskMode_ << " corner " << corner.first << ":"
2107                 << corner.second << " croner0 " << corner0.first << ":" << corner0.second;
2108 #endif
2109           } else {
2110             waferTypes[wl] = corner;
2111 #ifdef EDM_ML_DEBUG
2112             edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2113                                           << corner.first << ":" << corner.second;
2114 #endif
2115           }
2116         }
2117       }
2118     }
2119   }
2120   if (nwarn > 0)
2121     edm::LogWarning("HGCalGeom") << "HGCalGeomParameters::loadWafer8: there are " << nwarn
2122                                  << " wafers with non-matching partial- orientation types";
2123   php.waferUVMax_ = uvmax;
2124   php.waferUVMaxLayer_ = uvmx;
2125   php.wafersInLayers_ = wafersInLayers;
2126   php.typesInLayers_ = typesInLayers;
2127   php.waferTypes_ = waferTypes;
2128   php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2129   HGCalParameters::hgtrap mytr;
2130   mytr.lay = 1;
2131   mytr.bl = php.waferR_;
2132   mytr.tl = php.waferR_;
2133   mytr.h = php.waferR_;
2134   mytr.alpha = 0.0;
2135   mytr.cellSize = HGCalParameters::k_ScaleToDDD * php.waferSize_;
2136   for (auto const& dz : php.cellThickness_) {
2137     mytr.dz = 0.5 * HGCalParameters::k_ScaleToDDD * dz;
2138     php.fillModule(mytr, false);
2139   }
2140   for (unsigned k = 0; k < php.cellThickness_.size(); ++k) {
2141     HGCalParameters::hgtrap mytr = php.getModule(k, false);
2142     mytr.bl *= HGCalParameters::k_ScaleFromDDD;
2143     mytr.tl *= HGCalParameters::k_ScaleFromDDD;
2144     mytr.h *= HGCalParameters::k_ScaleFromDDD;
2145     mytr.dz *= HGCalParameters::k_ScaleFromDDD;
2146     mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2147     php.fillModule(mytr, true);
2148   }
2149 #ifdef EDM_ML_DEBUG
2150   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Total of " << php.waferCopy_.size() << " wafers";
2151   for (unsigned int k = 0; k < php.waferCopy_.size(); ++k) {
2152     int id = php.waferCopy_[k];
2153     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << std::hex << id << std::dec << ":"
2154                                   << HGCalWaferIndex::waferLayer(id) << ":" << HGCalWaferIndex::waferU(id) << ":"
2155                                   << HGCalWaferIndex::waferV(id) << " x " << php.waferPosX_[k] << " y "
2156                                   << php.waferPosY_[k] << " index " << php.wafersInLayers_[id];
2157   }
2158   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Total of " << php.waferTypeL_.size() << " wafer types";
2159   for (unsigned int k = 0; k < php.waferTypeL_.size(); ++k) {
2160     int id = indtypes[k];
2161     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.typesInLayers_[id] << ":" << php.waferTypeL_[k] << " ID "
2162                                   << std::hex << id << std::dec << ":" << HGCalWaferIndex::waferLayer(id) << ":"
2163                                   << HGCalWaferIndex::waferU(id) << ":" << HGCalWaferIndex::waferV(id);
2164   }
2165 #endif
2166 
2167   //Wafer offset
2168   php.xLayerHex_.clear();
2169   php.yLayerHex_.clear();
2170   double waferSize = php.waferSize_ + php.sensorSeparation_;
2171 #ifdef EDM_ML_DEBUG
2172   edm::LogVerbatim("HGCalGeom") << "WaferSize " << waferSize;
2173 #endif
2174   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2175     int copy = k + php.layerOffset_;
2176     std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], waferSize);
2177     php.xLayerHex_.emplace_back(xyoff.first);
2178     php.yLayerHex_.emplace_back(xyoff.second);
2179 #ifdef EDM_ML_DEBUG
2180     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Off " << copy << ":" << php.layerCenter_[copy] << " Shift "
2181                                   << xyoff.first << ":" << xyoff.second;
2182 #endif
2183   }
2184 }
2185 
2186 void HGCalGeomParameters::loadCellParsHexagon(const DDCompactView* cpv, HGCalParameters& php) {
2187   // Special parameters for cell parameters
2188   std::string attribute = "OnlyForHGCalNumbering";
2189   DDSpecificsHasNamedValueFilter filter1{attribute};
2190   DDFilteredView fv1(*cpv, filter1);
2191   bool ok = fv1.firstChild();
2192 
2193   if (ok) {
2194     php.cellFine_ = dbl_to_int(cpv->vector("waferFine"));
2195     php.cellCoarse_ = dbl_to_int(cpv->vector("waferCoarse"));
2196   }
2197 
2198   loadCellParsHexagon(php);
2199 }
2200 
2201 void HGCalGeomParameters::loadCellParsHexagon(const cms::DDVectorsMap& vmap, HGCalParameters& php) {
2202   for (auto const& it : vmap) {
2203     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferFine")) {
2204       for (const auto& i : it.second)
2205         php.cellFine_.emplace_back(std::round(i));
2206     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferCoarse")) {
2207       for (const auto& i : it.second)
2208         php.cellCoarse_.emplace_back(std::round(i));
2209     }
2210   }
2211 
2212   loadCellParsHexagon(php);
2213 }
2214 
2215 void HGCalGeomParameters::loadCellParsHexagon(const HGCalParameters& php) {
2216 #ifdef EDM_ML_DEBUG
2217   edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellFine_.size() << " rows for fine cells";
2218   for (unsigned int k = 0; k < php.cellFine_.size(); ++k)
2219     edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellFine_[k];
2220   edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellCoarse_.size() << " rows for coarse cells";
2221   for (unsigned int k = 0; k < php.cellCoarse_.size(); ++k)
2222     edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellCoarse_[k];
2223 #endif
2224 }
2225 
2226 void HGCalGeomParameters::loadCellTrapezoid(HGCalParameters& php) {
2227   php.xLayerHex_.resize(php.zLayerHex_.size(), 0);
2228   php.yLayerHex_.resize(php.zLayerHex_.size(), 0);
2229 #ifdef EDM_ML_DEBUG
2230   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: x|y|zLayerHex in array of size " << php.zLayerHex_.size();
2231   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k)
2232     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Shift " << php.xLayerHex_[k] << ":" << php.yLayerHex_[k] << ":"
2233                                   << php.zLayerHex_[k];
2234 #endif
2235   // Find the radius of each eta-partitions
2236 
2237   if ((php.mode_ == HGCalGeometryMode::TrapezoidFile) || (php.mode_ == HGCalGeometryMode::TrapezoidModule) ||
2238       (php.mode_ == HGCalGeometryMode::TrapezoidCassette)) {
2239     //Ring radii for each partition
2240     for (unsigned int k = 0; k < 2; ++k) {
2241       for (unsigned int kk = 0; kk < php.tileRingR_.size(); ++kk) {
2242         php.radiusLayer_[k].emplace_back(php.tileRingR_[kk].first);
2243 #ifdef EDM_ML_DEBUG
2244         double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2245                               : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2246         double rv = php.radiusLayer_[k].back();
2247         double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2248         edm::LogVerbatim("HGCalGeom") << "New [" << kk << "] new R = " << rv << " Eta = " << eta;
2249 #endif
2250       }
2251       php.radiusLayer_[k].emplace_back(php.tileRingR_[php.tileRingR_.size() - 1].second);
2252     }
2253     // Minimum and maximum radius index for each layer
2254     for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2255       php.iradMinBH_.emplace_back(1 + php.tileRingRange_[k].first);
2256       php.iradMaxBH_.emplace_back(1 + php.tileRingRange_[k].second);
2257 #ifdef EDM_ML_DEBUG
2258       int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2259       edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge " << php.iradMinBH_.back()
2260                                     << " Top edge " << php.iradMaxBH_.back();
2261 #endif
2262     }
2263   } else {
2264     //Ring radii for each partition
2265     for (unsigned int k = 0; k < 2; ++k) {
2266       double rmax = ((k == 0) ? (php.rMaxLayHex_[php.layerFrontBH_[1] - php.firstLayer_] - 1)
2267                               : (php.rMaxLayHex_[php.rMaxLayHex_.size() - 1]));
2268       double rv = php.rMinLayerBH_[k];
2269       double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2270                             : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2271       php.radiusLayer_[k].emplace_back(rv);
2272 #ifdef EDM_ML_DEBUG
2273       double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2274       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] rmax " << rmax << " Z = " << zv
2275                                     << " dEta = " << php.cellSize_[k] << "\nOld[0] new R = " << rv << " Eta = " << eta;
2276       int kount(1);
2277 #endif
2278       while (rv < rmax) {
2279         double eta = -(php.cellSize_[k] + std::log(std::tan(0.5 * std::atan(rv / zv))));
2280         rv = zv * std::tan(2.0 * std::atan(std::exp(-eta)));
2281         php.radiusLayer_[k].emplace_back(rv);
2282 #ifdef EDM_ML_DEBUG
2283         edm::LogVerbatim("HGCalGeom") << "Old [" << kount << "] new R = " << rv << " Eta = " << eta;
2284         ++kount;
2285 #endif
2286       }
2287     }
2288     // Find minimum and maximum radius index for each layer
2289     for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2290       int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2291       std::vector<double>::iterator low, high;
2292       low = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMinLayHex_[k]);
2293 #ifdef EDM_ML_DEBUG
2294       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RLow = " << php.rMinLayHex_[k] << " pos "
2295                                     << static_cast<int>(low - php.radiusLayer_[kk].begin());
2296 #endif
2297       if (low == php.radiusLayer_[kk].begin())
2298         ++low;
2299       int irlow = static_cast<int>(low - php.radiusLayer_[kk].begin());
2300       double drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2301 #ifdef EDM_ML_DEBUG
2302       edm::LogVerbatim("HGCalGeom") << "irlow " << irlow << " dr " << drlow << " min " << php.minTileSize_;
2303 #endif
2304       if (drlow < php.minTileSize_) {
2305         ++irlow;
2306 #ifdef EDM_ML_DEBUG
2307         drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2308         edm::LogVerbatim("HGCalGeom") << "Modified irlow " << irlow << " dr " << drlow;
2309 #endif
2310       }
2311       high = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMaxLayHex_[k]);
2312 #ifdef EDM_ML_DEBUG
2313       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RHigh = " << php.rMaxLayHex_[k] << " pos "
2314                                     << static_cast<int>(high - php.radiusLayer_[kk].begin());
2315 #endif
2316       if (high == php.radiusLayer_[kk].end())
2317         --high;
2318       int irhigh = static_cast<int>(high - php.radiusLayer_[kk].begin());
2319       double drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2320 #ifdef EDM_ML_DEBUG
2321       edm::LogVerbatim("HGCalGeom") << "irhigh " << irhigh << " dr " << drhigh << " min " << php.minTileSize_;
2322 #endif
2323       if (drhigh < php.minTileSize_) {
2324         --irhigh;
2325 #ifdef EDM_ML_DEBUG
2326         drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2327         edm::LogVerbatim("HGCalGeom") << "Modified irhigh " << irhigh << " dr " << drhigh;
2328 #endif
2329       }
2330       php.iradMinBH_.emplace_back(irlow);
2331       php.iradMaxBH_.emplace_back(irhigh);
2332 #ifdef EDM_ML_DEBUG
2333       edm::LogVerbatim("HGCalGeom") << "Old Layer " << k << " Type " << kk << " Low edge " << irlow << ":" << drlow
2334                                     << " Top edge " << irhigh << ":" << drhigh;
2335 #endif
2336     }
2337   }
2338 #ifdef EDM_ML_DEBUG
2339   for (unsigned int k = 0; k < 2; ++k) {
2340     edm::LogVerbatim("HGCalGeom") << "Type " << k << " with " << php.radiusLayer_[k].size() << " radii";
2341     for (unsigned int kk = 0; kk < php.radiusLayer_[k].size(); ++kk)
2342       edm::LogVerbatim("HGCalGeom") << "Ring[" << kk << "] " << php.radiusLayer_[k][kk];
2343   }
2344 #endif
2345 
2346   // Now define the volumes
2347   int im(0);
2348   php.waferUVMax_ = 0;
2349   HGCalParameters::hgtrap mytr;
2350   mytr.alpha = 0.0;
2351   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2352     if (php.iradMaxBH_[k] > php.waferUVMax_)
2353       php.waferUVMax_ = php.iradMaxBH_[k];
2354     int kk = ((php.firstLayer_ + static_cast<int>(k)) < php.layerFrontBH_[1]) ? 0 : 1;
2355     int irm = php.radiusLayer_[kk].size() - 1;
2356 #ifdef EDM_ML_DEBUG
2357     double rmin = php.radiusLayer_[kk][std::max((php.iradMinBH_[k] - 1), 0)];
2358     double rmax = php.radiusLayer_[kk][std::min(php.iradMaxBH_[k], irm)];
2359     edm::LogVerbatim("HGCalGeom") << "Layer " << php.firstLayer_ + k << ":" << kk << " Radius range "
2360                                   << php.iradMinBH_[k] << ":" << php.iradMaxBH_[k] << ":" << rmin << ":" << rmax;
2361 #endif
2362     mytr.lay = php.firstLayer_ + k;
2363     for (int irad = php.iradMinBH_[k]; irad <= php.iradMaxBH_[k]; ++irad) {
2364       double rmin = php.radiusLayer_[kk][std::max((irad - 1), 0)];
2365       double rmax = php.radiusLayer_[kk][std::min(irad, irm)];
2366       mytr.bl = 0.5 * rmin * php.scintCellSize(mytr.lay);
2367       mytr.tl = 0.5 * rmax * php.scintCellSize(mytr.lay);
2368       mytr.h = 0.5 * (rmax - rmin);
2369       mytr.dz = 0.5 * php.waferThick_;
2370       mytr.cellSize = 0.5 * (rmax + rmin) * php.scintCellSize(mytr.lay);
2371       php.fillModule(mytr, true);
2372       mytr.bl *= HGCalParameters::k_ScaleToDDD;
2373       mytr.tl *= HGCalParameters::k_ScaleToDDD;
2374       mytr.h *= HGCalParameters::k_ScaleToDDD;
2375       mytr.dz *= HGCalParameters::k_ScaleToDDD;
2376       mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2377       php.fillModule(mytr, false);
2378       if (irad == php.iradMinBH_[k])
2379         php.firstModule_.emplace_back(im);
2380       ++im;
2381       if (irad == php.iradMaxBH_[k] - 1)
2382         php.lastModule_.emplace_back(im);
2383     }
2384   }
2385   php.nSectors_ = php.waferUVMax_;
2386 #ifdef EDM_ML_DEBUG
2387   edm::LogVerbatim("HGCalGeom") << "Maximum radius index " << php.waferUVMax_;
2388   for (unsigned int k = 0; k < php.firstModule_.size(); ++k)
2389     edm::LogVerbatim("HGCalGeom") << "Layer " << k + php.firstLayer_ << " Modules " << php.firstModule_[k] << ":"
2390                                   << php.lastModule_[k];
2391 #endif
2392 }
2393 
2394 std::vector<double> HGCalGeomParameters::getDDDArray(const std::string& str, const DDsvalues_type& sv, const int nmin) {
2395   DDValue value(str);
2396   if (DDfetch(&sv, value)) {
2397     const std::vector<double>& fvec = value.doubles();
2398     int nval = fvec.size();
2399     if (nmin > 0) {
2400       if (nval < nmin) {
2401         throw cms::Exception("DDException")
2402             << "HGCalGeomParameters:  # of " << str << " bins " << nval << " < " << nmin << " ==> illegal";
2403       }
2404     } else {
2405       if (nval < 1 && nmin == 0) {
2406         throw cms::Exception("DDException")
2407             << "HGCalGeomParameters: # of " << str << " bins " << nval << " < 1 ==> illegal"
2408             << " (nmin=" << nmin << ")";
2409       }
2410     }
2411     return fvec;
2412   } else {
2413     if (nmin >= 0) {
2414       throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str;
2415     }
2416     std::vector<double> fvec;
2417     return fvec;
2418   }
2419 }
2420 
2421 std::pair<double, double> HGCalGeomParameters::cellPosition(
2422     const std::vector<HGCalGeomParameters::cellParameters>& wafers,
2423     std::vector<HGCalGeomParameters::cellParameters>::const_iterator& itrf,
2424     int wafer,
2425     double xx,
2426     double yy) {
2427   if (itrf == wafers.end()) {
2428     for (std::vector<HGCalGeomParameters::cellParameters>::const_iterator itr = wafers.begin(); itr != wafers.end();
2429          ++itr) {
2430       if (itr->wafer == wafer) {
2431         itrf = itr;
2432         break;
2433       }
2434     }
2435   }
2436   double dx(0), dy(0);
2437   if (itrf != wafers.end()) {
2438     dx = (xx - itrf->xyz.x());
2439     if (std::abs(dx) < tolerance)
2440       dx = 0;
2441     dy = (yy - itrf->xyz.y());
2442     if (std::abs(dy) < tolerance)
2443       dy = 0;
2444   }
2445   return std::make_pair(dx, dy);
2446 }
2447 
2448 void HGCalGeomParameters::rescale(std::vector<double>& v, const double s) {
2449   std::for_each(v.begin(), v.end(), [s](double& n) { n *= s; });
2450 }
2451 
2452 void HGCalGeomParameters::resetZero(std::vector<double>& v) {
2453   for (auto& n : v) {
2454     if (std::abs(n) < tolmin)
2455       n = 0;
2456   }
2457 }