Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-06-20 00:21:44 

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