Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-14 22:20:52 

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) || (php.mode_ == HGCalGeometryMode::Hexagon8FineCell)) {
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) || (php.mode_ == HGCalGeometryMode::Hexagon8FineCell)) {
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) || (php.mode_ == HGCalGeometryMode::Hexagon8FineCell)) {
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) || (php.mode_ == HGCalGeometryMode::Hexagon8FineCell)) {
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) || (php.mode_ == HGCalGeometryMode::Hexagon8FineCell)) {
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          (php.mode_ == HGCalGeometryMode::Hexagon8FineCell))
1596             ? HGCalProperty::waferOrient(waferProperties[k])
1597             : HGCalWaferMask::getRotation(php.waferZSide_, partial, HGCalProperty::waferOrient(waferProperties[k]));
1598     php.waferInfoMap_[waferIndex[k]] = HGCalParameters::waferInfo(HGCalProperty::waferThick(waferProperties[k]),
1599                                                                   partial,
1600                                                                   orient,
1601                                                                   HGCalProperty::waferCassette(waferProperties[k]));
1602 #ifdef EDM_ML_DEBUG
1603     edm::LogVerbatim("HGCalGeom") << "[" << k << ":" << waferIndex[k] << ":"
1604                                   << HGCalWaferIndex::waferLayer(waferIndex[k]) << ":"
1605                                   << HGCalWaferIndex::waferU(waferIndex[k]) << ":"
1606                                   << HGCalWaferIndex::waferV(waferIndex[k]) << "]  Thickness type "
1607                                   << HGCalProperty::waferThick(waferProperties[k]) << " Partial type " << partial
1608                                   << " Orientation " << HGCalProperty::waferOrient(waferProperties[k]) << ":" << orient;
1609 #endif
1610   }
1611 }
1612 
1613 void HGCalGeomParameters::loadSpecParsTrapezoid(const DDFilteredView& fv, HGCalParameters& php) {
1614   DDsvalues_type sv(fv.mergedSpecifics());
1615   php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1616   rescale(php.radiusMixBoundary_, HGCalParameters::k_ScaleFromDDD);
1617 
1618   php.nPhiBinBH_ = dbl_to_int(getDDDArray("NPhiBinBH", sv, 0));
1619   php.layerFrontBH_ = dbl_to_int(getDDDArray("LayerFrontBH", sv, 0));
1620   php.rMinLayerBH_ = getDDDArray("RMinLayerBH", sv, 0);
1621   rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDDD);
1622   assert(php.nPhiBinBH_.size() > 1);
1623   php.nCellsFine_ = php.nPhiBinBH_[1];
1624   php.nCellsCoarse_ = php.nPhiBinBH_[0];
1625   assert(0 != php.nCellsFine_);
1626   assert(0 != php.nCellsCoarse_);
1627   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1628   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1629 
1630   php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1631   php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1632   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDDD);
1633   php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1634   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDDD);
1635 
1636   php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1637   php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1638   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDDD);
1639   php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1640   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDDD);
1641 
1642   php.zRanges_ = fv.vector("ZRanges");
1643   rescale(php.zRanges_, HGCalParameters::k_ScaleFromDDD);
1644 
1645   // Offsets
1646   const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1647   php.layerOffset_ = dummy2[0];
1648   php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1649 
1650   loadSpecParsTrapezoid(php);
1651 #ifdef EDM_ML_DEBUG
1652   edm::LogVerbatim("HGCalGeom") << "WaferMaskMode " << php.waferMaskMode_ << " Compare " << scintillatorFile << ":"
1653                                 << scintillatorCassette << ":" << scintillatorFineCell;
1654 #endif
1655   // tile parameters from Katja's file
1656   if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette) ||
1657       (php.waferMaskMode_ == scintillatorFineCell)) {
1658     std::vector<int> tileIndx, tileProperty;
1659     std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4, tileHEX5, tileHEX6;
1660     std::vector<double> tileRMin, tileRMax, tileRMinFine, tileRMaxFine;
1661     std::vector<int> tileRingMin, tileRingMax, tileRingMinFine, tileRingMaxFine;
1662     std::vector<double> cassetteShift;
1663     php.nPhiLayer_ = dbl_to_int(fv.vector("NPhiLayer"));
1664     tileIndx = dbl_to_int(fv.vector("TileIndex"));
1665     tileProperty = dbl_to_int(fv.vector("TileProperty"));
1666     tileHEX1 = dbl_to_int(fv.vector("TileHEX1"));
1667     tileHEX2 = dbl_to_int(fv.vector("TileHEX2"));
1668     tileHEX3 = dbl_to_int(fv.vector("TileHEX3"));
1669     tileHEX4 = dbl_to_int(fv.vector("TileHEX4"));
1670     tileRMin = fv.vector("TileRMin");
1671     tileRMax = fv.vector("TileRMax");
1672     rescale(tileRMin, HGCalParameters::k_ScaleFromDDD);
1673     rescale(tileRMax, HGCalParameters::k_ScaleFromDDD);
1674     tileRingMin = dbl_to_int(fv.vector("TileRingMin"));
1675     tileRingMax = dbl_to_int(fv.vector("TileRingMax"));
1676     if (php.waferMaskMode_ == scintillatorFineCell) {
1677       tileHEX5 = dbl_to_int(fv.vector("TileHEX5"));
1678       tileHEX6 = dbl_to_int(fv.vector("TileHEX6"));
1679       tileRMinFine = fv.vector("TileRMin6");
1680       tileRMaxFine = fv.vector("TileRMax6");
1681       rescale(tileRMinFine, HGCalParameters::k_ScaleFromDDD);
1682       rescale(tileRMaxFine, HGCalParameters::k_ScaleFromDDD);
1683       tileRingMinFine = dbl_to_int(fv.vector("TileRingMin6"));
1684       tileRingMaxFine = dbl_to_int(fv.vector("TileRingMax6"));
1685       php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1686       php.nphiFineCassette_ = php.nCellsFine_ / php.cassettes_;
1687       std::vector<double> rectract = fv.vector("ScintRetract");
1688       rescale(rectract, HGCalParameters::k_ScaleFromDDD);
1689       for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1690         php.cassetteRetractTile_.emplace_back(rectract[k1]);
1691 #ifdef EDM_ML_DEBUG
1692         edm::LogVerbatim("HGCalGeom") << "cassetteRetractTile_[" << k1 << "] " << rectract[k1];
1693 #endif
1694       }
1695       int n = 2 * php.cassettes_ * (php.firstLayer_ - 1);
1696       for (int k1 = 0; k1 < n; ++k1)
1697         cassetteShift.emplace_back(0.);
1698 #ifdef EDM_ML_DEBUG
1699       edm::LogVerbatim("HGCalGeom") << "First Layer " << php.firstLayer_ << " Rett size " << rectract.size() << " N "
1700                                     << n;
1701 #endif
1702       double dphi = M_PI / php.cassettes_;
1703       for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1704         for (int k2 = 0; k2 < php.cassettes_; ++k2) {
1705           double phi = (2 * k2 + 1) * dphi;
1706           cassetteShift.emplace_back(rectract[k1] * cos(phi));
1707           cassetteShift.emplace_back(rectract[k1] * sin(phi));
1708         }
1709       }
1710     } else if (php.waferMaskMode_ == scintillatorCassette) {
1711       if (php.cassettes_ > 0)
1712         php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1713       cassetteShift = fv.vector("CassetteShiftHE");
1714     }
1715 
1716     rescale(cassetteShift, HGCalParameters::k_ScaleFromDDD);
1717     if (php.waferMaskMode_ == scintillatorFineCell)
1718       php.cassetteShiftTile_ = cassetteShift;
1719     else
1720       php.cassetteShift_ = cassetteShift;
1721     loadSpecParsTrapezoid(php,
1722                           tileIndx,
1723                           tileProperty,
1724                           tileHEX1,
1725                           tileHEX2,
1726                           tileHEX3,
1727                           tileHEX4,
1728                           tileHEX5,
1729                           tileHEX6,
1730                           tileRMin,
1731                           tileRMax,
1732                           tileRMinFine,
1733                           tileRMaxFine,
1734                           tileRingMin,
1735                           tileRingMax,
1736                           tileRingMinFine,
1737                           tileRingMaxFine);
1738   }
1739 }
1740 
1741 void HGCalGeomParameters::loadSpecParsTrapezoid(const cms::DDFilteredView& fv,
1742                                                 const cms::DDVectorsMap& vmap,
1743                                                 HGCalParameters& php,
1744                                                 const std::string& sdTag1) {
1745   for (auto const& it : vmap) {
1746     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1747       for (const auto& i : it.second)
1748         php.radiusMixBoundary_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1749     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1750       for (const auto& i : it.second)
1751         php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1752     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1753       for (const auto& i : it.second)
1754         php.layerCenter_.emplace_back(std::round(i));
1755     }
1756   }
1757 
1758   php.nPhiBinBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "NPhiBinBH"));
1759   php.layerFrontBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "LayerFrontBH"));
1760   php.rMinLayerBH_ = fv.get<std::vector<double> >(sdTag1, "RMinLayerBH");
1761   rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDD4hep);
1762   assert(php.nPhiBinBH_.size() > 1);
1763   php.nCellsFine_ = php.nPhiBinBH_[1];
1764   php.nCellsCoarse_ = php.nPhiBinBH_[0];
1765   assert(0 != php.nCellsFine_);
1766   assert(0 != php.nCellsCoarse_);
1767   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1768   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1769 
1770   php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1771   php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1772   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDD4hep);
1773   php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1774   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDD4hep);
1775 
1776   php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1777   php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1778   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDD4hep);
1779   php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1780   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDD4hep);
1781   unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1782   for (unsigned int k = 0; k < kmax; ++k)
1783     php.slopeTop_.emplace_back(0);
1784 
1785   const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1786   php.layerOffset_ = dummy2[0];
1787 
1788   loadSpecParsTrapezoid(php);
1789 
1790   // tile parameters from Katja's file
1791   if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette) ||
1792       (php.waferMaskMode_ == scintillatorFineCell)) {
1793     std::vector<int> tileIndx, tileProperty;
1794     std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4, tileHEX5, tileHEX6;
1795     std::vector<double> tileRMin, tileRMax, tileRMinFine, tileRMaxFine;
1796     std::vector<int> tileRingMin, tileRingMax, tileRingMinFine, tileRingMaxFine;
1797     std::vector<double> cassetteShift;
1798     for (auto const& it : vmap) {
1799       if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileIndex")) {
1800         for (const auto& i : it.second)
1801           tileIndx.emplace_back(std::round(i));
1802       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileProperty")) {
1803         for (const auto& i : it.second)
1804           tileProperty.emplace_back(std::round(i));
1805       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "NPhiLayer")) {
1806         for (const auto& i : it.second)
1807           php.nPhiLayer_.emplace_back(std::round(i));
1808       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX1")) {
1809         for (const auto& i : it.second)
1810           tileHEX1.emplace_back(std::round(i));
1811       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX2")) {
1812         for (const auto& i : it.second)
1813           tileHEX2.emplace_back(std::round(i));
1814       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX3")) {
1815         for (const auto& i : it.second)
1816           tileHEX3.emplace_back(std::round(i));
1817       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX4")) {
1818         for (const auto& i : it.second)
1819           tileHEX4.emplace_back(std::round(i));
1820       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin")) {
1821         for (const auto& i : it.second)
1822           tileRMin.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1823       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax")) {
1824         for (const auto& i : it.second)
1825           tileRMax.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1826       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin")) {
1827         for (const auto& i : it.second)
1828           tileRingMin.emplace_back(std::round(i));
1829       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax")) {
1830         for (const auto& i : it.second)
1831           tileRingMax.emplace_back(std::round(i));
1832       }
1833     }
1834     if (php.waferMaskMode_ == scintillatorFineCell) {
1835       php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1836       php.nphiFineCassette_ = php.nCellsFine_ / php.cassettes_;
1837       std::vector<double> rectract;
1838       for (auto const& it : vmap) {
1839         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX5")) {
1840           for (const auto& i : it.second)
1841             tileHEX5.emplace_back(std::round(i));
1842         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX6")) {
1843           for (const auto& i : it.second)
1844             tileHEX6.emplace_back(std::round(i));
1845         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin6")) {
1846           for (const auto& i : it.second)
1847             tileRMinFine.emplace_back(i);
1848           rescale(tileRMinFine, HGCalParameters::k_ScaleFromDDD);
1849         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax6")) {
1850           for (const auto& i : it.second)
1851             tileRMaxFine.emplace_back(i);
1852           rescale(tileRMaxFine, HGCalParameters::k_ScaleFromDDD);
1853         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin6")) {
1854           for (const auto& i : it.second)
1855             tileRingMinFine.emplace_back(std::round(i));
1856         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax6")) {
1857           for (const auto& i : it.second)
1858             tileRingMaxFine.emplace_back(std::round(i));
1859         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ScintRetract")) {
1860           for (const auto& i : it.second)
1861             rectract.emplace_back(i);
1862           rescale(rectract, HGCalParameters::k_ScaleFromDDD);
1863           for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1864             php.cassetteRetractTile_.emplace_back(rectract[k1]);
1865 #ifdef EDM_ML_DEBUG
1866             edm::LogVerbatim("HGCalGeom") << "cassetteRetractTile_[" << k1 << "] " << rectract[k1];
1867 #endif
1868           }
1869           int n = 2 * php.cassettes_ * (php.firstLayer_ - 1);
1870           for (int k1 = 0; k1 < n; ++k1)
1871             cassetteShift.emplace_back(0.);
1872 #ifdef EDM_ML_DEBUG
1873           edm::LogVerbatim("HGCalGeom") << "First Layer " << php.firstLayer_ << " Rett size " << rectract.size()
1874                                         << " N " << n;
1875 #endif
1876           double dphi = M_PI / php.cassettes_;
1877           for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1878             for (int k2 = 0; k2 < php.cassettes_; ++k2) {
1879               double phi = (2 * k2 + 1) * dphi;
1880               cassetteShift.emplace_back(rectract[k1] * cos(phi));
1881               cassetteShift.emplace_back(rectract[k1] * sin(phi));
1882             }
1883           }
1884         }
1885       }
1886     } else if (php.waferMaskMode_ == scintillatorCassette) {
1887       for (auto const& it : vmap) {
1888         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1889           for (const auto& i : it.second)
1890             cassetteShift.emplace_back(i);
1891         }
1892       }
1893     }
1894 
1895     rescale(cassetteShift, HGCalParameters::k_ScaleFromDD4hep);
1896     if (php.waferMaskMode_ == scintillatorFineCell)
1897       php.cassetteShiftTile_ = cassetteShift;
1898     else
1899       php.cassetteShift_ = cassetteShift;
1900     loadSpecParsTrapezoid(php,
1901                           tileIndx,
1902                           tileProperty,
1903                           tileHEX1,
1904                           tileHEX2,
1905                           tileHEX3,
1906                           tileHEX4,
1907                           tileHEX5,
1908                           tileHEX6,
1909                           tileRMin,
1910                           tileRMax,
1911                           tileRMinFine,
1912                           tileRMaxFine,
1913                           tileRingMin,
1914                           tileRingMax,
1915                           tileRingMinFine,
1916                           tileRingMaxFine);
1917   }
1918 }
1919 
1920 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php) {
1921 #ifdef EDM_ML_DEBUG
1922   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters:nCells " << php.nCellsFine_ << ":" << php.nCellsCoarse_
1923                                 << " cellSize: " << php.cellSize_[0] << ":" << php.cellSize_[1];
1924   for (unsigned int k = 0; k < php.layerFrontBH_.size(); ++k)
1925     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Type[" << k << "] Front Layer = " << php.layerFrontBH_[k]
1926                                   << " rMin = " << php.rMinLayerBH_[k];
1927   for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k) {
1928     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k]
1929                                   << " Nphi = " << php.scintCells(k + php.firstLayer_)
1930                                   << " dPhi = " << php.scintCellSize(k + php.firstLayer_);
1931   }
1932 
1933   for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1934     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1935                                   << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1936 
1937   for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1938     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1939                                   << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1940 
1941   std::ostringstream st1;
1942   for (unsigned int k = 0; k < php.zRanges_.size(); ++k)
1943     st1 << ":" << php.zRanges_[k];
1944   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary[" << php.zRanges_.size() << "] " << st1.str();
1945 
1946   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1947                                 << php.layerCenter_.size();
1948   for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1949     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1950 #endif
1951 }
1952 
1953 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php,
1954                                                 const std::vector<int>& tileIndx,
1955                                                 const std::vector<int>& tileProperty,
1956                                                 const std::vector<int>& tileHEX1,
1957                                                 const std::vector<int>& tileHEX2,
1958                                                 const std::vector<int>& tileHEX3,
1959                                                 const std::vector<int>& tileHEX4,
1960                                                 const std::vector<int>& tileHEX5,
1961                                                 const std::vector<int>& tileHEX6,
1962                                                 const std::vector<double>& tileRMin,
1963                                                 const std::vector<double>& tileRMax,
1964                                                 const std::vector<double>& tileRMinFine,
1965                                                 const std::vector<double>& tileRMaxFine,
1966                                                 const std::vector<int>& tileRingMin,
1967                                                 const std::vector<int>& tileRingMax,
1968                                                 const std::vector<int>& tileRingMinFine,
1969                                                 const std::vector<int>& tileRingMaxFine) {
1970   // tile parameters from Katja's file
1971   for (unsigned int k = 0; k < tileIndx.size(); ++k) {
1972     int hex5 = (k < tileHEX5.size()) ? tileHEX5[k] : 0;
1973     int hex6 = (k < tileHEX6.size()) ? tileHEX6[k] : 0;
1974     php.tileInfoMap_[tileIndx[k]] = HGCalParameters::tileInfo(HGCalProperty::tileType(tileProperty[k]),
1975                                                               HGCalProperty::tileSiPM(tileProperty[k]),
1976                                                               tileHEX1[k],
1977                                                               tileHEX2[k],
1978                                                               tileHEX3[k],
1979                                                               tileHEX4[k],
1980                                                               hex5,
1981                                                               hex6);
1982 #ifdef EDM_ML_DEBUG
1983     edm::LogVerbatim("HGCalGeom") << "Tile[" << k << ":" << tileIndx[k] << "] "
1984                                   << " Type " << HGCalProperty::tileType(tileProperty[k]) << " SiPM "
1985                                   << HGCalProperty::tileSiPM(tileProperty[k]) << " HEX " << std::hex << tileHEX1[k]
1986                                   << ":" << tileHEX2[k] << ":" << tileHEX3[k] << ":" << tileHEX4[k] << ":" << hex5
1987                                   << ":" << hex6 << std::dec;
1988 #endif
1989   }
1990 
1991   for (unsigned int k = 0; k < tileRMinFine.size(); ++k) {
1992     php.tileRingFineR_.emplace_back(tileRMinFine[k], tileRMaxFine[k]);
1993 #ifdef EDM_ML_DEBUG
1994     edm::LogVerbatim("HGCalGeom") << "TileRingFineR[" << k << "] " << tileRMinFine[k] << ":" << tileRMaxFine[k];
1995 #endif
1996   }
1997   for (unsigned int k = 0; k < tileRMin.size(); ++k) {
1998     php.tileRingR_.emplace_back(tileRMin[k], tileRMax[k]);
1999 #ifdef EDM_ML_DEBUG
2000     edm::LogVerbatim("HGCalGeom") << "TileRingR[" << k << "] " << tileRMin[k] << ":" << tileRMax[k];
2001 #endif
2002   }
2003 
2004   for (unsigned int k = 0; k < tileRingMinFine.size(); ++k) {
2005     php.tileRingFineRange_.emplace_back(tileRingMinFine[k], tileRingMaxFine[k]);
2006 #ifdef EDM_ML_DEBUG
2007     edm::LogVerbatim("HGCalGeom") << "TileRingFineRange[" << k << "] " << tileRingMinFine[k] << ":"
2008                                   << tileRingMaxFine[k];
2009 #endif
2010   }
2011   for (unsigned k = 0; k < tileRingMin.size(); ++k) {
2012     php.tileRingRange_.emplace_back(tileRingMin[k], tileRingMax[k]);
2013 #ifdef EDM_ML_DEBUG
2014     edm::LogVerbatim("HGCalGeom") << "TileRingRange[" << k << "] " << tileRingMin[k] << ":" << tileRingMax[k];
2015 #endif
2016   }
2017 #ifdef EDM_ML_DEBUG
2018   std::ostringstream st1, st2;
2019   st1 << "Scintillator CassetteShift for " << php.cassetteShiftTile_.size() << " tiles:";
2020   for (unsigned int k = 0; k < php.cassetteShiftTile_.size(); ++k)
2021     st1 << " " << php.cassetteShiftTile_[k];
2022   edm::LogVerbatim("HGCalGeom") << st1.str();
2023   st2 << "NPhi for scntillator " << php.nPhiLayer_.size() << " layers:";
2024   for (unsigned int k = 0; k < php.nPhiLayer_.size(); ++k)
2025     st2 << " " << php.nPhiLayer_[k];
2026   edm::LogVerbatim("HGCalGeom") << st2.str();
2027 #endif
2028 }
2029 
2030 void HGCalGeomParameters::loadWaferHexagon(HGCalParameters& php) {
2031   double waferW(HGCalParameters::k_ScaleFromDDD * waferSize_), rmin(HGCalParameters::k_ScaleFromDDD * php.waferR_);
2032   double rin(php.rLimit_[0]), rout(php.rLimit_[1]), rMaxFine(php.boundR_[1]);
2033 #ifdef EDM_ML_DEBUG
2034   edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << rmin << " R Limits: " << rin << ":" << rout
2035                                 << " Fine " << rMaxFine;
2036 #endif
2037   // Clear the vectors
2038   php.waferCopy_.clear();
2039   php.waferTypeL_.clear();
2040   php.waferTypeT_.clear();
2041   php.waferPosX_.clear();
2042   php.waferPosY_.clear();
2043   double dx = 0.5 * waferW;
2044   double dy = 3.0 * dx * tan(30._deg);
2045   double rr = 2.0 * dx * tan(30._deg);
2046   int ncol = static_cast<int>(2.0 * rout / waferW) + 1;
2047   int nrow = static_cast<int>(rout / (waferW * tan(30._deg))) + 1;
2048   int ns2 = (2 * ncol + 1) * (2 * nrow + 1) * php.layer_.size();
2049   int incm(0), inrm(0);
2050   HGCalParameters::layer_map copiesInLayers(php.layer_.size() + 1);
2051   HGCalParameters::waferT_map waferTypes(ns2 + 1);
2052 #ifdef EDM_ML_DEBUG
2053   int kount(0), ntot(0);
2054   edm::LogVerbatim("HGCalGeom") << "Row " << nrow << " Column " << ncol;
2055 #endif
2056   for (int nr = -nrow; nr <= nrow; ++nr) {
2057     int inr = (nr >= 0) ? nr : -nr;
2058     for (int nc = -ncol; nc <= ncol; ++nc) {
2059       int inc = (nc >= 0) ? nc : -nc;
2060       if (inr % 2 == inc % 2) {
2061         double xpos = nc * dx;
2062         double ypos = nr * dy;
2063         std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, rin, rout, true);
2064         double rpos = std::sqrt(xpos * xpos + ypos * ypos);
2065         int typet = (rpos < rMaxFine) ? 1 : 2;
2066         int typel(3);
2067         for (int k = 1; k < 4; ++k) {
2068           if ((rpos + rmin) <= php.boundR_[k]) {
2069             typel = k;
2070             break;
2071           }
2072         }
2073 #ifdef EDM_ML_DEBUG
2074         ++ntot;
2075 #endif
2076         if (corner.first > 0) {
2077           int copy = HGCalTypes::packTypeUV(typel, nc, nr);
2078           if (inc > incm)
2079             incm = inc;
2080           if (inr > inrm)
2081             inrm = inr;
2082 #ifdef EDM_ML_DEBUG
2083           kount++;
2084           edm::LogVerbatim("HGCalGeom") << kount << ":" << ntot << " Copy " << copy << " Type " << typel << ":" << typet
2085                                         << " Location " << corner.first << " Position " << xpos << ":" << ypos
2086                                         << " Layers " << php.layer_.size();
2087 #endif
2088           php.waferCopy_.emplace_back(copy);
2089           php.waferTypeL_.emplace_back(typel);
2090           php.waferTypeT_.emplace_back(typet);
2091           php.waferPosX_.emplace_back(xpos);
2092           php.waferPosY_.emplace_back(ypos);
2093           for (unsigned int il = 0; il < php.layer_.size(); ++il) {
2094             std::pair<int, int> corner =
2095                 HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, php.rMinLayHex_[il], php.rMaxLayHex_[il], true);
2096             if (corner.first > 0) {
2097               auto cpy = copiesInLayers[php.layer_[il]].find(copy);
2098               if (cpy == copiesInLayers[php.layer_[il]].end())
2099                 copiesInLayers[php.layer_[il]][copy] =
2100                     ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ? php.waferCopy_.size() : -1);
2101             }
2102             if ((corner.first > 0) && (corner.first < static_cast<int>(HGCalParameters::k_CornerSize))) {
2103               int wl = HGCalWaferIndex::waferIndex(php.layer_[il], copy, 0, true);
2104               waferTypes[wl] = corner;
2105             }
2106           }
2107         }
2108       }
2109     }
2110   }
2111   php.copiesInLayers_ = copiesInLayers;
2112   php.waferTypes_ = waferTypes;
2113   php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2114   php.waferUVMax_ = 0;
2115 #ifdef EDM_ML_DEBUG
2116   edm::LogVerbatim("HGCalGeom") << "HGCalWaferHexagon: # of columns " << incm << " # of rows " << inrm << " and "
2117                                 << kount << ":" << ntot << " wafers; R " << rin << ":" << rout;
2118   edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
2119   for (unsigned int k = 0; k < copiesInLayers.size(); ++k) {
2120     const auto& theModules = copiesInLayers[k];
2121     edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
2122     int k2(0);
2123     for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
2124       edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
2125     }
2126   }
2127 #endif
2128 }
2129 
2130 void HGCalGeomParameters::loadWaferHexagon8(HGCalParameters& php) {
2131   double waferW(php.waferSize_);
2132   double waferS(php.sensorSeparation_);
2133   auto wType = std::make_unique<HGCalWaferType>(php.radius100to200_,
2134                                                 php.radius200to300_,
2135                                                 HGCalParameters::k_ScaleToDDD * (waferW + waferS),
2136                                                 HGCalParameters::k_ScaleToDDD * php.zMinForRad_,
2137                                                 php.choiceType_,
2138                                                 php.nCornerCut_,
2139                                                 php.fracAreaMin_);
2140 
2141   double rout(php.rLimit_[1]);
2142 #ifdef EDM_ML_DEBUG
2143   edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << waferS << " R Max: " << rout;
2144 #endif
2145   // Clear the vectors
2146   php.waferCopy_.clear();
2147   php.waferTypeL_.clear();
2148   php.waferTypeT_.clear();
2149   php.waferPosX_.clear();
2150   php.waferPosY_.clear();
2151   double r = 0.5 * (waferW + waferS);
2152   double R = 2.0 * r / sqrt3_;
2153   double dy = 0.75 * R;
2154   double r1 = 0.5 * waferW;
2155   double R1 = 2.0 * r1 / sqrt3_;
2156   int N = (r == 0) ? 3 : (static_cast<int>(0.5 * rout / r) + 3);
2157   int ns1 = (2 * N + 1) * (2 * N + 1);
2158   int ns2 = ns1 * php.zLayerHex_.size();
2159 #ifdef EDM_ML_DEBUG
2160   edm::LogVerbatim("HGCalGeom") << "wafer " << waferW << ":" << waferS << " r " << r << " dy " << dy << " N " << N
2161                                 << " sizes " << ns1 << ":" << ns2;
2162   std::vector<int> indtypes(ns1 + 1);
2163   indtypes.clear();
2164 #endif
2165   HGCalParameters::wafer_map wafersInLayers(ns1 + 1);
2166   HGCalParameters::wafer_map typesInLayers(ns2 + 1);
2167   HGCalParameters::waferT_map waferTypes(ns2 + 1);
2168   int ipos(0), lpos(0), uvmax(0), nwarn(0);
2169   std::vector<int> uvmx(php.zLayerHex_.size(), 0);
2170   for (int v = -N; v <= N; ++v) {
2171     for (int u = -N; u <= N; ++u) {
2172       int nr = 2 * v;
2173       int nc = -2 * u + v;
2174       double xpos = nc * r;
2175       double ypos = nr * dy;
2176       int indx = HGCalWaferIndex::waferIndex(0, u, v);
2177       php.waferCopy_.emplace_back(indx);
2178       php.waferPosX_.emplace_back(xpos);
2179       php.waferPosY_.emplace_back(ypos);
2180       wafersInLayers[indx] = ipos;
2181       ++ipos;
2182       std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, r1, R1, 0, rout, false);
2183       if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2184           ((corner.first > 0) && php.defineFull_)) {
2185         uvmax = std::max(uvmax, std::max(std::abs(u), std::abs(v)));
2186       }
2187       for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
2188         int copy = i + php.layerOffset_;
2189         std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], (waferW + waferS));
2190         int lay = php.layer_[php.layerIndex_[i]];
2191         double xpos0 = xpos + xyoff.first;
2192         double ypos0 = ypos + xyoff.second;
2193         double zpos = php.zLayerHex_[i];
2194         int kndx = HGCalWaferIndex::waferIndex(lay, u, v);
2195         int type(-1);
2196         if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2197             (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) || (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell) ||
2198             (php.mode_ == HGCalGeometryMode::Hexagon8FineCell))
2199           type = wType->getType(kndx, php.waferInfoMap_);
2200         if (type < 0)
2201           type = wType->getType(HGCalParameters::k_ScaleToDDD * xpos0,
2202                                 HGCalParameters::k_ScaleToDDD * ypos0,
2203                                 HGCalParameters::k_ScaleToDDD * zpos);
2204         php.waferTypeL_.emplace_back(type);
2205         typesInLayers[kndx] = lpos;
2206         ++lpos;
2207 #ifdef EDM_ML_DEBUG
2208         indtypes.emplace_back(kndx);
2209 #endif
2210         std::pair<int, int> corner =
2211             HGCalGeomTools::waferCorner(xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], false);
2212 #ifdef EDM_ML_DEBUG
2213         if (((corner.first == 0) && std::abs(u) < 5 && std::abs(v) < 5) || (std::abs(u) < 2 && std::abs(v) < 2)) {
2214           edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " R " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2215                                         << " u " << u << " v " << v << " with " << corner.first << " corners";
2216         }
2217 #endif
2218         if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2219             ((corner.first > 0) && php.defineFull_)) {
2220           uvmx[i] = std::max(uvmx[i], std::max(std::abs(u), std::abs(v)));
2221         }
2222         if ((corner.first < static_cast<int>(HGCalParameters::k_CornerSize)) && (corner.first > 0)) {
2223 #ifdef EDM_ML_DEBUG
2224           edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2225                                         << corner.first << ":" << corner.second;
2226 #endif
2227           int wl = HGCalWaferIndex::waferIndex(lay, u, v);
2228           if (php.waferMaskMode_ > 0) {
2229             bool v17OrLess = (php.mode_ < HGCalGeometryMode::Hexagon8CalibCell);
2230             std::pair<int, int> corner0 = HGCalWaferMask::getTypeMode(
2231                 xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], type, php.waferMaskMode_, v17OrLess);
2232             if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2233                 (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2234                 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell) ||
2235                 (php.mode_ == HGCalGeometryMode::Hexagon8FineCell)) {
2236               auto itr = php.waferInfoMap_.find(wl);
2237               if (itr != php.waferInfoMap_.end()) {
2238                 int part = (itr->second).part;
2239                 int orient = (itr->second).orient;
2240                 bool ok = ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2241                            (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell) ||
2242                            (php.mode_ == HGCalGeometryMode::Hexagon8FineCell))
2243                               ? true
2244                               : HGCalWaferMask::goodTypeMode(
2245                                     xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], part, orient, false);
2246                 if (ok)
2247                   corner0 = std::make_pair(part, (HGCalTypes::k_OffsetRotation + orient));
2248 #ifdef EDM_ML_DEBUG
2249                 edm::LogVerbatim("HGCalGeom")
2250                     << "Layer:u:v " << i << ":" << lay << ":" << u << ":" << v << " Part " << corner0.first << ":"
2251                     << part << " Orient " << corner0.second << ":" << orient << " Position " << xpos0 << ":" << ypos0
2252                     << " delta " << r1 << ":" << R1 << " Limit " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2253                     << " Compatibiliety Flag " << ok;
2254 #endif
2255                 if (!ok)
2256                   ++nwarn;
2257               }
2258             }
2259             waferTypes[wl] = corner0;
2260 #ifdef EDM_ML_DEBUG
2261             edm::LogVerbatim("HGCalGeom")
2262                 << "Layer " << lay << " u|v " << u << ":" << v << " Index " << std::hex << wl << std::dec << " pos "
2263                 << xpos0 << ":" << ypos0 << " R " << r1 << ":" << R1 << " Range " << php.rMinLayHex_[i] << ":"
2264                 << php.rMaxLayHex_[i] << type << ":" << php.waferMaskMode_ << " corner " << corner.first << ":"
2265                 << corner.second << " croner0 " << corner0.first << ":" << corner0.second;
2266 #endif
2267           } else {
2268             waferTypes[wl] = corner;
2269 #ifdef EDM_ML_DEBUG
2270             edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2271                                           << corner.first << ":" << corner.second;
2272 #endif
2273           }
2274         }
2275       }
2276     }
2277   }
2278   if (nwarn > 0)
2279     edm::LogWarning("HGCalGeom") << "HGCalGeomParameters::loadWafer8: there are " << nwarn
2280                                  << " wafers with non-matching partial- orientation types";
2281   php.waferUVMax_ = uvmax;
2282   php.waferUVMaxLayer_ = uvmx;
2283   php.wafersInLayers_ = wafersInLayers;
2284   php.typesInLayers_ = typesInLayers;
2285   php.waferTypes_ = waferTypes;
2286   php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2287   HGCalParameters::hgtrap mytr;
2288   mytr.lay = 1;
2289   mytr.bl = php.waferR_;
2290   mytr.tl = php.waferR_;
2291   mytr.h = php.waferR_;
2292   mytr.alpha = 0.0;
2293   mytr.cellSize = HGCalParameters::k_ScaleToDDD * php.waferSize_;
2294   for (auto const& dz : php.cellThickness_) {
2295     mytr.dz = 0.5 * HGCalParameters::k_ScaleToDDD * dz;
2296     php.fillModule(mytr, false);
2297   }
2298   for (unsigned k = 0; k < php.cellThickness_.size(); ++k) {
2299     HGCalParameters::hgtrap mytr = php.getModule(k, false);
2300     mytr.bl *= HGCalParameters::k_ScaleFromDDD;
2301     mytr.tl *= HGCalParameters::k_ScaleFromDDD;
2302     mytr.h *= HGCalParameters::k_ScaleFromDDD;
2303     mytr.dz *= HGCalParameters::k_ScaleFromDDD;
2304     mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2305     php.fillModule(mytr, true);
2306   }
2307 #ifdef EDM_ML_DEBUG
2308   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Total of " << php.waferCopy_.size() << " wafers";
2309   for (unsigned int k = 0; k < php.waferCopy_.size(); ++k) {
2310     int id = php.waferCopy_[k];
2311     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << std::hex << id << std::dec << ":"
2312                                   << HGCalWaferIndex::waferLayer(id) << ":" << HGCalWaferIndex::waferU(id) << ":"
2313                                   << HGCalWaferIndex::waferV(id) << " x " << php.waferPosX_[k] << " y "
2314                                   << php.waferPosY_[k] << " index " << php.wafersInLayers_[id];
2315   }
2316   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Total of " << php.waferTypeL_.size() << " wafer types";
2317   for (unsigned int k = 0; k < php.waferTypeL_.size(); ++k) {
2318     int id = indtypes[k];
2319     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.typesInLayers_[id] << ":" << php.waferTypeL_[k] << " ID "
2320                                   << std::hex << id << std::dec << ":" << HGCalWaferIndex::waferLayer(id) << ":"
2321                                   << HGCalWaferIndex::waferU(id) << ":" << HGCalWaferIndex::waferV(id);
2322   }
2323 #endif
2324 
2325   //Wafer offset
2326   php.xLayerHex_.clear();
2327   php.yLayerHex_.clear();
2328   double waferSize = php.waferSize_ + php.sensorSeparation_;
2329 #ifdef EDM_ML_DEBUG
2330   edm::LogVerbatim("HGCalGeom") << "WaferSize " << waferSize;
2331 #endif
2332   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2333     int copy = k + php.layerOffset_;
2334     std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], waferSize);
2335     php.xLayerHex_.emplace_back(xyoff.first);
2336     php.yLayerHex_.emplace_back(xyoff.second);
2337 #ifdef EDM_ML_DEBUG
2338     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Off " << copy << ":" << php.layerCenter_[copy] << " Shift "
2339                                   << xyoff.first << ":" << xyoff.second;
2340 #endif
2341   }
2342 }
2343 
2344 void HGCalGeomParameters::loadCellParsHexagon(const DDCompactView* cpv, HGCalParameters& php) {
2345   // Special parameters for cell parameters
2346   std::string attribute = "OnlyForHGCalNumbering";
2347   DDSpecificsHasNamedValueFilter filter1{attribute};
2348   DDFilteredView fv1(*cpv, filter1);
2349   bool ok = fv1.firstChild();
2350 
2351   if (ok) {
2352     php.cellFine_ = dbl_to_int(cpv->vector("waferFine"));
2353     php.cellCoarse_ = dbl_to_int(cpv->vector("waferCoarse"));
2354   }
2355 
2356   loadCellParsHexagon(php);
2357 }
2358 
2359 void HGCalGeomParameters::loadCellParsHexagon(const cms::DDVectorsMap& vmap, HGCalParameters& php) {
2360   for (auto const& it : vmap) {
2361     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferFine")) {
2362       for (const auto& i : it.second)
2363         php.cellFine_.emplace_back(std::round(i));
2364     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferCoarse")) {
2365       for (const auto& i : it.second)
2366         php.cellCoarse_.emplace_back(std::round(i));
2367     }
2368   }
2369 
2370   loadCellParsHexagon(php);
2371 }
2372 
2373 void HGCalGeomParameters::loadCellParsHexagon(const HGCalParameters& php) {
2374 #ifdef EDM_ML_DEBUG
2375   edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellFine_.size() << " rows for fine cells";
2376   for (unsigned int k = 0; k < php.cellFine_.size(); ++k)
2377     edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellFine_[k];
2378   edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellCoarse_.size() << " rows for coarse cells";
2379   for (unsigned int k = 0; k < php.cellCoarse_.size(); ++k)
2380     edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellCoarse_[k];
2381 #endif
2382 }
2383 
2384 void HGCalGeomParameters::loadCellTrapezoid(HGCalParameters& php) {
2385   php.xLayerHex_.resize(php.zLayerHex_.size(), 0);
2386   php.yLayerHex_.resize(php.zLayerHex_.size(), 0);
2387 #ifdef EDM_ML_DEBUG
2388   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: x|y|zLayerHex in array of size " << php.zLayerHex_.size();
2389   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k)
2390     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Shift " << php.xLayerHex_[k] << ":" << php.yLayerHex_[k] << ":"
2391                                   << php.zLayerHex_[k];
2392 #endif
2393   // Find the radius of each eta-partitions
2394   if ((php.mode_ == HGCalGeometryMode::TrapezoidFile) || (php.mode_ == HGCalGeometryMode::TrapezoidModule) ||
2395       (php.mode_ == HGCalGeometryMode::TrapezoidCassette) || (php.mode_ == HGCalGeometryMode::TrapezoidFineCell)) {
2396     //Ring radii for each partition
2397 #ifdef EDM_ML_DEBUG
2398     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Mode " << php.mode_ << ":"
2399                                   << HGCalGeometryMode::TrapezoidFineCell << " Sizes " << php.tileRingFineR_.size()
2400                                   << ":" << php.tileRingR_.size();
2401 #endif
2402     for (unsigned int k = 0; k < 2; ++k) {
2403       bool fine = ((k == 1) && (php.mode_ == HGCalGeometryMode::TrapezoidFineCell));
2404       unsigned int sizeR = (fine) ? php.tileRingFineR_.size() : php.tileRingR_.size();
2405       for (unsigned int kk = 0; kk < sizeR; ++kk) {
2406         if (fine)
2407           php.radiusLayer_[k].emplace_back(php.tileRingFineR_[kk].first);
2408         else
2409           php.radiusLayer_[k].emplace_back(php.tileRingR_[kk].first);
2410 #ifdef EDM_ML_DEBUG
2411         double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2412                               : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2413         double rv = php.radiusLayer_[k].back();
2414         double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2415         edm::LogVerbatim("HGCalGeom") << "New [" << kk << "] new R = " << rv << " Eta = " << eta;
2416 #endif
2417       }
2418       if (fine)
2419         php.radiusLayer_[k].emplace_back(php.tileRingFineR_[php.tileRingFineR_.size() - 1].second);
2420       else
2421         php.radiusLayer_[k].emplace_back(php.tileRingR_[php.tileRingR_.size() - 1].second);
2422     }
2423     // Minimum and maximum radius index for each layer
2424 #ifdef EDM_ML_DEBUG
2425     edm::LogVerbatim("HGCalGeom") << "Tile Ring Range size " << php.tileRingFineRange_.size() << ":"
2426                                   << php.tileRingRange_.size() << ":" << php.nPhiBinBH_.size() << ":"
2427                                   << php.zLayerHex_.size() << ":" << php.nPhiLayer_.size();
2428 #endif
2429     unsigned int k1(0), k2(0);
2430     for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2431       if (!php.tileRingFineRange_.empty()) {
2432         unsigned int k0 = (k1 < (php.tileRingFineRange_.size() - 1)) ? k1 : (php.tileRingFineRange_.size() - 1);
2433 #ifdef EDM_ML_DEBUG
2434         edm::LogVerbatim("HGCalGeom") << "Layer " << k << " Fine " << k1 << ":" << k0 << ":"
2435                                       << php.tileRingFineRange_.size();
2436 #endif
2437         php.iradMinBHFine_.emplace_back(1 + php.tileRingFineRange_[k0].first);
2438         php.iradMaxBHFine_.emplace_back(1 + php.tileRingFineRange_[k0].second);
2439       }
2440       if (!php.tileRingRange_.empty()) {
2441         unsigned int k0 = (k2 < (php.tileRingRange_.size() - 1)) ? k2 : (php.tileRingRange_.size() - 1);
2442 #ifdef EDM_ML_DEBUG
2443         edm::LogVerbatim("HGCalGeom") << "Layer " << k << " Coarse " << k2 << ":" << k0 << ":"
2444                                       << php.tileRingRange_.size();
2445 #endif
2446         php.iradMinBH_.emplace_back(1 + php.tileRingRange_[k0].first);
2447         php.iradMaxBH_.emplace_back(1 + php.tileRingRange_[k0].second);
2448       }
2449       if (php.nPhiLayer_[k] > 288) {
2450         ++k1;
2451       } else {
2452         ++k2;
2453       }
2454 #ifdef EDM_ML_DEBUG
2455       int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2456       if (php.nPhiLayer_[k] > 288)
2457         edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge "
2458                                       << php.iradMinBHFine_.back() << " Top edge " << php.iradMaxBHFine_.back();
2459       else
2460         edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge " << php.iradMinBH_.back()
2461                                       << " Top edge " << php.iradMaxBH_.back();
2462 #endif
2463     }
2464   } else {
2465     //Ring radii for each partition
2466     for (unsigned int k = 0; k < 2; ++k) {
2467       double rmax = ((k == 0) ? (php.rMaxLayHex_[php.layerFrontBH_[1] - php.firstLayer_] - 1)
2468                               : (php.rMaxLayHex_[php.rMaxLayHex_.size() - 1]));
2469       double rv = php.rMinLayerBH_[k];
2470       double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2471                             : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2472       php.radiusLayer_[k].emplace_back(rv);
2473 #ifdef EDM_ML_DEBUG
2474       double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2475       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] rmax " << rmax << " Z = " << zv
2476                                     << " dEta = " << php.cellSize_[k] << "\nOld[0] new R = " << rv << " Eta = " << eta;
2477       int kount(1);
2478 #endif
2479       while (rv < rmax) {
2480         double eta = -(php.cellSize_[k] + std::log(std::tan(0.5 * std::atan(rv / zv))));
2481         rv = zv * std::tan(2.0 * std::atan(std::exp(-eta)));
2482         php.radiusLayer_[k].emplace_back(rv);
2483 #ifdef EDM_ML_DEBUG
2484         edm::LogVerbatim("HGCalGeom") << "Old [" << kount << "] new R = " << rv << " Eta = " << eta;
2485         ++kount;
2486 #endif
2487       }
2488     }
2489     // Find minimum and maximum radius index for each layer
2490     for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2491       int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2492       std::vector<double>::iterator low, high;
2493       low = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMinLayHex_[k]);
2494 #ifdef EDM_ML_DEBUG
2495       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RLow = " << php.rMinLayHex_[k] << " pos "
2496                                     << static_cast<int>(low - php.radiusLayer_[kk].begin());
2497 #endif
2498       if (low == php.radiusLayer_[kk].begin())
2499         ++low;
2500       int irlow = static_cast<int>(low - php.radiusLayer_[kk].begin());
2501       double drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2502 #ifdef EDM_ML_DEBUG
2503       edm::LogVerbatim("HGCalGeom") << "irlow " << irlow << " dr " << drlow << " min " << php.minTileSize_;
2504 #endif
2505       if (drlow < php.minTileSize_) {
2506         ++irlow;
2507 #ifdef EDM_ML_DEBUG
2508         drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2509         edm::LogVerbatim("HGCalGeom") << "Modified irlow " << irlow << " dr " << drlow;
2510 #endif
2511       }
2512       high = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMaxLayHex_[k]);
2513 #ifdef EDM_ML_DEBUG
2514       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RHigh = " << php.rMaxLayHex_[k] << " pos "
2515                                     << static_cast<int>(high - php.radiusLayer_[kk].begin());
2516 #endif
2517       if (high == php.radiusLayer_[kk].end())
2518         --high;
2519       int irhigh = static_cast<int>(high - php.radiusLayer_[kk].begin());
2520       double drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2521 #ifdef EDM_ML_DEBUG
2522       edm::LogVerbatim("HGCalGeom") << "irhigh " << irhigh << " dr " << drhigh << " min " << php.minTileSize_;
2523 #endif
2524       if (drhigh < php.minTileSize_) {
2525         --irhigh;
2526 #ifdef EDM_ML_DEBUG
2527         drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2528         edm::LogVerbatim("HGCalGeom") << "Modified irhigh " << irhigh << " dr " << drhigh;
2529 #endif
2530       }
2531       php.iradMinBHFine_.emplace_back(irlow);
2532       php.iradMaxBHFine_.emplace_back(irhigh);
2533       php.iradMinBH_.emplace_back(irlow);
2534       php.iradMaxBH_.emplace_back(irhigh);
2535 #ifdef EDM_ML_DEBUG
2536       edm::LogVerbatim("HGCalGeom") << "Old Layer " << k << " Type " << kk << " Low edge " << irlow << ":" << drlow
2537                                     << " Top edge " << irhigh << ":" << drhigh;
2538 #endif
2539     }
2540   }
2541 #ifdef EDM_ML_DEBUG
2542   for (unsigned int k = 0; k < 2; ++k) {
2543     edm::LogVerbatim("HGCalGeom") << "Type " << k << " with " << php.radiusLayer_[k].size() << " radii";
2544     for (unsigned int kk = 0; kk < php.radiusLayer_[k].size(); ++kk)
2545       edm::LogVerbatim("HGCalGeom") << "Ring[" << k << "][" << kk << "] " << php.radiusLayer_[k][kk];
2546   }
2547 #endif
2548 
2549   // Now define the volumes
2550   int im(0);
2551   php.tileUVMax_ = php.tileUVMaxFine_ = 00;
2552   HGCalParameters::hgtrap mytr;
2553   mytr.alpha = 0.0;
2554   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2555     if (php.nPhiLayer_[k] > 288) {
2556       if (php.iradMaxBHFine_[k] > php.tileUVMaxFine_)
2557         php.tileUVMaxFine_ = php.iradMaxBHFine_[k];
2558     } else {
2559       if (php.iradMaxBH_[k] > php.tileUVMax_)
2560         php.tileUVMax_ = php.iradMaxBH_[k];
2561     }
2562     int kk = (php.nPhiLayer_[k] > 288) ? 1 : 0;
2563     int irm = php.radiusLayer_[kk].size() - 1;
2564     int irmin = (php.nPhiLayer_[k] > 288) ? php.iradMinBHFine_[k] : php.iradMinBH_[k];
2565     int irmax = (php.nPhiLayer_[k] > 288) ? php.iradMaxBHFine_[k] : php.iradMaxBH_[k];
2566 #ifdef EDM_ML_DEBUG
2567     double rmin = php.radiusLayer_[kk][std::max((irmin - 1), 0)];
2568     double rmax = php.radiusLayer_[kk][std::min(irmax, irm)];
2569     edm::LogVerbatim("HGCalGeom") << "Layer " << php.firstLayer_ + k << ":" << kk << " Radius range " << irmin << ":"
2570                                   << irmax << ":" << rmin << ":" << rmax << " Size " << php.radiusLayer_[kk].size();
2571 #endif
2572     mytr.lay = php.firstLayer_ + k;
2573     for (int irad = irmin; irad <= irmax; ++irad) {
2574       double rmin = php.radiusLayer_[kk][std::max((irad - 1), 0)];
2575       double rmax = php.radiusLayer_[kk][std::min(irad, irm)];
2576 #ifdef EDM_ML_DEBUG
2577       edm::LogVerbatim("HGCalGeom") << "irad " << irad << ":" << php.radiusLayer_[kk].size() << " R " << rmin << ":"
2578                                     << rmax;
2579 #endif
2580       mytr.bl = 0.5 * rmin * php.scintCellSize(mytr.lay);
2581       mytr.tl = 0.5 * rmax * php.scintCellSize(mytr.lay);
2582       mytr.h = 0.5 * (rmax - rmin);
2583       mytr.dz = 0.5 * php.waferThick_;
2584       mytr.cellSize = 0.5 * (rmax + rmin) * php.scintCellSize(mytr.lay);
2585       php.fillModule(mytr, true);
2586       mytr.bl *= HGCalParameters::k_ScaleToDDD;
2587       mytr.tl *= HGCalParameters::k_ScaleToDDD;
2588       mytr.h *= HGCalParameters::k_ScaleToDDD;
2589       mytr.dz *= HGCalParameters::k_ScaleToDDD;
2590       mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2591       php.fillModule(mytr, false);
2592       if (irad == irmin)
2593         php.firstModule_.emplace_back(im);
2594       ++im;
2595       if (irad == irmax)
2596         php.lastModule_.emplace_back(im);
2597     }
2598   }
2599   php.nSectors_ = std::max(php.tileUVMax_, php.tileUVMaxFine_);
2600 #ifdef EDM_ML_DEBUG
2601   edm::LogVerbatim("HGCalGeom") << "Maximum radius index " << php.tileUVMax_ << ":" << php.tileUVMaxFine_;
2602   for (unsigned int k = 0; k < php.firstModule_.size(); ++k)
2603     edm::LogVerbatim("HGCalGeom") << "Layer " << k + php.firstLayer_ << " Modules " << php.firstModule_[k] << ":"
2604                                   << php.lastModule_[k];
2605 
2606 #endif
2607 }
2608 
2609 std::vector<double> HGCalGeomParameters::getDDDArray(const std::string& str, const DDsvalues_type& sv, const int nmin) {
2610   DDValue value(str);
2611   if (DDfetch(&sv, value)) {
2612     const std::vector<double>& fvec = value.doubles();
2613     int nval = fvec.size();
2614     if (nmin > 0) {
2615       if (nval < nmin) {
2616         throw cms::Exception("DDException")
2617             << "HGCalGeomParameters:  # of " << str << " bins " << nval << " < " << nmin << " ==> illegal";
2618       }
2619     } else {
2620       if (nval < 1 && nmin == 0) {
2621         throw cms::Exception("DDException")
2622             << "HGCalGeomParameters: # of " << str << " bins " << nval << " < 1 ==> illegal"
2623             << " (nmin=" << nmin << ")";
2624       }
2625     }
2626     return fvec;
2627   } else {
2628     if (nmin >= 0) {
2629       throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str;
2630     }
2631     std::vector<double> fvec;
2632     return fvec;
2633   }
2634 }
2635 
2636 void HGCalGeomParameters::rescale(std::vector<double>& v, const double s) {
2637   std::for_each(v.begin(), v.end(), [s](double& n) { n *= s; });
2638 }
2639 
2640 std::pair<double, double> HGCalGeomParameters::cellPosition(
2641     const std::vector<HGCalGeomParameters::cellParameters>& wafers,
2642     std::vector<HGCalGeomParameters::cellParameters>::const_iterator& itrf,
2643     int wafer,
2644     double xx,
2645     double yy) {
2646   if (itrf == wafers.end()) {
2647     for (std::vector<HGCalGeomParameters::cellParameters>::const_iterator itr = wafers.begin(); itr != wafers.end();
2648          ++itr) {
2649       if (itr->wafer == wafer) {
2650         itrf = itr;
2651         break;
2652       }
2653     }
2654   }
2655   double dx(0), dy(0);
2656   if (itrf != wafers.end()) {
2657     dx = (xx - itrf->xyz.x());
2658     if (std::abs(dx) < tolerance)
2659       dx = 0;
2660     dy = (yy - itrf->xyz.y());
2661     if (std::abs(dy) < tolerance)
2662       dy = 0;
2663   }
2664   return std::make_pair(dx, dy);
2665 }
2666 
2667 void HGCalGeomParameters::resetZero(std::vector<double>& v) {
2668   for (auto& n : v) {
2669     if (std::abs(n) < tolmin)
2670       n = 0;
2671   }
2672 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!