Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-02 02:40:40 

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