Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-08 23:18:59 

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 #ifdef EDM_ML_DEBUG
0885     ++ntot2;
0886 #endif
0887     std::vector<int> copy = fv2.copyNumbers();
0888     int nsiz = static_cast<int>(copy.size());
0889     if (levelTop < nsiz) {
0890       int lay = copy[levelTop];
0891       int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
0892       if (zside != 1)
0893         zside = -1;
0894       const DDSolid& sol = fv2.logicalPart().solid();
0895 #ifdef EDM_ML_DEBUG
0896       std::ostringstream st2;
0897       st2 << "Name1 " << sol.name() << " shape " << sol.shape() << " LTop " << levelTop << ":" << lay << " ZSide "
0898           << zside << ":" << php.levelZSide_ << " # of levels " << nsiz;
0899       for (const auto& c : copy)
0900         st2 << ":" << c;
0901       edm::LogVerbatim("HGCalGeom") << st2.str();
0902 #endif
0903       if (lay == 0) {
0904         throw cms::Exception("DDException")
0905             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0906       } else if (sol.shape() == DDSolidShape::ddtubs) {
0907         if (zvals.find(std::make_pair(lay, zside)) != zvals.end()) {
0908           if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0909             php.layer_.emplace_back(lay);
0910           auto itr = layers.find(lay);
0911           if (itr == layers.end()) {
0912             const DDTubs& tube = static_cast<DDTubs>(sol);
0913             double rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
0914             double rout = (php.firstMixedLayer_ > 0 && lay >= php.firstMixedLayer_)
0915                               ? php.radiusMixBoundary_[lay - php.firstMixedLayer_]
0916                               : HGCalParameters::k_ScaleFromDDD * tube.rOut();
0917             double zp = zvals[std::make_pair(lay, 1)];
0918             HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
0919             layers[lay] = laypar;
0920 #ifdef EDM_ML_DEBUG
0921             std::ostringstream st3;
0922             st3 << "Name1 " << fv2.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
0923                 << nsiz;
0924             for (const auto& c : copy)
0925               st3 << ":" << c;
0926             st3 << " R " << rin << ":" << rout;
0927             edm::LogVerbatim("HGCalGeom") << st3.str();
0928 #endif
0929           }
0930 
0931           if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
0932             DD3Vector x, y, z;
0933             fv2.rotation().GetComponents(x, y, z);
0934             const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0935             const CLHEP::HepRotation hr(rotation);
0936             double xx = ((std::abs(fv2.translation().X()) < tolerance)
0937                              ? 0
0938                              : HGCalParameters::k_ScaleFromDDD * fv2.translation().X());
0939             double yy = ((std::abs(fv2.translation().Y()) < tolerance)
0940                              ? 0
0941                              : HGCalParameters::k_ScaleFromDDD * fv2.translation().Y());
0942             const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
0943             HGCalParameters::hgtrform mytrf;
0944             mytrf.zp = zside;
0945             mytrf.lay = lay;
0946             mytrf.sec = 0;
0947             mytrf.subsec = 0;
0948             mytrf.h3v = h3v;
0949             mytrf.hr = hr;
0950             trforms[std::make_pair(lay, zside)] = mytrf;
0951 #ifdef EDM_ML_DEBUG
0952             edm::LogVerbatim("HGCalGeom") << "Translation " << h3v;
0953 #endif
0954           }
0955         }
0956       }
0957     }
0958     dodet = fv2.next();
0959   }
0960 #ifdef EDM_ML_DEBUG
0961   edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot1 << ":" << ntot2;
0962 #endif
0963   loadGeometryHexagon8(layers, trforms, firstLayer, php);
0964 }
0965 
0966 void HGCalGeomParameters::loadGeometryHexagonModule(const cms::DDCompactView* cpv,
0967                                                     HGCalParameters& php,
0968                                                     const std::string& sdTag1,
0969                                                     const std::string& sdTag2,
0970                                                     int firstLayer) {
0971 #ifdef EDM_ML_DEBUG
0972   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters (DD4hep)::loadGeometryHexagonModule called with tags " << sdTag1
0973                                 << ":" << sdTag2 << " firstLayer " << firstLayer;
0974   int ntot1(0), ntot2(0);
0975 #endif
0976   std::map<int, HGCalGeomParameters::layerParameters> layers;
0977   std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0978   std::map<std::pair<int, int>, double> zvals;
0979   int levelTop = php.levelT_[0];
0980 
0981   const cms::DDFilter filter1("Volume", sdTag2);
0982   cms::DDFilteredView fv1((*cpv), filter1);
0983   while (fv1.firstChild()) {
0984 #ifdef EDM_ML_DEBUG
0985     ++ntot1;
0986 #endif
0987     int nsiz = static_cast<int>(fv1.level());
0988     if (nsiz > levelTop) {
0989       std::vector<int> copy = fv1.copyNos();
0990       int lay = copy[nsiz - levelTop - 1];
0991       int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
0992       if (zside != 1)
0993         zside = -1;
0994       if (lay == 0) {
0995         throw cms::Exception("DDException")
0996             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0997       } else {
0998         if (zvals.find(std::make_pair(lay, zside)) == zvals.end()) {
0999           zvals[std::make_pair(lay, zside)] = HGCalParameters::k_ScaleFromDD4hep * fv1.translation().Z();
1000 #ifdef EDM_ML_DEBUG
1001           std::ostringstream st1;
1002           st1 << "Name0 " << fv1.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
1003               << nsiz;
1004           for (const auto& c : copy)
1005             st1 << ":" << c;
1006           st1 << " Z " << zvals[std::make_pair(lay, zside)];
1007           edm::LogVerbatim("HGCalGeom") << st1.str();
1008 #endif
1009         }
1010       }
1011     }
1012   }
1013 
1014   const cms::DDFilter filter2("Volume", sdTag1);
1015   cms::DDFilteredView fv2((*cpv), filter2);
1016   while (fv2.firstChild()) {
1017     // Layers first
1018     int nsiz = static_cast<int>(fv2.level());
1019 #ifdef EDM_ML_DEBUG
1020     ++ntot2;
1021 #endif
1022     if (nsiz > levelTop) {
1023       std::vector<int> copy = fv2.copyNos();
1024       int lay = copy[nsiz - levelTop - 1];
1025       int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
1026       if (zside != 1)
1027         zside = -1;
1028 #ifdef EDM_ML_DEBUG
1029       std::ostringstream st2;
1030       st2 << "Name1 " << fv2.name() << "Shape " << cms::dd::name(cms::DDSolidShapeMap, fv2.shape()) << " LTop "
1031           << levelTop << ":" << lay << " ZSide " << zside << ":" << php.levelZSide_ << " # of levels " << nsiz;
1032       for (const auto& c : copy)
1033         st2 << ":" << c;
1034       edm::LogVerbatim("HGCalGeom") << st2.str();
1035 #endif
1036       if (lay == 0) {
1037         throw cms::Exception("DDException")
1038             << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
1039       } else {
1040         if (zvals.find(std::make_pair(lay, zside)) != zvals.end()) {
1041           if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
1042             php.layer_.emplace_back(lay);
1043           auto itr = layers.find(lay);
1044           if (itr == layers.end()) {
1045             const std::vector<double>& pars = fv2.parameters();
1046             double rin = HGCalParameters::k_ScaleFromDD4hep * pars[0];
1047             double rout = (php.firstMixedLayer_ > 0 && lay >= php.firstMixedLayer_)
1048                               ? php.radiusMixBoundary_[lay - php.firstMixedLayer_]
1049                               : HGCalParameters::k_ScaleFromDD4hep * pars[1];
1050             double zp = zvals[std::make_pair(lay, 1)];
1051             HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
1052             layers[lay] = laypar;
1053 #ifdef EDM_ML_DEBUG
1054             std::ostringstream st3;
1055             st3 << "Name2 " << fv2.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
1056                 << nsiz;
1057             for (const auto& c : copy)
1058               st3 << ":" << c;
1059             st3 << " R " << rin << ":" << rout;
1060             edm::LogVerbatim("HGCalGeom") << st3.str();
1061 #endif
1062           }
1063 
1064           if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
1065             DD3Vector x, y, z;
1066             fv2.rotation().GetComponents(x, y, z);
1067             const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
1068             const CLHEP::HepRotation hr(rotation);
1069             double xx = ((std::abs(fv2.translation().X()) < tolerance)
1070                              ? 0
1071                              : HGCalParameters::k_ScaleFromDD4hep * fv2.translation().X());
1072             double yy = ((std::abs(fv2.translation().Y()) < tolerance)
1073                              ? 0
1074                              : HGCalParameters::k_ScaleFromDD4hep * fv2.translation().Y());
1075             const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
1076             HGCalParameters::hgtrform mytrf;
1077             mytrf.zp = zside;
1078             mytrf.lay = lay;
1079             mytrf.sec = 0;
1080             mytrf.subsec = 0;
1081             mytrf.h3v = h3v;
1082             mytrf.hr = hr;
1083             trforms[std::make_pair(lay, zside)] = mytrf;
1084 #ifdef EDM_ML_DEBUG
1085             edm::LogVerbatim("HGCalGeom") << "Translation " << h3v;
1086 #endif
1087           }
1088         }
1089       }
1090     }
1091   }
1092 #ifdef EDM_ML_DEBUG
1093   edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot1 << ":" << ntot2;
1094 #endif
1095   loadGeometryHexagon8(layers, trforms, firstLayer, php);
1096 }
1097 
1098 void HGCalGeomParameters::loadGeometryHexagon8(const std::map<int, HGCalGeomParameters::layerParameters>& layers,
1099                                                std::map<std::pair<int, int>, HGCalParameters::hgtrform>& trforms,
1100                                                const int& firstLayer,
1101                                                HGCalParameters& php) {
1102   double rmin(0), rmax(0);
1103   for (unsigned int i = 0; i < layers.size(); ++i) {
1104     for (auto& layer : layers) {
1105       if (layer.first == static_cast<int>(i + firstLayer)) {
1106         php.layerIndex_.emplace_back(i);
1107         php.rMinLayHex_.emplace_back(layer.second.rmin);
1108         php.rMaxLayHex_.emplace_back(layer.second.rmax);
1109         php.zLayerHex_.emplace_back(layer.second.zpos);
1110         if (i == 0) {
1111           rmin = layer.second.rmin;
1112           rmax = layer.second.rmax;
1113         } else {
1114           if (rmin > layer.second.rmin)
1115             rmin = layer.second.rmin;
1116           if (rmax < layer.second.rmax)
1117             rmax = layer.second.rmax;
1118         }
1119         break;
1120       }
1121     }
1122   }
1123   php.rLimit_.emplace_back(rmin);
1124   php.rLimit_.emplace_back(rmax);
1125   php.depth_ = php.layer_;
1126   php.depthIndex_ = php.layerIndex_;
1127   php.depthLayerF_ = php.layerIndex_;
1128 
1129   for (unsigned int i = 0; i < php.layer_.size(); ++i) {
1130     for (auto& trform : trforms) {
1131       if (trform.first.first == static_cast<int>(i + firstLayer)) {
1132         php.fillTrForm(trform.second);
1133       }
1134     }
1135   }
1136 
1137 #ifdef EDM_ML_DEBUG
1138   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Minimum/maximum R " << php.rLimit_[0] << ":" << php.rLimit_[1];
1139   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.zLayerHex_.size() << " layers";
1140   for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
1141     int k = php.layerIndex_[i];
1142     edm::LogVerbatim("HGCalGeom") << "Layer[" << i << ":" << k << ":" << php.layer_[k]
1143                                   << "] with r = " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
1144                                   << " at z = " << php.zLayerHex_[i];
1145   }
1146   edm::LogVerbatim("HGCalGeom") << "Obtained " << php.trformIndex_.size() << " transformation matrices";
1147   for (unsigned int k = 0; k < php.trformIndex_.size(); ++k) {
1148     edm::LogVerbatim("HGCalGeom") << "Matrix[" << k << "] (" << std::hex << php.trformIndex_[k] << std::dec
1149                                   << ") Translation (" << php.trformTranX_[k] << ", " << php.trformTranY_[k] << ", "
1150                                   << php.trformTranZ_[k] << " Rotation (" << php.trformRotXX_[k] << ", "
1151                                   << php.trformRotYX_[k] << ", " << php.trformRotZX_[k] << ", " << php.trformRotXY_[k]
1152                                   << ", " << php.trformRotYY_[k] << ", " << php.trformRotZY_[k] << ", "
1153                                   << php.trformRotXZ_[k] << ", " << php.trformRotYZ_[k] << ", " << php.trformRotZZ_[k]
1154                                   << ")";
1155   }
1156 #endif
1157 }
1158 
1159 void HGCalGeomParameters::loadSpecParsHexagon(const DDFilteredView& fv,
1160                                               HGCalParameters& php,
1161                                               const DDCompactView* cpv,
1162                                               const std::string& sdTag1,
1163                                               const std::string& sdTag2) {
1164   DDsvalues_type sv(fv.mergedSpecifics());
1165   php.boundR_ = getDDDArray("RadiusBound", sv, 4);
1166   rescale(php.boundR_, HGCalParameters::k_ScaleFromDDD);
1167   php.rLimit_ = getDDDArray("RadiusLimits", sv, 2);
1168   rescale(php.rLimit_, HGCalParameters::k_ScaleFromDDD);
1169   php.levelT_ = dbl_to_int(getDDDArray("LevelTop", sv, 0));
1170 
1171   // Grouping of layers
1172   php.layerGroup_ = dbl_to_int(getDDDArray("GroupingZFine", sv, 0));
1173   php.layerGroupM_ = dbl_to_int(getDDDArray("GroupingZMid", sv, 0));
1174   php.layerGroupO_ = dbl_to_int(getDDDArray("GroupingZOut", sv, 0));
1175   php.slopeMin_ = getDDDArray("Slope", sv, 1);
1176   const auto& dummy2 = getDDDArray("LayerOffset", sv, 0);
1177   if (!dummy2.empty())
1178     php.layerOffset_ = dummy2[0];
1179   else
1180     php.layerOffset_ = 0;
1181 
1182   // Wafer size
1183   std::string attribute = "Volume";
1184   DDSpecificsMatchesValueFilter filter1{DDValue(attribute, sdTag1, 0.0)};
1185   DDFilteredView fv1(*cpv, filter1);
1186   if (fv1.firstChild()) {
1187     DDsvalues_type sv(fv1.mergedSpecifics());
1188     const auto& dummy = getDDDArray("WaferSize", sv, 0);
1189     waferSize_ = dummy[0];
1190   }
1191 
1192   // Cell size
1193   DDSpecificsMatchesValueFilter filter2{DDValue(attribute, sdTag2, 0.0)};
1194   DDFilteredView fv2(*cpv, filter2);
1195   if (fv2.firstChild()) {
1196     DDsvalues_type sv(fv2.mergedSpecifics());
1197     php.cellSize_ = getDDDArray("CellSize", sv, 0);
1198   }
1199 
1200   loadSpecParsHexagon(php);
1201 }
1202 
1203 void HGCalGeomParameters::loadSpecParsHexagon(const cms::DDFilteredView& fv,
1204                                               HGCalParameters& php,
1205                                               const std::string& sdTag1,
1206                                               const std::string& sdTag2,
1207                                               const std::string& sdTag3,
1208                                               const std::string& sdTag4) {
1209   php.boundR_ = fv.get<std::vector<double> >(sdTag4, "RadiusBound");
1210   rescale(php.boundR_, HGCalParameters::k_ScaleFromDD4hep);
1211   php.rLimit_ = fv.get<std::vector<double> >(sdTag4, "RadiusLimits");
1212   rescale(php.rLimit_, HGCalParameters::k_ScaleFromDD4hep);
1213   php.levelT_ = dbl_to_int(fv.get<std::vector<double> >(sdTag4, "LevelTop"));
1214 
1215   // Grouping of layers
1216   php.layerGroup_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZFine"));
1217   php.layerGroupM_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZMid"));
1218   php.layerGroupO_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZOut"));
1219   php.slopeMin_ = fv.get<std::vector<double> >(sdTag4, "Slope");
1220   if (php.slopeMin_.empty())
1221     php.slopeMin_.emplace_back(0);
1222 
1223   // Wafer size
1224   const auto& dummy = fv.get<std::vector<double> >(sdTag2, "WaferSize");
1225   waferSize_ = dummy[0] * HGCalParameters::k_ScaleFromDD4hepToG4;
1226 
1227   // Cell size
1228   php.cellSize_ = fv.get<std::vector<double> >(sdTag3, "CellSize");
1229   rescale(php.cellSize_, HGCalParameters::k_ScaleFromDD4hepToG4);
1230 
1231   // Layer Offset
1232   const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1233   if (!dummy2.empty()) {
1234     php.layerOffset_ = dummy2[0];
1235   } else {
1236     php.layerOffset_ = 0;
1237   }
1238 
1239   loadSpecParsHexagon(php);
1240 }
1241 
1242 void HGCalGeomParameters::loadSpecParsHexagon(const HGCalParameters& php) {
1243 #ifdef EDM_ML_DEBUG
1244   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer radius ranges"
1245                                 << " for cell grouping " << php.boundR_[0] << ":" << php.boundR_[1] << ":"
1246                                 << php.boundR_[2] << ":" << php.boundR_[3];
1247   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Minimum/maximum R " << php.rLimit_[0] << ":" << php.rLimit_[1];
1248   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: LevelTop " << php.levelT_[0];
1249   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: minimum slope " << php.slopeMin_[0] << " and layer groupings "
1250                                 << "for the 3 ranges:";
1251   for (unsigned int k = 0; k < php.layerGroup_.size(); ++k)
1252     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerGroup_[k] << ":" << php.layerGroupM_[k] << ":"
1253                                   << php.layerGroupO_[k];
1254   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Wafer Size: " << waferSize_;
1255   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: " << php.cellSize_.size() << " cells of sizes:";
1256   for (unsigned int k = 0; k < php.cellSize_.size(); ++k)
1257     edm::LogVerbatim("HGCalGeom") << " [" << k << "] " << php.cellSize_[k];
1258   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: First Layer " << php.firstLayer_ << " and layer offset "
1259                                 << php.layerOffset_;
1260 #endif
1261 }
1262 
1263 void HGCalGeomParameters::loadSpecParsHexagon8(const DDFilteredView& fv, HGCalParameters& php) {
1264   DDsvalues_type sv(fv.mergedSpecifics());
1265   php.cellThickness_ = getDDDArray("CellThickness", sv, 3);
1266   rescale(php.cellThickness_, HGCalParameters::k_ScaleFromDDD);
1267   if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1268       (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1269     php.waferThickness_ = getDDDArray("WaferThickness", sv, 3);
1270     rescale(php.waferThickness_, HGCalParameters::k_ScaleFromDDD);
1271   } else {
1272     for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1273       php.waferThickness_.emplace_back(php.waferThick_);
1274   }
1275 
1276   php.radius100to200_ = getDDDArray("Radius100to200", sv, 5);
1277   php.radius200to300_ = getDDDArray("Radius200to300", sv, 5);
1278 
1279   const auto& dummy = getDDDArray("RadiusCuts", sv, 4);
1280   php.choiceType_ = static_cast<int>(dummy[0]);
1281   php.nCornerCut_ = static_cast<int>(dummy[1]);
1282   php.fracAreaMin_ = dummy[2];
1283   php.zMinForRad_ = HGCalParameters::k_ScaleFromDDD * dummy[3];
1284 
1285   php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1286   rescale(php.radiusMixBoundary_, HGCalParameters::k_ScaleFromDDD);
1287 
1288   php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1289   php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1290   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDDD);
1291   php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1292   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDDD);
1293 
1294   php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1295   php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1296   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDDD);
1297   php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1298   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDDD);
1299 
1300   php.zRanges_ = fv.vector("ZRanges");
1301   rescale(php.zRanges_, HGCalParameters::k_ScaleFromDDD);
1302 
1303   const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1304   php.layerOffset_ = dummy2[0];
1305   php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1306 
1307   if (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell) {
1308     const auto& dummy3 = fv.vector("CalibCellRadius");
1309     php.calibCellRHD_ = HGCalParameters::k_ScaleFromDDD * dummy3[0];
1310     php.calibCellFullHD_ = dbl_to_int(fv.vector("CalibCellFullHD"));
1311     php.calibCellPartHD_ = dbl_to_int(fv.vector("CalibCellPartHD"));
1312     php.calibCellRLD_ = HGCalParameters::k_ScaleFromDDD * dummy3[1];
1313     php.calibCellFullLD_ = dbl_to_int(fv.vector("CalibCellFullLD"));
1314     php.calibCellPartLD_ = dbl_to_int(fv.vector("CalibCellPartLD"));
1315   }
1316 
1317   loadSpecParsHexagon8(php);
1318 
1319   // Read in parameters from Philip's file
1320   if (php.waferMaskMode_ > 1) {
1321     std::vector<int> layerType, waferIndex, waferProperties;
1322     std::vector<double> cassetteShift;
1323     if (php.waferMaskMode_ == siliconFileEE) {
1324       waferIndex = dbl_to_int(fv.vector("WaferIndexEE"));
1325       waferProperties = dbl_to_int(fv.vector("WaferPropertiesEE"));
1326     } else if (php.waferMaskMode_ == siliconCassetteEE) {
1327       waferIndex = dbl_to_int(fv.vector("WaferIndexEE"));
1328       waferProperties = dbl_to_int(fv.vector("WaferPropertiesEE"));
1329       cassetteShift = fv.vector("CassetteShiftEE");
1330     } else if (php.waferMaskMode_ == siliconFileHE) {
1331       waferIndex = dbl_to_int(fv.vector("WaferIndexHE"));
1332       waferProperties = dbl_to_int(fv.vector("WaferPropertiesHE"));
1333     } else if (php.waferMaskMode_ == siliconCassetteHE) {
1334       waferIndex = dbl_to_int(fv.vector("WaferIndexHE"));
1335       waferProperties = dbl_to_int(fv.vector("WaferPropertiesHE"));
1336       cassetteShift = fv.vector("CassetteShiftHE");
1337     }
1338     if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1339         (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1340       if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1341         layerType = dbl_to_int(fv.vector("LayerTypesEE"));
1342       } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1343         layerType = dbl_to_int(fv.vector("LayerTypesHE"));
1344       }
1345     }
1346 
1347     php.cassetteShift_ = cassetteShift;
1348     rescale(php.cassetteShift_, HGCalParameters::k_ScaleFromDDD);
1349     loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1350   }
1351 }
1352 
1353 void HGCalGeomParameters::loadSpecParsHexagon8(const cms::DDFilteredView& fv,
1354                                                const cms::DDVectorsMap& vmap,
1355                                                HGCalParameters& php,
1356                                                const std::string& sdTag1) {
1357   php.cellThickness_ = fv.get<std::vector<double> >(sdTag1, "CellThickness");
1358   rescale(php.cellThickness_, HGCalParameters::k_ScaleFromDD4hep);
1359   if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1360       (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1361     php.waferThickness_ = fv.get<std::vector<double> >(sdTag1, "WaferThickness");
1362     rescale(php.waferThickness_, HGCalParameters::k_ScaleFromDD4hep);
1363   } else {
1364     for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1365       php.waferThickness_.emplace_back(php.waferThick_);
1366   }
1367 
1368   php.radius100to200_ = fv.get<std::vector<double> >(sdTag1, "Radius100to200");
1369   php.radius200to300_ = fv.get<std::vector<double> >(sdTag1, "Radius200to300");
1370 
1371   const auto& dummy = fv.get<std::vector<double> >(sdTag1, "RadiusCuts");
1372   if (dummy.size() > 3) {
1373     php.choiceType_ = static_cast<int>(dummy[0]);
1374     php.nCornerCut_ = static_cast<int>(dummy[1]);
1375     php.fracAreaMin_ = dummy[2];
1376     php.zMinForRad_ = HGCalParameters::k_ScaleFromDD4hep * dummy[3];
1377   } else {
1378     php.choiceType_ = php.nCornerCut_ = php.fracAreaMin_ = php.zMinForRad_ = 0;
1379   }
1380 
1381   php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1382   php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1383   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDD4hep);
1384   php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1385   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDD4hep);
1386 
1387   php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1388   php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1389   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDD4hep);
1390   php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1391   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDD4hep);
1392   unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1393   for (unsigned int k = 0; k < kmax; ++k)
1394     php.slopeTop_.emplace_back(0);
1395 
1396   const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1397   if (!dummy2.empty()) {
1398     php.layerOffset_ = dummy2[0];
1399   } else {
1400     php.layerOffset_ = 0;
1401   }
1402 
1403   for (auto const& it : vmap) {
1404     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1405       for (const auto& i : it.second)
1406         php.radiusMixBoundary_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1407     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1408       for (const auto& i : it.second)
1409         php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1410     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1411       for (const auto& i : it.second)
1412         php.layerCenter_.emplace_back(std::round(i));
1413     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellRadius")) {
1414       php.calibCellRHD_ = HGCalParameters::k_ScaleFromDD4hep * it.second[0];
1415       php.calibCellRLD_ = HGCalParameters::k_ScaleFromDDD * it.second[1];
1416     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellFullHD")) {
1417       for (const auto& i : it.second)
1418         php.calibCellFullHD_.emplace_back(std::round(i));
1419     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellPartHD")) {
1420       for (const auto& i : it.second)
1421         php.calibCellPartHD_.emplace_back(std::round(i));
1422     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellFullLD")) {
1423       for (const auto& i : it.second)
1424         php.calibCellFullLD_.emplace_back(std::round(i));
1425     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellPartLD")) {
1426       for (const auto& i : it.second)
1427         php.calibCellPartLD_.emplace_back(std::round(i));
1428     }
1429   }
1430 
1431   loadSpecParsHexagon8(php);
1432 
1433   // Read in parameters from Philip's file
1434   if (php.waferMaskMode_ > 1) {
1435     std::vector<int> layerType, waferIndex, waferProperties;
1436     std::vector<double> cassetteShift;
1437     if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1438       for (auto const& it : vmap) {
1439         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferIndexEE")) {
1440           for (const auto& i : it.second)
1441             waferIndex.emplace_back(std::round(i));
1442         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferPropertiesEE")) {
1443           for (const auto& i : it.second)
1444             waferProperties.emplace_back(std::round(i));
1445         }
1446       }
1447       if (php.waferMaskMode_ == siliconCassetteEE) {
1448         for (auto const& it : vmap) {
1449           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftEE")) {
1450             for (const auto& i : it.second)
1451               cassetteShift.emplace_back(i);
1452           }
1453         }
1454       }
1455     } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1456       for (auto const& it : vmap) {
1457         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferIndexHE")) {
1458           for (const auto& i : it.second)
1459             waferIndex.emplace_back(std::round(i));
1460         } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferPropertiesHE")) {
1461           for (const auto& i : it.second)
1462             waferProperties.emplace_back(std::round(i));
1463         }
1464       }
1465       if (php.waferMaskMode_ == siliconCassetteHE) {
1466         for (auto const& it : vmap) {
1467           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1468             for (const auto& i : it.second)
1469               cassetteShift.emplace_back(i);
1470           }
1471         }
1472       }
1473     }
1474     if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1475         (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1476       if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconFileHE)) {
1477         for (auto const& it : vmap) {
1478           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerTypesEE")) {
1479             for (const auto& i : it.second)
1480               layerType.emplace_back(std::round(i));
1481           }
1482         }
1483       } else if (php.waferMaskMode_ == siliconFileHE) {
1484         for (auto const& it : vmap) {
1485           if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerTypesHE")) {
1486             for (const auto& i : it.second)
1487               layerType.emplace_back(std::round(i));
1488           }
1489         }
1490       }
1491     }
1492 
1493     php.cassetteShift_ = cassetteShift;
1494     rescale(php.cassetteShift_, HGCalParameters::k_ScaleFromDD4hep);
1495     loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1496   }
1497 }
1498 
1499 void HGCalGeomParameters::loadSpecParsHexagon8(HGCalParameters& php) {
1500 #ifdef EDM_ML_DEBUG
1501   for (unsigned int k = 0; k < php.waferThickness_.size(); ++k)
1502     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer[" << k << "] Thickness " << php.waferThickness_[k];
1503   for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1504     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: cell[" << k << "] Thickness " << php.cellThickness_[k];
1505   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Polynomial "
1506                                 << "parameters for 120 to 200 micron "
1507                                 << "transition with" << php.radius100to200_.size() << " elements";
1508   for (unsigned int k = 0; k < php.radius100to200_.size(); ++k)
1509     edm::LogVerbatim("HGCalGeom") << "Element [" << k << "] " << php.radius100to200_[k];
1510   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Polynomial "
1511                                 << "parameters for 200 to 300 micron "
1512                                 << "transition with " << php.radius200to300_.size() << " elements";
1513   for (unsigned int k = 0; k < php.radius200to300_.size(); ++k)
1514     edm::LogVerbatim("HGCalGeom") << "Element [" << k << "] " << php.radius200to300_[k];
1515   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Parameters for the"
1516                                 << " transition " << php.choiceType_ << ":" << php.nCornerCut_ << ":"
1517                                 << php.fracAreaMin_ << ":" << php.zMinForRad_;
1518   for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k)
1519     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k];
1520   for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1521     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1522                                   << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1523   for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1524     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1525                                   << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1526   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary " << php.zRanges_[0] << ":" << php.zRanges_[1] << ":"
1527                                 << php.zRanges_[2] << ":" << php.zRanges_[3];
1528   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1529                                 << php.layerCenter_.size();
1530   for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1531     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1532   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Calibration cells in HD having radius " << php.calibCellRHD_
1533                                 << " in arrays of size " << php.calibCellFullHD_.size() << ":"
1534                                 << php.calibCellPartHD_.size();
1535   for (unsigned int k = 0; k < php.calibCellFullHD_.size(); ++k)
1536     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.calibCellFullHD_[k] << ":" << php.calibCellPartHD_[k];
1537   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Calibration cells in LD having radius " << php.calibCellRLD_
1538                                 << " in arrays of size " << php.calibCellFullLD_.size() << ":"
1539                                 << php.calibCellPartLD_.size();
1540   for (unsigned int k = 0; k < php.calibCellFullLD_.size(); ++k)
1541     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.calibCellFullLD_[k] << ":" << php.calibCellPartLD_[k];
1542 #endif
1543 }
1544 
1545 void HGCalGeomParameters::loadSpecParsHexagon8(HGCalParameters& php,
1546                                                const std::vector<int>& layerType,
1547                                                const std::vector<int>& waferIndex,
1548                                                const std::vector<int>& waferProperties) {
1549   // Store parameters from Philip's file
1550   for (unsigned int k = 0; k < layerType.size(); ++k) {
1551     php.layerType_.emplace_back(HGCalTypes::layerType(layerType[k]));
1552 #ifdef EDM_ML_DEBUG
1553     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Type " << layerType[k] << ":" << php.layerType_[k];
1554 #endif
1555   }
1556   for (unsigned int k = 0; k < php.layerType_.size(); ++k) {
1557     double cth = (php.layerType_[k] == HGCalTypes::WaferCenterR) ? cos(php.layerRotation_) : 1.0;
1558     double sth = (php.layerType_[k] == HGCalTypes::WaferCenterR) ? sin(php.layerRotation_) : 0.0;
1559     php.layerRotV_.emplace_back(std::make_pair(cth, sth));
1560 #ifdef EDM_ML_DEBUG
1561     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Type " << php.layerType_[k] << " cos|sin(Theta) "
1562                                   << php.layerRotV_[k].first << ":" << php.layerRotV_[k].second;
1563 #endif
1564   }
1565   for (unsigned int k = 0; k < waferIndex.size(); ++k) {
1566     int partial = HGCalProperty::waferPartial(waferProperties[k]);
1567     int orient =
1568         ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) || (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
1569             ? HGCalProperty::waferOrient(waferProperties[k])
1570             : HGCalWaferMask::getRotation(php.waferZSide_, partial, HGCalProperty::waferOrient(waferProperties[k]));
1571     php.waferInfoMap_[waferIndex[k]] = HGCalParameters::waferInfo(HGCalProperty::waferThick(waferProperties[k]),
1572                                                                   partial,
1573                                                                   orient,
1574                                                                   HGCalProperty::waferCassette(waferProperties[k]));
1575 #ifdef EDM_ML_DEBUG
1576     edm::LogVerbatim("HGCalGeom") << "[" << k << ":" << waferIndex[k] << ":"
1577                                   << HGCalWaferIndex::waferLayer(waferIndex[k]) << ":"
1578                                   << HGCalWaferIndex::waferU(waferIndex[k]) << ":"
1579                                   << HGCalWaferIndex::waferV(waferIndex[k]) << "]  Thickness type "
1580                                   << HGCalProperty::waferThick(waferProperties[k]) << " Partial type " << partial
1581                                   << " Orientation " << HGCalProperty::waferOrient(waferProperties[k]) << ":" << orient;
1582 #endif
1583   }
1584 }
1585 
1586 void HGCalGeomParameters::loadSpecParsTrapezoid(const DDFilteredView& fv, HGCalParameters& php) {
1587   DDsvalues_type sv(fv.mergedSpecifics());
1588   php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1589   rescale(php.radiusMixBoundary_, HGCalParameters::k_ScaleFromDDD);
1590 
1591   php.nPhiBinBH_ = dbl_to_int(getDDDArray("NPhiBinBH", sv, 0));
1592   php.layerFrontBH_ = dbl_to_int(getDDDArray("LayerFrontBH", sv, 0));
1593   php.rMinLayerBH_ = getDDDArray("RMinLayerBH", sv, 0);
1594   rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDDD);
1595   php.nCellsFine_ = php.nPhiBinBH_[0];
1596   php.nCellsCoarse_ = php.nPhiBinBH_[1];
1597   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1598   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1599 
1600   php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1601   php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1602   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDDD);
1603   php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1604   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDDD);
1605 
1606   php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1607   php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1608   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDDD);
1609   php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1610   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDDD);
1611 
1612   php.zRanges_ = fv.vector("ZRanges");
1613   rescale(php.zRanges_, HGCalParameters::k_ScaleFromDDD);
1614 
1615   // Offsets
1616   const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1617   php.layerOffset_ = dummy2[0];
1618   php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1619 
1620   loadSpecParsTrapezoid(php);
1621 
1622   // tile parameters from Katja's file
1623   if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette)) {
1624     std::vector<int> tileIndx, tileProperty;
1625     std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4;
1626     std::vector<double> tileRMin, tileRMax;
1627     std::vector<int> tileRingMin, tileRingMax;
1628     std::vector<double> cassetteShift;
1629     tileIndx = dbl_to_int(fv.vector("TileIndex"));
1630     tileProperty = dbl_to_int(fv.vector("TileProperty"));
1631     tileHEX1 = dbl_to_int(fv.vector("TileHEX1"));
1632     tileHEX2 = dbl_to_int(fv.vector("TileHEX2"));
1633     tileHEX3 = dbl_to_int(fv.vector("TileHEX3"));
1634     tileHEX4 = dbl_to_int(fv.vector("TileHEX4"));
1635     tileRMin = fv.vector("TileRMin");
1636     tileRMax = fv.vector("TileRMax");
1637     rescale(tileRMin, HGCalParameters::k_ScaleFromDDD);
1638     rescale(tileRMax, HGCalParameters::k_ScaleFromDDD);
1639     tileRingMin = dbl_to_int(fv.vector("TileRingMin"));
1640     tileRingMax = dbl_to_int(fv.vector("TileRingMax"));
1641     if (php.waferMaskMode_ == scintillatorCassette) {
1642       if (php.cassettes_ > 0)
1643         php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1644       cassetteShift = fv.vector("CassetteShiftHE");
1645       rescale(cassetteShift, HGCalParameters::k_ScaleFromDDD);
1646     }
1647 
1648     php.cassetteShift_ = cassetteShift;
1649     loadSpecParsTrapezoid(php,
1650                           tileIndx,
1651                           tileProperty,
1652                           tileHEX1,
1653                           tileHEX2,
1654                           tileHEX3,
1655                           tileHEX4,
1656                           tileRMin,
1657                           tileRMax,
1658                           tileRingMin,
1659                           tileRingMax);
1660   }
1661 }
1662 
1663 void HGCalGeomParameters::loadSpecParsTrapezoid(const cms::DDFilteredView& fv,
1664                                                 const cms::DDVectorsMap& vmap,
1665                                                 HGCalParameters& php,
1666                                                 const std::string& sdTag1) {
1667   for (auto const& it : vmap) {
1668     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1669       for (const auto& i : it.second)
1670         php.radiusMixBoundary_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1671     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1672       for (const auto& i : it.second)
1673         php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1674     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1675       for (const auto& i : it.second)
1676         php.layerCenter_.emplace_back(std::round(i));
1677     }
1678   }
1679 
1680   php.nPhiBinBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "NPhiBinBH"));
1681   php.layerFrontBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "LayerFrontBH"));
1682   php.rMinLayerBH_ = fv.get<std::vector<double> >(sdTag1, "RMinLayerBH");
1683   rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDD4hep);
1684   php.nCellsFine_ = php.nPhiBinBH_[0];
1685   php.nCellsCoarse_ = php.nPhiBinBH_[1];
1686   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1687   php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1688 
1689   php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1690   php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1691   rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDD4hep);
1692   php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1693   rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDD4hep);
1694 
1695   php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1696   php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1697   rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDD4hep);
1698   php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1699   rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDD4hep);
1700   unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1701   for (unsigned int k = 0; k < kmax; ++k)
1702     php.slopeTop_.emplace_back(0);
1703 
1704   const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1705   php.layerOffset_ = dummy2[0];
1706 
1707   loadSpecParsTrapezoid(php);
1708 
1709   // tile parameters from Katja's file
1710   if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette)) {
1711     std::vector<int> tileIndx, tileProperty;
1712     std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4;
1713     std::vector<double> tileRMin, tileRMax;
1714     std::vector<int> tileRingMin, tileRingMax;
1715     std::vector<double> cassetteShift;
1716     for (auto const& it : vmap) {
1717       if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileIndex")) {
1718         for (const auto& i : it.second)
1719           tileIndx.emplace_back(std::round(i));
1720       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileProperty")) {
1721         for (const auto& i : it.second)
1722           tileProperty.emplace_back(std::round(i));
1723       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX1")) {
1724         for (const auto& i : it.second)
1725           tileHEX1.emplace_back(std::round(i));
1726       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX2")) {
1727         for (const auto& i : it.second)
1728           tileHEX2.emplace_back(std::round(i));
1729       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX3")) {
1730         for (const auto& i : it.second)
1731           tileHEX3.emplace_back(std::round(i));
1732       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX4")) {
1733         for (const auto& i : it.second)
1734           tileHEX4.emplace_back(std::round(i));
1735       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin")) {
1736         for (const auto& i : it.second)
1737           tileRMin.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1738       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax")) {
1739         for (const auto& i : it.second)
1740           tileRMax.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1741       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin")) {
1742         for (const auto& i : it.second)
1743           tileRingMin.emplace_back(std::round(i));
1744       } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax")) {
1745         for (const auto& i : it.second)
1746           tileRingMax.emplace_back(std::round(i));
1747       }
1748     }
1749     if (php.waferMaskMode_ == scintillatorCassette) {
1750       for (auto const& it : vmap) {
1751         if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1752           for (const auto& i : it.second)
1753             cassetteShift.emplace_back(i);
1754         }
1755       }
1756     }
1757 
1758     rescale(cassetteShift, HGCalParameters::k_ScaleFromDD4hep);
1759     php.cassetteShift_ = cassetteShift;
1760     loadSpecParsTrapezoid(php,
1761                           tileIndx,
1762                           tileProperty,
1763                           tileHEX1,
1764                           tileHEX2,
1765                           tileHEX3,
1766                           tileHEX4,
1767                           tileRMin,
1768                           tileRMax,
1769                           tileRingMin,
1770                           tileRingMax);
1771   }
1772 }
1773 
1774 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php) {
1775 #ifdef EDM_ML_DEBUG
1776   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters:nCells " << php.nCellsFine_ << ":" << php.nCellsCoarse_
1777                                 << " cellSize: " << php.cellSize_[0] << ":" << php.cellSize_[1];
1778   for (unsigned int k = 0; k < php.layerFrontBH_.size(); ++k)
1779     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Type[" << k << "] Front Layer = " << php.layerFrontBH_[k]
1780                                   << " rMin = " << php.rMinLayerBH_[k];
1781   for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k) {
1782     edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k]
1783                                   << " Nphi = " << php.scintCells(k + php.firstLayer_)
1784                                   << " dPhi = " << php.scintCellSize(k + php.firstLayer_);
1785   }
1786 
1787   for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1788     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1789                                   << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1790 
1791   for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1792     edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1793                                   << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1794 
1795   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary " << php.zRanges_[0] << ":" << php.zRanges_[1] << ":"
1796                                 << php.zRanges_[2] << ":" << php.zRanges_[3];
1797 
1798   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1799                                 << php.layerCenter_.size();
1800   for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1801     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1802 #endif
1803 }
1804 
1805 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php,
1806                                                 const std::vector<int>& tileIndx,
1807                                                 const std::vector<int>& tileProperty,
1808                                                 const std::vector<int>& tileHEX1,
1809                                                 const std::vector<int>& tileHEX2,
1810                                                 const std::vector<int>& tileHEX3,
1811                                                 const std::vector<int>& tileHEX4,
1812                                                 const std::vector<double>& tileRMin,
1813                                                 const std::vector<double>& tileRMax,
1814                                                 const std::vector<int>& tileRingMin,
1815                                                 const std::vector<int>& tileRingMax) {
1816   // tile parameters from Katja's file
1817   for (unsigned int k = 0; k < tileIndx.size(); ++k) {
1818     php.tileInfoMap_[tileIndx[k]] = HGCalParameters::tileInfo(HGCalProperty::tileType(tileProperty[k]),
1819                                                               HGCalProperty::tileSiPM(tileProperty[k]),
1820                                                               tileHEX1[k],
1821                                                               tileHEX2[k],
1822                                                               tileHEX3[k],
1823                                                               tileHEX4[k]);
1824 #ifdef EDM_ML_DEBUG
1825     edm::LogVerbatim("HGCalGeom") << "Tile[" << k << ":" << tileIndx[k] << "] "
1826                                   << " Type " << HGCalProperty::tileType(tileProperty[k]) << " SiPM "
1827                                   << HGCalProperty::tileSiPM(tileProperty[k]) << " HEX " << std::hex << tileHEX1[k]
1828                                   << ":" << tileHEX2[k] << ":" << tileHEX3[k] << ":" << tileHEX4[k] << std::dec;
1829 #endif
1830   }
1831 
1832   for (unsigned int k = 0; k < tileRMin.size(); ++k) {
1833     php.tileRingR_.emplace_back(tileRMin[k], tileRMax[k]);
1834 #ifdef EDM_ML_DEBUG
1835     edm::LogVerbatim("HGCalGeom") << "TileRingR[" << k << "] " << tileRMin[k] << ":" << tileRMax[k];
1836 #endif
1837   }
1838 
1839   for (unsigned k = 0; k < tileRingMin.size(); ++k) {
1840     php.tileRingRange_.emplace_back(tileRingMin[k], tileRingMax[k]);
1841 #ifdef EDM_ML_DEBUG
1842     edm::LogVerbatim("HGCalGeom") << "TileRingRange[" << k << "] " << tileRingMin[k] << ":" << tileRingMax[k];
1843 #endif
1844   }
1845 }
1846 
1847 void HGCalGeomParameters::loadWaferHexagon(HGCalParameters& php) {
1848   double waferW(HGCalParameters::k_ScaleFromDDD * waferSize_), rmin(HGCalParameters::k_ScaleFromDDD * php.waferR_);
1849   double rin(php.rLimit_[0]), rout(php.rLimit_[1]), rMaxFine(php.boundR_[1]);
1850 #ifdef EDM_ML_DEBUG
1851   edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << rmin << " R Limits: " << rin << ":" << rout
1852                                 << " Fine " << rMaxFine;
1853 #endif
1854   // Clear the vectors
1855   php.waferCopy_.clear();
1856   php.waferTypeL_.clear();
1857   php.waferTypeT_.clear();
1858   php.waferPosX_.clear();
1859   php.waferPosY_.clear();
1860   double dx = 0.5 * waferW;
1861   double dy = 3.0 * dx * tan(30._deg);
1862   double rr = 2.0 * dx * tan(30._deg);
1863   int ncol = static_cast<int>(2.0 * rout / waferW) + 1;
1864   int nrow = static_cast<int>(rout / (waferW * tan(30._deg))) + 1;
1865   int ns2 = (2 * ncol + 1) * (2 * nrow + 1) * php.layer_.size();
1866   int incm(0), inrm(0);
1867   HGCalParameters::layer_map copiesInLayers(php.layer_.size() + 1);
1868   HGCalParameters::waferT_map waferTypes(ns2 + 1);
1869 #ifdef EDM_ML_DEBUG
1870   int kount(0), ntot(0);
1871   edm::LogVerbatim("HGCalGeom") << "Row " << nrow << " Column " << ncol;
1872 #endif
1873   for (int nr = -nrow; nr <= nrow; ++nr) {
1874     int inr = (nr >= 0) ? nr : -nr;
1875     for (int nc = -ncol; nc <= ncol; ++nc) {
1876       int inc = (nc >= 0) ? nc : -nc;
1877       if (inr % 2 == inc % 2) {
1878         double xpos = nc * dx;
1879         double ypos = nr * dy;
1880         std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, rin, rout, true);
1881         double rpos = std::sqrt(xpos * xpos + ypos * ypos);
1882         int typet = (rpos < rMaxFine) ? 1 : 2;
1883         int typel(3);
1884         for (int k = 1; k < 4; ++k) {
1885           if ((rpos + rmin) <= php.boundR_[k]) {
1886             typel = k;
1887             break;
1888           }
1889         }
1890 #ifdef EDM_ML_DEBUG
1891         ++ntot;
1892 #endif
1893         if (corner.first > 0) {
1894           int copy = HGCalTypes::packTypeUV(typel, nc, nr);
1895           if (inc > incm)
1896             incm = inc;
1897           if (inr > inrm)
1898             inrm = inr;
1899 #ifdef EDM_ML_DEBUG
1900           kount++;
1901           edm::LogVerbatim("HGCalGeom") << kount << ":" << ntot << " Copy " << copy << " Type " << typel << ":" << typet
1902                                         << " Location " << corner.first << " Position " << xpos << ":" << ypos
1903                                         << " Layers " << php.layer_.size();
1904 #endif
1905           php.waferCopy_.emplace_back(copy);
1906           php.waferTypeL_.emplace_back(typel);
1907           php.waferTypeT_.emplace_back(typet);
1908           php.waferPosX_.emplace_back(xpos);
1909           php.waferPosY_.emplace_back(ypos);
1910           for (unsigned int il = 0; il < php.layer_.size(); ++il) {
1911             std::pair<int, int> corner =
1912                 HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, php.rMinLayHex_[il], php.rMaxLayHex_[il], true);
1913             if (corner.first > 0) {
1914               auto cpy = copiesInLayers[php.layer_[il]].find(copy);
1915               if (cpy == copiesInLayers[php.layer_[il]].end())
1916                 copiesInLayers[php.layer_[il]][copy] =
1917                     ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ? php.waferCopy_.size() : -1);
1918             }
1919             if ((corner.first > 0) && (corner.first < static_cast<int>(HGCalParameters::k_CornerSize))) {
1920               int wl = HGCalWaferIndex::waferIndex(php.layer_[il], copy, 0, true);
1921               waferTypes[wl] = corner;
1922             }
1923           }
1924         }
1925       }
1926     }
1927   }
1928   php.copiesInLayers_ = copiesInLayers;
1929   php.waferTypes_ = waferTypes;
1930   php.nSectors_ = static_cast<int>(php.waferCopy_.size());
1931   php.waferUVMax_ = 0;
1932 #ifdef EDM_ML_DEBUG
1933   edm::LogVerbatim("HGCalGeom") << "HGCalWaferHexagon: # of columns " << incm << " # of rows " << inrm << " and "
1934                                 << kount << ":" << ntot << " wafers; R " << rin << ":" << rout;
1935   edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
1936   for (unsigned int k = 0; k < copiesInLayers.size(); ++k) {
1937     const auto& theModules = copiesInLayers[k];
1938     edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
1939     int k2(0);
1940     for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
1941       edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
1942     }
1943   }
1944 #endif
1945 }
1946 
1947 void HGCalGeomParameters::loadWaferHexagon8(HGCalParameters& php) {
1948   double waferW(php.waferSize_);
1949   double waferS(php.sensorSeparation_);
1950   auto wType = std::make_unique<HGCalWaferType>(php.radius100to200_,
1951                                                 php.radius200to300_,
1952                                                 HGCalParameters::k_ScaleToDDD * (waferW + waferS),
1953                                                 HGCalParameters::k_ScaleToDDD * php.zMinForRad_,
1954                                                 php.choiceType_,
1955                                                 php.nCornerCut_,
1956                                                 php.fracAreaMin_);
1957 
1958   double rout(php.rLimit_[1]);
1959 #ifdef EDM_ML_DEBUG
1960   edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << waferS << " R Max: " << rout;
1961 #endif
1962   // Clear the vectors
1963   php.waferCopy_.clear();
1964   php.waferTypeL_.clear();
1965   php.waferTypeT_.clear();
1966   php.waferPosX_.clear();
1967   php.waferPosY_.clear();
1968   double r = 0.5 * (waferW + waferS);
1969   double R = 2.0 * r / sqrt3_;
1970   double dy = 0.75 * R;
1971   double r1 = 0.5 * waferW;
1972   double R1 = 2.0 * r1 / sqrt3_;
1973   int N = (r == 0) ? 3 : (static_cast<int>(0.5 * rout / r) + 3);
1974   int ns1 = (2 * N + 1) * (2 * N + 1);
1975   int ns2 = ns1 * php.zLayerHex_.size();
1976 #ifdef EDM_ML_DEBUG
1977   edm::LogVerbatim("HGCalGeom") << "wafer " << waferW << ":" << waferS << " r " << r << " dy " << dy << " N " << N
1978                                 << " sizes " << ns1 << ":" << ns2;
1979   std::vector<int> indtypes(ns1 + 1);
1980   indtypes.clear();
1981 #endif
1982   HGCalParameters::wafer_map wafersInLayers(ns1 + 1);
1983   HGCalParameters::wafer_map typesInLayers(ns2 + 1);
1984   HGCalParameters::waferT_map waferTypes(ns2 + 1);
1985   int ipos(0), lpos(0), uvmax(0), nwarn(0);
1986   std::vector<int> uvmx(php.zLayerHex_.size(), 0);
1987   for (int v = -N; v <= N; ++v) {
1988     for (int u = -N; u <= N; ++u) {
1989       int nr = 2 * v;
1990       int nc = -2 * u + v;
1991       double xpos = nc * r;
1992       double ypos = nr * dy;
1993       int indx = HGCalWaferIndex::waferIndex(0, u, v);
1994       php.waferCopy_.emplace_back(indx);
1995       php.waferPosX_.emplace_back(xpos);
1996       php.waferPosY_.emplace_back(ypos);
1997       wafersInLayers[indx] = ipos;
1998       ++ipos;
1999       std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, r1, R1, 0, rout, false);
2000       if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2001           ((corner.first > 0) && php.defineFull_)) {
2002         uvmax = std::max(uvmax, std::max(std::abs(u), std::abs(v)));
2003       }
2004       for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
2005         int copy = i + php.layerOffset_;
2006         std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], (waferW + waferS));
2007         int lay = php.layer_[php.layerIndex_[i]];
2008         double xpos0 = xpos + xyoff.first;
2009         double ypos0 = ypos + xyoff.second;
2010         double zpos = php.zLayerHex_[i];
2011         int kndx = HGCalWaferIndex::waferIndex(lay, u, v);
2012         int type(-1);
2013         if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2014             (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) || (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
2015           type = wType->getType(kndx, php.waferInfoMap_);
2016         if (type < 0)
2017           type = wType->getType(HGCalParameters::k_ScaleToDDD * xpos0,
2018                                 HGCalParameters::k_ScaleToDDD * ypos0,
2019                                 HGCalParameters::k_ScaleToDDD * zpos);
2020         php.waferTypeL_.emplace_back(type);
2021         typesInLayers[kndx] = lpos;
2022         ++lpos;
2023 #ifdef EDM_ML_DEBUG
2024         indtypes.emplace_back(kndx);
2025 #endif
2026         std::pair<int, int> corner =
2027             HGCalGeomTools::waferCorner(xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], false);
2028 #ifdef EDM_ML_DEBUG
2029         if (((corner.first == 0) && std::abs(u) < 5 && std::abs(v) < 5) || (std::abs(u) < 2 && std::abs(v) < 2)) {
2030           edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " R " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2031                                         << " u " << u << " v " << v << " with " << corner.first << " corners";
2032         }
2033 #endif
2034         if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2035             ((corner.first > 0) && php.defineFull_)) {
2036           uvmx[i] = std::max(uvmx[i], std::max(std::abs(u), std::abs(v)));
2037         }
2038         if ((corner.first < static_cast<int>(HGCalParameters::k_CornerSize)) && (corner.first > 0)) {
2039 #ifdef EDM_ML_DEBUG
2040           edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2041                                         << corner.first << ":" << corner.second;
2042 #endif
2043           int wl = HGCalWaferIndex::waferIndex(lay, u, v);
2044           if (php.waferMaskMode_ > 0) {
2045             std::pair<int, int> corner0 = HGCalWaferMask::getTypeMode(
2046                 xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], type, php.waferMaskMode_);
2047             if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2048                 (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2049                 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
2050               auto itr = php.waferInfoMap_.find(wl);
2051               if (itr != php.waferInfoMap_.end()) {
2052                 int part = (itr->second).part;
2053                 int orient = (itr->second).orient;
2054                 bool ok = ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2055                            (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
2056                               ? true
2057                               : HGCalWaferMask::goodTypeMode(
2058                                     xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], part, orient, false);
2059                 if (ok)
2060                   corner0 = std::make_pair(part, (HGCalTypes::k_OffsetRotation + orient));
2061 #ifdef EDM_ML_DEBUG
2062                 edm::LogVerbatim("HGCalGeom")
2063                     << "Layer:u:v " << i << ":" << lay << ":" << u << ":" << v << " Part " << corner0.first << ":"
2064                     << part << " Orient " << corner0.second << ":" << orient << " Position " << xpos0 << ":" << ypos0
2065                     << " delta " << r1 << ":" << R1 << " Limit " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2066                     << " Compatibiliety Flag " << ok;
2067 #endif
2068                 if (!ok)
2069                   ++nwarn;
2070               }
2071             }
2072             waferTypes[wl] = corner0;
2073 #ifdef EDM_ML_DEBUG
2074             edm::LogVerbatim("HGCalGeom")
2075                 << "Layer " << lay << " u|v " << u << ":" << v << " Index " << std::hex << wl << std::dec << " pos "
2076                 << xpos0 << ":" << ypos0 << " R " << r1 << ":" << R1 << " Range " << php.rMinLayHex_[i] << ":"
2077                 << php.rMaxLayHex_[i] << type << ":" << php.waferMaskMode_ << " corner " << corner.first << ":"
2078                 << corner.second << " croner0 " << corner0.first << ":" << corner0.second;
2079 #endif
2080           } else {
2081             waferTypes[wl] = corner;
2082 #ifdef EDM_ML_DEBUG
2083             edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2084                                           << corner.first << ":" << corner.second;
2085 #endif
2086           }
2087         }
2088       }
2089     }
2090   }
2091   if (nwarn > 0)
2092     edm::LogWarning("HGCalGeom") << "HGCalGeomParameters::loadWafer8: there are " << nwarn
2093                                  << " wafers with non-matching partial- orientation types";
2094   php.waferUVMax_ = uvmax;
2095   php.waferUVMaxLayer_ = uvmx;
2096   php.wafersInLayers_ = wafersInLayers;
2097   php.typesInLayers_ = typesInLayers;
2098   php.waferTypes_ = waferTypes;
2099   php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2100   HGCalParameters::hgtrap mytr;
2101   mytr.lay = 1;
2102   mytr.bl = php.waferR_;
2103   mytr.tl = php.waferR_;
2104   mytr.h = php.waferR_;
2105   mytr.alpha = 0.0;
2106   mytr.cellSize = HGCalParameters::k_ScaleToDDD * php.waferSize_;
2107   for (auto const& dz : php.cellThickness_) {
2108     mytr.dz = 0.5 * HGCalParameters::k_ScaleToDDD * dz;
2109     php.fillModule(mytr, false);
2110   }
2111   for (unsigned k = 0; k < php.cellThickness_.size(); ++k) {
2112     HGCalParameters::hgtrap mytr = php.getModule(k, false);
2113     mytr.bl *= HGCalParameters::k_ScaleFromDDD;
2114     mytr.tl *= HGCalParameters::k_ScaleFromDDD;
2115     mytr.h *= HGCalParameters::k_ScaleFromDDD;
2116     mytr.dz *= HGCalParameters::k_ScaleFromDDD;
2117     mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2118     php.fillModule(mytr, true);
2119   }
2120 #ifdef EDM_ML_DEBUG
2121   edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Total of " << php.waferCopy_.size() << " wafers";
2122   for (unsigned int k = 0; k < php.waferCopy_.size(); ++k) {
2123     int id = php.waferCopy_[k];
2124     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << std::hex << id << std::dec << ":"
2125                                   << HGCalWaferIndex::waferLayer(id) << ":" << HGCalWaferIndex::waferU(id) << ":"
2126                                   << HGCalWaferIndex::waferV(id) << " x " << php.waferPosX_[k] << " y "
2127                                   << php.waferPosY_[k] << " index " << php.wafersInLayers_[id];
2128   }
2129   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Total of " << php.waferTypeL_.size() << " wafer types";
2130   for (unsigned int k = 0; k < php.waferTypeL_.size(); ++k) {
2131     int id = indtypes[k];
2132     edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.typesInLayers_[id] << ":" << php.waferTypeL_[k] << " ID "
2133                                   << std::hex << id << std::dec << ":" << HGCalWaferIndex::waferLayer(id) << ":"
2134                                   << HGCalWaferIndex::waferU(id) << ":" << HGCalWaferIndex::waferV(id);
2135   }
2136 #endif
2137 
2138   //Wafer offset
2139   php.xLayerHex_.clear();
2140   php.yLayerHex_.clear();
2141   double waferSize = php.waferSize_ + php.sensorSeparation_;
2142 #ifdef EDM_ML_DEBUG
2143   edm::LogVerbatim("HGCalGeom") << "WaferSize " << waferSize;
2144 #endif
2145   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2146     int copy = k + php.layerOffset_;
2147     std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], waferSize);
2148     php.xLayerHex_.emplace_back(xyoff.first);
2149     php.yLayerHex_.emplace_back(xyoff.second);
2150 #ifdef EDM_ML_DEBUG
2151     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Off " << copy << ":" << php.layerCenter_[copy] << " Shift "
2152                                   << xyoff.first << ":" << xyoff.second;
2153 #endif
2154   }
2155 }
2156 
2157 void HGCalGeomParameters::loadCellParsHexagon(const DDCompactView* cpv, HGCalParameters& php) {
2158   // Special parameters for cell parameters
2159   std::string attribute = "OnlyForHGCalNumbering";
2160   DDSpecificsHasNamedValueFilter filter1{attribute};
2161   DDFilteredView fv1(*cpv, filter1);
2162   bool ok = fv1.firstChild();
2163 
2164   if (ok) {
2165     php.cellFine_ = dbl_to_int(cpv->vector("waferFine"));
2166     php.cellCoarse_ = dbl_to_int(cpv->vector("waferCoarse"));
2167   }
2168 
2169   loadCellParsHexagon(php);
2170 }
2171 
2172 void HGCalGeomParameters::loadCellParsHexagon(const cms::DDVectorsMap& vmap, HGCalParameters& php) {
2173   for (auto const& it : vmap) {
2174     if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferFine")) {
2175       for (const auto& i : it.second)
2176         php.cellFine_.emplace_back(std::round(i));
2177     } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferCoarse")) {
2178       for (const auto& i : it.second)
2179         php.cellCoarse_.emplace_back(std::round(i));
2180     }
2181   }
2182 
2183   loadCellParsHexagon(php);
2184 }
2185 
2186 void HGCalGeomParameters::loadCellParsHexagon(const HGCalParameters& php) {
2187 #ifdef EDM_ML_DEBUG
2188   edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellFine_.size() << " rows for fine cells";
2189   for (unsigned int k = 0; k < php.cellFine_.size(); ++k)
2190     edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellFine_[k];
2191   edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellCoarse_.size() << " rows for coarse cells";
2192   for (unsigned int k = 0; k < php.cellCoarse_.size(); ++k)
2193     edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellCoarse_[k];
2194 #endif
2195 }
2196 
2197 void HGCalGeomParameters::loadCellTrapezoid(HGCalParameters& php) {
2198   php.xLayerHex_.resize(php.zLayerHex_.size(), 0);
2199   php.yLayerHex_.resize(php.zLayerHex_.size(), 0);
2200 #ifdef EDM_ML_DEBUG
2201   edm::LogVerbatim("HGCalGeom") << "HGCalParameters: x|y|zLayerHex in array of size " << php.zLayerHex_.size();
2202   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k)
2203     edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Shift " << php.xLayerHex_[k] << ":" << php.yLayerHex_[k] << ":"
2204                                   << php.zLayerHex_[k];
2205 #endif
2206   // Find the radius of each eta-partitions
2207 
2208   if ((php.mode_ == HGCalGeometryMode::TrapezoidFile) || (php.mode_ == HGCalGeometryMode::TrapezoidModule) ||
2209       (php.mode_ == HGCalGeometryMode::TrapezoidCassette)) {
2210     //Ring radii for each partition
2211     for (unsigned int k = 0; k < 2; ++k) {
2212       for (unsigned int kk = 0; kk < php.tileRingR_.size(); ++kk) {
2213         php.radiusLayer_[k].emplace_back(php.tileRingR_[kk].first);
2214 #ifdef EDM_ML_DEBUG
2215         double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2216                               : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2217         double rv = php.radiusLayer_[k].back();
2218         double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2219         edm::LogVerbatim("HGCalGeom") << "New [" << kk << "] new R = " << rv << " Eta = " << eta;
2220 #endif
2221       }
2222       php.radiusLayer_[k].emplace_back(php.tileRingR_[php.tileRingR_.size() - 1].second);
2223     }
2224     // Minimum and maximum radius index for each layer
2225     for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2226       php.iradMinBH_.emplace_back(1 + php.tileRingRange_[k].first);
2227       php.iradMaxBH_.emplace_back(1 + php.tileRingRange_[k].second);
2228 #ifdef EDM_ML_DEBUG
2229       int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2230       edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge " << php.iradMinBH_.back()
2231                                     << " Top edge " << php.iradMaxBH_.back();
2232 #endif
2233     }
2234   } else {
2235     //Ring radii for each partition
2236     for (unsigned int k = 0; k < 2; ++k) {
2237       double rmax = ((k == 0) ? (php.rMaxLayHex_[php.layerFrontBH_[1] - php.firstLayer_] - 1)
2238                               : (php.rMaxLayHex_[php.rMaxLayHex_.size() - 1]));
2239       double rv = php.rMinLayerBH_[k];
2240       double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2241                             : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2242       php.radiusLayer_[k].emplace_back(rv);
2243 #ifdef EDM_ML_DEBUG
2244       double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2245       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] rmax " << rmax << " Z = " << zv
2246                                     << " dEta = " << php.cellSize_[k] << "\nOld[0] new R = " << rv << " Eta = " << eta;
2247       int kount(1);
2248 #endif
2249       while (rv < rmax) {
2250         double eta = -(php.cellSize_[k] + std::log(std::tan(0.5 * std::atan(rv / zv))));
2251         rv = zv * std::tan(2.0 * std::atan(std::exp(-eta)));
2252         php.radiusLayer_[k].emplace_back(rv);
2253 #ifdef EDM_ML_DEBUG
2254         edm::LogVerbatim("HGCalGeom") << "Old [" << kount << "] new R = " << rv << " Eta = " << eta;
2255         ++kount;
2256 #endif
2257       }
2258     }
2259     // Find minimum and maximum radius index for each layer
2260     for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2261       int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2262       std::vector<double>::iterator low, high;
2263       low = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMinLayHex_[k]);
2264 #ifdef EDM_ML_DEBUG
2265       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RLow = " << php.rMinLayHex_[k] << " pos "
2266                                     << static_cast<int>(low - php.radiusLayer_[kk].begin());
2267 #endif
2268       if (low == php.radiusLayer_[kk].begin())
2269         ++low;
2270       int irlow = static_cast<int>(low - php.radiusLayer_[kk].begin());
2271       double drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2272 #ifdef EDM_ML_DEBUG
2273       edm::LogVerbatim("HGCalGeom") << "irlow " << irlow << " dr " << drlow << " min " << php.minTileSize_;
2274 #endif
2275       if (drlow < php.minTileSize_) {
2276         ++irlow;
2277 #ifdef EDM_ML_DEBUG
2278         drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2279         edm::LogVerbatim("HGCalGeom") << "Modified irlow " << irlow << " dr " << drlow;
2280 #endif
2281       }
2282       high = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMaxLayHex_[k]);
2283 #ifdef EDM_ML_DEBUG
2284       edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RHigh = " << php.rMaxLayHex_[k] << " pos "
2285                                     << static_cast<int>(high - php.radiusLayer_[kk].begin());
2286 #endif
2287       if (high == php.radiusLayer_[kk].end())
2288         --high;
2289       int irhigh = static_cast<int>(high - php.radiusLayer_[kk].begin());
2290       double drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2291 #ifdef EDM_ML_DEBUG
2292       edm::LogVerbatim("HGCalGeom") << "irhigh " << irhigh << " dr " << drhigh << " min " << php.minTileSize_;
2293 #endif
2294       if (drhigh < php.minTileSize_) {
2295         --irhigh;
2296 #ifdef EDM_ML_DEBUG
2297         drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2298         edm::LogVerbatim("HGCalGeom") << "Modified irhigh " << irhigh << " dr " << drhigh;
2299 #endif
2300       }
2301       php.iradMinBH_.emplace_back(irlow);
2302       php.iradMaxBH_.emplace_back(irhigh);
2303 #ifdef EDM_ML_DEBUG
2304       edm::LogVerbatim("HGCalGeom") << "Old Layer " << k << " Type " << kk << " Low edge " << irlow << ":" << drlow
2305                                     << " Top edge " << irhigh << ":" << drhigh;
2306 #endif
2307     }
2308   }
2309 #ifdef EDM_ML_DEBUG
2310   for (unsigned int k = 0; k < 2; ++k) {
2311     edm::LogVerbatim("HGCalGeom") << "Type " << k << " with " << php.radiusLayer_[k].size() << " radii";
2312     for (unsigned int kk = 0; kk < php.radiusLayer_[k].size(); ++kk)
2313       edm::LogVerbatim("HGCalGeom") << "Ring[" << kk << "] " << php.radiusLayer_[k][kk];
2314   }
2315 #endif
2316 
2317   // Now define the volumes
2318   int im(0);
2319   php.waferUVMax_ = 0;
2320   HGCalParameters::hgtrap mytr;
2321   mytr.alpha = 0.0;
2322   for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2323     if (php.iradMaxBH_[k] > php.waferUVMax_)
2324       php.waferUVMax_ = php.iradMaxBH_[k];
2325     int kk = ((php.firstLayer_ + static_cast<int>(k)) < php.layerFrontBH_[1]) ? 0 : 1;
2326     int irm = php.radiusLayer_[kk].size() - 1;
2327 #ifdef EDM_ML_DEBUG
2328     double rmin = php.radiusLayer_[kk][std::max((php.iradMinBH_[k] - 1), 0)];
2329     double rmax = php.radiusLayer_[kk][std::min(php.iradMaxBH_[k], irm)];
2330     edm::LogVerbatim("HGCalGeom") << "Layer " << php.firstLayer_ + k << ":" << kk << " Radius range "
2331                                   << php.iradMinBH_[k] << ":" << php.iradMaxBH_[k] << ":" << rmin << ":" << rmax;
2332 #endif
2333     mytr.lay = php.firstLayer_ + k;
2334     for (int irad = php.iradMinBH_[k]; irad <= php.iradMaxBH_[k]; ++irad) {
2335       double rmin = php.radiusLayer_[kk][std::max((irad - 1), 0)];
2336       double rmax = php.radiusLayer_[kk][std::min(irad, irm)];
2337       mytr.bl = 0.5 * rmin * php.scintCellSize(mytr.lay);
2338       mytr.tl = 0.5 * rmax * php.scintCellSize(mytr.lay);
2339       mytr.h = 0.5 * (rmax - rmin);
2340       mytr.dz = 0.5 * php.waferThick_;
2341       mytr.cellSize = 0.5 * (rmax + rmin) * php.scintCellSize(mytr.lay);
2342       php.fillModule(mytr, true);
2343       mytr.bl *= HGCalParameters::k_ScaleToDDD;
2344       mytr.tl *= HGCalParameters::k_ScaleToDDD;
2345       mytr.h *= HGCalParameters::k_ScaleToDDD;
2346       mytr.dz *= HGCalParameters::k_ScaleToDDD;
2347       mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2348       php.fillModule(mytr, false);
2349       if (irad == php.iradMinBH_[k])
2350         php.firstModule_.emplace_back(im);
2351       ++im;
2352       if (irad == php.iradMaxBH_[k] - 1)
2353         php.lastModule_.emplace_back(im);
2354     }
2355   }
2356   php.nSectors_ = php.waferUVMax_;
2357 #ifdef EDM_ML_DEBUG
2358   edm::LogVerbatim("HGCalGeom") << "Maximum radius index " << php.waferUVMax_;
2359   for (unsigned int k = 0; k < php.firstModule_.size(); ++k)
2360     edm::LogVerbatim("HGCalGeom") << "Layer " << k + php.firstLayer_ << " Modules " << php.firstModule_[k] << ":"
2361                                   << php.lastModule_[k];
2362 #endif
2363 }
2364 
2365 std::vector<double> HGCalGeomParameters::getDDDArray(const std::string& str, const DDsvalues_type& sv, const int nmin) {
2366   DDValue value(str);
2367   if (DDfetch(&sv, value)) {
2368     const std::vector<double>& fvec = value.doubles();
2369     int nval = fvec.size();
2370     if (nmin > 0) {
2371       if (nval < nmin) {
2372         throw cms::Exception("DDException")
2373             << "HGCalGeomParameters:  # of " << str << " bins " << nval << " < " << nmin << " ==> illegal";
2374       }
2375     } else {
2376       if (nval < 1 && nmin == 0) {
2377         throw cms::Exception("DDException")
2378             << "HGCalGeomParameters: # of " << str << " bins " << nval << " < 1 ==> illegal"
2379             << " (nmin=" << nmin << ")";
2380       }
2381     }
2382     return fvec;
2383   } else {
2384     if (nmin >= 0) {
2385       throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str;
2386     }
2387     std::vector<double> fvec;
2388     return fvec;
2389   }
2390 }
2391 
2392 std::pair<double, double> HGCalGeomParameters::cellPosition(
2393     const std::vector<HGCalGeomParameters::cellParameters>& wafers,
2394     std::vector<HGCalGeomParameters::cellParameters>::const_iterator& itrf,
2395     int wafer,
2396     double xx,
2397     double yy) {
2398   if (itrf == wafers.end()) {
2399     for (std::vector<HGCalGeomParameters::cellParameters>::const_iterator itr = wafers.begin(); itr != wafers.end();
2400          ++itr) {
2401       if (itr->wafer == wafer) {
2402         itrf = itr;
2403         break;
2404       }
2405     }
2406   }
2407   double dx(0), dy(0);
2408   if (itrf != wafers.end()) {
2409     dx = (xx - itrf->xyz.x());
2410     if (std::abs(dx) < tolerance)
2411       dx = 0;
2412     dy = (yy - itrf->xyz.y());
2413     if (std::abs(dy) < tolerance)
2414       dy = 0;
2415   }
2416   return std::make_pair(dx, dy);
2417 }
2418 
2419 void HGCalGeomParameters::rescale(std::vector<double>& v, const double s) {
2420   std::for_each(v.begin(), v.end(), [s](double& n) { n *= s; });
2421 }
2422 
2423 void HGCalGeomParameters::resetZero(std::vector<double>& v) {
2424   for (auto& n : v) {
2425     if (std::abs(n) < tolmin)
2426       n = 0;
2427   }
2428 }
[ 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!