Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:15:10

0001 /* for High Granularity Calorimeter
0002  * This geometry is essentially driven by topology, 
0003  * which is thus encapsulated in this class. 
0004  * This makes this geometry not suitable to be loaded
0005  * by regular CaloGeometryLoader<T>
0006  */
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "DataFormats/GeometrySurface/interface/Plane.h"
0010 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0011 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0012 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0014 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0015 
0016 #include <cmath>
0017 
0018 #include <Math/Transform3D.h>
0019 #include <Math/EulerAngles.h>
0020 
0021 typedef CaloCellGeometry::Tr3D Tr3D;
0022 typedef std::vector<float> ParmVec;
0023 
0024 //#define EDM_ML_DEBUG
0025 
0026 const bool debugLocate = false;
0027 
0028 HGCalGeometry::HGCalGeometry(const HGCalTopology& topology_)
0029     : m_topology(topology_),
0030       m_validGeomIds(topology_.totalGeomModules()),
0031       m_det(topology_.detector()),
0032       m_subdet(topology_.subDetector()),
0033       twoBysqrt3_(2.0 / std::sqrt(3.0)) {
0034   if (m_det == DetId::HGCalHSc) {
0035     m_cellVec2 = CellVec2(topology_.totalGeomModules());
0036   } else {
0037     m_cellVec = CellVec(topology_.totalGeomModules());
0038   }
0039   m_validIds.reserve(m_topology.totalModules());
0040 #ifdef EDM_ML_DEBUG
0041   edm::LogVerbatim("HGCalGeom") << "Expected total # of Geometry Modules " << m_topology.totalGeomModules();
0042 #endif
0043 }
0044 
0045 void HGCalGeometry::fillNamedParams(DDFilteredView fv) {}
0046 
0047 void HGCalGeometry::initializeParms() {}
0048 
0049 void HGCalGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int i, Pt3D& ref) {
0050   if (m_det == DetId::HGCalHSc) {
0051     FlatTrd::localCorners(lc, pv, ref);
0052   } else {
0053     FlatHexagon::localCorners(lc, pv, ref);
0054   }
0055 }
0056 
0057 void HGCalGeometry::newCell(
0058     const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0059   DetId geomId = getGeometryDetId(detId);
0060   int cells(0);
0061   HGCalTopology::DecodedDetId id = m_topology.decode(detId);
0062   if (m_topology.waferHexagon6()) {
0063     cells = m_topology.dddConstants().numberCellsHexagon(id.iSec1);
0064 #ifdef EDM_ML_DEBUG
0065     edm::LogVerbatim("HGCalGeom") << "NewCell " << HGCalDetId(detId) << " GEOM " << HGCalDetId(geomId);
0066 #endif
0067   } else if (m_topology.tileTrapezoid()) {
0068     cells = 1;
0069 #ifdef EDM_ML_DEBUG
0070     edm::LogVerbatim("HGCalGeom") << "NewCell " << HGCScintillatorDetId(detId) << " GEOM "
0071                                   << HGCScintillatorDetId(geomId);
0072 #endif
0073   } else {
0074     cells = m_topology.dddConstants().numberCellsHexagon(id.iLay, id.iSec1, id.iSec2, false);
0075 #ifdef EDM_ML_DEBUG
0076     if (detId.det() == DetId::Forward)
0077       edm::LogVerbatim("HGCalGeom") << "NewCell " << HFNoseDetId(detId) << " GEOM " << HFNoseDetId(geomId);
0078     else
0079       edm::LogVerbatim("HGCalGeom") << "NewCell " << HGCSiliconDetId(detId) << " GEOM " << HGCSiliconDetId(geomId);
0080 #endif
0081   }
0082   const uint32_t cellIndex(m_topology.detId2denseGeomId(geomId));
0083 
0084   if (m_det == DetId::HGCalHSc) {
0085     m_cellVec2.at(cellIndex) = FlatTrd(cornersMgr(), f1, f2, f3, parm);
0086   } else {
0087     m_cellVec.at(cellIndex) = FlatHexagon(cornersMgr(), f1, f2, f3, parm);
0088   }
0089   m_validGeomIds.at(cellIndex) = geomId;
0090 
0091 #ifdef EDM_ML_DEBUG
0092   edm::LogVerbatim("HGCalGeom") << "Store for DetId " << std::hex << detId.rawId() << " GeomId " << geomId.rawId()
0093                                 << std::dec << " Index " << cellIndex << " cells " << cells;
0094   unsigned int nOld = m_validIds.size();
0095 #endif
0096   if (m_topology.waferHexagon6()) {
0097     for (int cell = 0; cell < cells; ++cell) {
0098       id.iCell1 = cell;
0099       DetId idc = m_topology.encode(id);
0100       if (m_topology.valid(idc)) {
0101         m_validIds.emplace_back(idc);
0102 #ifdef EDM_ML_DEBUG
0103         edm::LogVerbatim("HGCalGeom") << "Valid Id [" << cell << "] " << HGCalDetId(idc);
0104 #endif
0105       }
0106     }
0107   } else if (m_topology.tileTrapezoid()) {
0108     DetId idc = m_topology.encode(id);
0109     if (m_topology.valid(idc)) {
0110       HGCScintillatorDetId hid(idc);
0111       std::pair<int, int> typm = m_topology.dddConstants().tileType(hid.layer(), hid.ring(), 0);
0112       if (typm.first >= 0) {
0113         hid.setType(typm.first);
0114         hid.setSiPM(typm.second);
0115         idc = static_cast<DetId>(hid);
0116       }
0117       m_validIds.emplace_back(idc);
0118 #ifdef EDM_ML_DEBUG
0119       edm::LogVerbatim("HGCalGeom") << "Valid Id [0] " << HGCScintillatorDetId(idc);
0120 #endif
0121     } else {
0122       edm::LogWarning("HGCalGeom") << "Check " << HGCScintillatorDetId(idc) << " from " << HGCScintillatorDetId(detId)
0123                                    << " ERROR ???";
0124     }
0125   } else {
0126 #ifdef EDM_ML_DEBUG
0127     unsigned int cellAll(0), cellSelect(0);
0128 #endif
0129     for (int u = 0; u < 2 * cells; ++u) {
0130       for (int v = 0; v < 2 * cells; ++v) {
0131         if (((v - u) < cells) && (u - v) <= cells) {
0132           id.iCell1 = u;
0133           id.iCell2 = v;
0134           DetId idc = m_topology.encode(id);
0135 #ifdef EDM_ML_DEBUG
0136           ++cellAll;
0137 #endif
0138           if (m_topology.dddConstants().cellInLayer(id.iSec1, id.iSec2, u, v, id.iLay, id.zSide, true)) {
0139             m_validIds.emplace_back(idc);
0140 #ifdef EDM_ML_DEBUG
0141             ++cellSelect;
0142             if (idc.det() == DetId::Forward)
0143               edm::LogVerbatim("HGCalGeom") << "Valid Id [" << u << ", " << v << "] " << HFNoseDetId(idc);
0144             else
0145               edm::LogVerbatim("HGCalGeom") << "Valid Id [" << u << ", " << v << "] " << HGCSiliconDetId(idc);
0146 #endif
0147           }
0148         }
0149       }
0150     }
0151 #ifdef EDM_ML_DEBUG
0152     edm::LogVerbatim("HGCalGeom") << "HGCalGeometry keeps " << cellSelect << " out of " << cellAll << " for wafer "
0153                                   << id.iSec1 << ":" << id.iSec2 << " in "
0154                                   << " layer " << id.iLay;
0155 #endif
0156   }
0157 #ifdef EDM_ML_DEBUG
0158   if (m_det == DetId::HGCalHSc) {
0159     edm::LogVerbatim("HGCalGeom") << "HGCalGeometry::newCell-> [" << cellIndex << "]"
0160                                   << " front:" << f1.x() << '/' << f1.y() << '/' << f1.z() << " back:" << f2.x() << '/'
0161                                   << f2.y() << '/' << f2.z() << " eta|phi " << m_cellVec2[cellIndex].etaPos() << ":"
0162                                   << m_cellVec2[cellIndex].phiPos();
0163   } else {
0164     edm::LogVerbatim("HGCalGeom") << "HGCalGeometry::newCell-> [" << cellIndex << "]"
0165                                   << " front:" << f1.x() << '/' << f1.y() << '/' << f1.z() << " back:" << f2.x() << '/'
0166                                   << f2.y() << '/' << f2.z() << " eta|phi " << m_cellVec[cellIndex].etaPos() << ":"
0167                                   << m_cellVec[cellIndex].phiPos();
0168   }
0169   unsigned int nNew = m_validIds.size();
0170   if (m_topology.waferHexagon6()) {
0171     edm::LogVerbatim("HGCalGeom") << "ID: " << HGCalDetId(detId) << " with valid DetId from " << nOld << " to " << nNew;
0172   } else if (m_topology.tileTrapezoid()) {
0173     edm::LogVerbatim("HGCalGeom") << "ID: " << HGCScintillatorDetId(detId) << " with valid DetId from " << nOld
0174                                   << " to " << nNew;
0175   } else if (m_topology.isHFNose()) {
0176     edm::LogVerbatim("HGCalGeom") << "ID: " << HFNoseDetId(detId) << " with valid DetId from " << nOld << " to "
0177                                   << nNew;
0178   } else {
0179     edm::LogVerbatim("HGCalGeom") << "ID: " << HGCSiliconDetId(detId) << " with valid DetId from " << nOld << " to "
0180                                   << nNew;
0181   }
0182   edm::LogVerbatim("HGCalGeom") << "Cell[" << cellIndex << "] " << std::hex << geomId.rawId() << ":"
0183                                 << m_validGeomIds[cellIndex].rawId() << std::dec;
0184 #endif
0185 }
0186 
0187 std::shared_ptr<const CaloCellGeometry> HGCalGeometry::getGeometry(const DetId& detId) const {
0188   if (detId == DetId())
0189     return nullptr;  // nothing to get
0190   DetId geomId = getGeometryDetId(detId);
0191   const uint32_t cellIndex(m_topology.detId2denseGeomId(geomId));
0192   const GlobalPoint pos = (detId != geomId) ? getPosition(detId, false) : GlobalPoint();
0193   return cellGeomPtr(cellIndex, pos);
0194 }
0195 
0196 bool HGCalGeometry::present(const DetId& detId) const {
0197   if (detId == DetId())
0198     return false;
0199   DetId geomId = getGeometryDetId(detId);
0200   const uint32_t index(m_topology.detId2denseGeomId(geomId));
0201   return (nullptr != getGeometryRawPtr(index));
0202 }
0203 
0204 GlobalPoint HGCalGeometry::getPosition(const DetId& detid, bool debug) const {
0205   unsigned int cellIndex = indexFor(detid);
0206   GlobalPoint glob;
0207   unsigned int maxSize = (m_topology.tileTrapezoid() ? m_cellVec2.size() : m_cellVec.size());
0208   if (cellIndex < maxSize) {
0209     HGCalTopology::DecodedDetId id = m_topology.decode(detid);
0210     std::pair<float, float> xy;
0211     if (m_topology.waferHexagon6()) {
0212       xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0213       const HepGeom::Point3D<float> lcoord(xy.first, xy.second, 0);
0214       glob = m_cellVec[cellIndex].getPosition(lcoord);
0215       if (debug)
0216         edm::LogVerbatim("HGCalGeom") << "getPosition:: index " << cellIndex << " Local " << lcoord.x() << ":"
0217                                       << lcoord.y() << " ID " << id.iCell1 << ":" << id.iSec1 << " Global " << glob;
0218     } else if (m_topology.tileTrapezoid()) {
0219       const HepGeom::Point3D<float> lcoord(0, 0, 0);
0220       glob = m_cellVec2[cellIndex].getPosition(lcoord);
0221       if (debug)
0222         edm::LogVerbatim("HGCalGeom") << "getPositionTrap:: index " << cellIndex << " Local " << lcoord.x() << ":"
0223                                       << lcoord.y() << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iCell1
0224                                       << " Global " << glob;
0225     } else {
0226       if (debug) {
0227         if (detid.det() == DetId::Forward)
0228           edm::LogVerbatim("HGCalGeom") << "getPosition for " << HFNoseDetId(detid) << " Layer " << id.iLay << " Wafer "
0229                                         << id.iSec1 << ":" << id.iSec2 << " Cell " << id.iCell1 << ":" << id.iCell2;
0230         else
0231           edm::LogVerbatim("HGCalGeom") << "getPosition for " << HGCSiliconDetId(detid) << " Layer " << id.iLay
0232                                         << " Wafer " << id.iSec1 << ":" << id.iSec2 << " Cell " << id.iCell1 << ":"
0233                                         << id.iCell2;
0234       }
0235       xy = m_topology.dddConstants().locateCell(
0236           id.zSide, id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, true, false, debug);
0237       double xx = id.zSide * xy.first;
0238       double zz = id.zSide * m_topology.dddConstants().waferZ(id.iLay, true);
0239       glob = GlobalPoint(xx, xy.second, zz);
0240       if (debug)
0241         edm::LogVerbatim("HGCalGeom") << "getPositionWafer:: index " << cellIndex << " Local " << xy.first << ":"
0242                                       << xy.second << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
0243                                       << id.iCell1 << ":" << id.iCell2 << " Global " << glob;
0244     }
0245   } else {
0246     edm::LogVerbatim("HGCalGeom") << "Cannot recognize " << std::hex << detid.rawId() << " cellIndex " << cellIndex
0247                                   << ":" << maxSize << " Type " << m_topology.tileTrapezoid();
0248   }
0249   return glob;
0250 }
0251 
0252 GlobalPoint HGCalGeometry::getWaferPosition(const DetId& detid) const {
0253   unsigned int cellIndex = indexFor(detid);
0254   GlobalPoint glob;
0255   unsigned int maxSize = (m_topology.tileTrapezoid() ? m_cellVec2.size() : m_cellVec.size());
0256   if (cellIndex < maxSize) {
0257     const HepGeom::Point3D<float> lcoord(0, 0, 0);
0258     if (m_topology.tileTrapezoid()) {
0259       glob = m_cellVec2[cellIndex].getPosition(lcoord);
0260     } else {
0261       glob = m_cellVec[cellIndex].getPosition(lcoord);
0262     }
0263 #ifdef EDM_ML_DEBUG
0264     edm::LogVerbatim("HGCalGeom") << "getPositionTrap:: ID " << std::hex << detid.rawId() << std::dec << " index "
0265                                   << cellIndex << " Global " << glob;
0266 #endif
0267   }
0268   return glob;
0269 }
0270 
0271 double HGCalGeometry::getArea(const DetId& detid) const {
0272   HGCalGeometry::CornersVec corners = getNewCorners(detid);
0273   double area(0);
0274   if (corners.size() > 1) {
0275     int n = corners.size() - 1;
0276     int j = n - 1;
0277     for (int i = 0; i < n; ++i) {
0278       area += ((corners[j].x() + corners[i].x()) * (corners[i].y() - corners[j].y()));
0279       j = i;
0280     }
0281   }
0282   return std::abs(0.5 * area);
0283 }
0284 
0285 HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& detid) const {
0286   unsigned int ncorner = ((m_det == DetId::HGCalHSc) ? FlatTrd::ncorner_ : FlatHexagon::ncorner_);
0287   HGCalGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
0288   unsigned int cellIndex = indexFor(detid);
0289   HGCalTopology::DecodedDetId id = m_topology.decode(detid);
0290   if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
0291     GlobalPoint v = getPosition(detid, false);
0292     int type = std::min(id.iType, 1);
0293     std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(type, id.iSec1);
0294     float dr = k_half * (rr.second - rr.first);
0295     float dfi = m_cellVec2[cellIndex].param()[FlatTrd::k_Cell];
0296     float dz = id.zSide * m_cellVec2[cellIndex].param()[FlatTrd::k_dZ];
0297     float r = v.perp();
0298     float fi = v.phi();
0299     static const int signr[] = {1, 1, -1, -1, 1, 1, -1, -1};
0300     static const int signf[] = {-1, 1, 1, -1, -1, 1, 1, -1};
0301     static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
0302     for (unsigned int i = 0; i < ncorner; ++i) {
0303       co[i] = GlobalPoint((r + signr[i] * dr) * cos(fi + signf[i] * dfi),
0304                           (r + signr[i] * dr) * sin(fi + signf[i] * dfi),
0305                           (v.z() + signz[i] * dz));
0306     }
0307   } else if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
0308     std::pair<float, float> xy;
0309     if (m_topology.waferHexagon6()) {
0310       xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0311       float dx = m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0312       float dy = k_half * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
0313       float dz = m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0314       static const int signx[] = {0, -1, -1, 0, 1, 1, 0, -1, -1, 0, 1, 1};
0315       static const int signy[] = {-2, -1, 1, 2, 1, -1, -2, -1, 1, 2, 1, -1};
0316       static const int signz[] = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
0317       for (unsigned int i = 0; i < ncorner; ++i) {
0318         const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, signz[i] * dz);
0319         co[i] = m_cellVec[cellIndex].getPosition(lcoord);
0320       }
0321     } else {
0322       xy = m_topology.dddConstants().locateCell(
0323           id.zSide, id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true, debugLocate);
0324       float zz = m_topology.dddConstants().waferZ(id.iLay, true);
0325       float dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0326       float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
0327       float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0328       static const int signx[] = {1, -1, -2, -1, 1, 2, 1, -1, -2, -1, 1, 2};
0329       static const int signy[] = {1, 1, 0, -1, -1, 0, 1, 1, 0, -1, -1, 0};
0330       static const int signz[] = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
0331       for (unsigned int i = 0; i < ncorner; ++i) {
0332         auto xyglob = m_topology.dddConstants().localToGlobal8(
0333             id.zSide, id.iLay, id.iSec1, id.iSec2, (xy.first + signx[i] * dx), (xy.second + signy[i] * dy), true, false);
0334         double xx = id.zSide * xyglob.first;
0335         co[i] = GlobalPoint(xx, xyglob.second, id.zSide * (zz + signz[i] * dz));
0336       }
0337     }
0338   }
0339   return co;
0340 }
0341 
0342 HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& detid) const {
0343   unsigned int ncorner = FlatTrd::ncorner_;
0344   HGCalGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
0345   unsigned int cellIndex = indexFor(detid);
0346   HGCalTopology::DecodedDetId id = m_topology.decode(detid);
0347   if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
0348     GlobalPoint v = getPosition(detid, false);
0349     int type = std::min(id.iType, 1);
0350     std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(type, id.iSec1);
0351     float dr = k_half * (rr.second - rr.first);
0352     float dfi = m_cellVec2[cellIndex].param()[FlatTrd::k_Cell];
0353     float dz = id.zSide * m_cellVec2[cellIndex].param()[FlatTrd::k_dZ];
0354     float r = v.perp();
0355     float fi = v.phi();
0356     static const int signr[] = {1, 1, -1, -1, 1, 1, -1, -1};
0357     static const int signf[] = {-1, 1, 1, -1, -1, 1, 1, -1};
0358     static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
0359     for (unsigned int i = 0; i < ncorner; ++i) {
0360       co[i] = GlobalPoint((r + signr[i] * dr) * cos(fi + signf[i] * dfi),
0361                           (r + signr[i] * dr) * sin(fi + signf[i] * dfi),
0362                           (v.z() + signz[i] * dz));
0363     }
0364   } else if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
0365     std::pair<float, float> xy;
0366     float dx(0);
0367     static const int signx[] = {-1, -1, 1, 1, -1, -1, 1, 1};
0368     static const int signy[] = {-1, 1, 1, -1, -1, 1, 1, -1};
0369     static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
0370     if (m_topology.waferHexagon6()) {
0371       xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0372       dx = m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0373       float dz = m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0374       for (unsigned int i = 0; i < ncorner; ++i) {
0375         const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dx, signz[i] * dz);
0376         co[i] = m_cellVec[cellIndex].getPosition(lcoord);
0377       }
0378     } else {
0379       xy = m_topology.dddConstants().locateCell(
0380           id.zSide, id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true, debugLocate);
0381       dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0382       float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
0383       float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0384       float zz = m_topology.dddConstants().waferZ(id.iLay, true);
0385       for (unsigned int i = 0; i < ncorner; ++i) {
0386         auto xyglob = m_topology.dddConstants().localToGlobal8(
0387             id.zSide, id.iLay, id.iSec1, id.iSec2, (xy.first + signx[i] * dx), (xy.second + signy[i] * dy), true, false);
0388         double xx = id.zSide * xyglob.first;
0389         co[i] = GlobalPoint(xx, xyglob.second, id.zSide * (zz + signz[i] * dz));
0390       }
0391     }
0392   }
0393   return co;
0394 }
0395 
0396 HGCalGeometry::CornersVec HGCalGeometry::getNewCorners(const DetId& detid, bool debug) const {
0397   unsigned int ncorner = (m_det == DetId::HGCalHSc) ? 5 : 7;
0398   HGCalGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
0399   unsigned int cellIndex = indexFor(detid);
0400   HGCalTopology::DecodedDetId id = m_topology.decode(detid);
0401   if (debug)
0402     edm::LogVerbatim("HGCalGeom") << "NewCorners for Layer " << id.iLay << " Wafer " << id.iSec1 << ":" << id.iSec2
0403                                   << " Cell " << id.iCell1 << ":" << id.iCell2;
0404   if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
0405     GlobalPoint v = getPosition(detid, false);
0406     int type = std::min(id.iType, 1);
0407     std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(type, id.iSec1);
0408     float dr = k_half * (rr.second - rr.first);
0409     float dfi = m_cellVec2[cellIndex].param()[FlatTrd::k_Cell];
0410     float dz = -id.zSide * m_cellVec2[cellIndex].param()[FlatTrd::k_dZ];
0411     float r = v.perp();
0412     float fi = v.phi();
0413     static const int signr[] = {1, 1, -1, -1};
0414     static const int signf[] = {-1, 1, 1, -1};
0415     for (unsigned int i = 0; i < ncorner - 1; ++i) {
0416       co[i] = GlobalPoint(
0417           (r + signr[i] * dr) * cos(fi + signf[i] * dfi), (r + signr[i] * dr) * sin(fi + signf[i] * dfi), (v.z() + dz));
0418     }
0419     // Used to pass downstream the thickness of this cell
0420     co[ncorner - 1] = GlobalPoint(0, 0, -2 * dz);
0421   } else if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
0422     std::pair<float, float> xy;
0423     float dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0424     float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
0425     float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0426     static const int signx[] = {1, -1, -2, -1, 1, 2};
0427     static const int signy[] = {1, 1, 0, -1, -1, 0};
0428 #ifdef EDM_ML_DEBUG
0429     if (debug)
0430       edm::LogVerbatim("HGCalGeom") << "kfac " << k_fac1 << ":" << k_fac2 << " dx:dy:dz " << dx << ":" << dy << ":"
0431                                     << dz;
0432 #endif
0433     if (m_topology.waferHexagon6()) {
0434       xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0435       for (unsigned int i = 0; i < ncorner - 1; ++i) {
0436         const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, dz);
0437         co[i] = m_cellVec[cellIndex].getPosition(lcoord);
0438       }
0439     } else {
0440       xy = m_topology.dddConstants().locateCell(
0441           id.zSide, id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true, debug);
0442       float zz = m_topology.dddConstants().waferZ(id.iLay, true);
0443       for (unsigned int i = 0; i < ncorner; ++i) {
0444         double xloc = xy.first + signx[i] * dx;
0445         double yloc = xy.second + signy[i] * dy;
0446 #ifdef EDM_ML_DEBUG
0447         if (debug)
0448           edm::LogVerbatim("HGCalGeom") << "Corner " << i << " x " << xy.first << ":" << xloc << " y " << xy.second
0449                                         << ":" << yloc << " z " << zz << ":" << id.zSide * (zz + dz);
0450 #endif
0451         auto xyglob =
0452             m_topology.dddConstants().localToGlobal8(id.zSide, id.iLay, id.iSec1, id.iSec2, xloc, yloc, true, debug);
0453         double xx = id.zSide * xyglob.first;
0454         co[i] = GlobalPoint(xx, xyglob.second, id.zSide * (zz + dz));
0455       }
0456     }
0457     // Used to pass downstream the thickness of this cell
0458     co[ncorner - 1] = GlobalPoint(0, 0, -2 * dz);
0459   }
0460   return co;
0461 }
0462 
0463 DetId HGCalGeometry::neighborZ(const DetId& idin, const GlobalVector& momentum) const {
0464   DetId idnew;
0465   HGCalTopology::DecodedDetId id = m_topology.decode(idin);
0466   int lay = ((momentum.z() * id.zSide > 0) ? (id.iLay + 1) : (id.iLay - 1));
0467 #ifdef EDM_ML_DEBUG
0468   edm::LogVerbatim("HGCalGeom") << "neighborz1:: ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
0469                                 << id.iCell1 << ":" << id.iCell2 << " New Layer " << lay << " Range "
0470                                 << m_topology.dddConstants().firstLayer() << ":"
0471                                 << m_topology.dddConstants().lastLayer(true) << " pz " << momentum.z();
0472 #endif
0473   if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
0474       (momentum.z() != 0.0)) {
0475     GlobalPoint v = getPosition(idin, false);
0476     double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
0477     double grad = (z - v.z()) / momentum.z();
0478     GlobalPoint p(v.x() + grad * momentum.x(), v.y() + grad * momentum.y(), z);
0479     double r = p.perp();
0480     auto rlimit = topology().dddConstants().rangeR(z, true);
0481     if (r >= rlimit.first && r <= rlimit.second)
0482       idnew = getClosestCell(p);
0483 #ifdef EDM_ML_DEBUG
0484     edm::LogVerbatim("HGCalGeom") << "neighborz1:: Position " << v << " New Z " << z << ":" << grad << " new position "
0485                                   << p << " r-limit " << rlimit.first << ":" << rlimit.second;
0486 #endif
0487   }
0488   return idnew;
0489 }
0490 
0491 DetId HGCalGeometry::neighborZ(const DetId& idin,
0492                                const MagneticField* bField,
0493                                int charge,
0494                                const GlobalVector& momentum) const {
0495   DetId idnew;
0496   HGCalTopology::DecodedDetId id = m_topology.decode(idin);
0497   int lay = ((momentum.z() * id.zSide > 0) ? (id.iLay + 1) : (id.iLay - 1));
0498 #ifdef EDM_ML_DEBUG
0499   edm::LogVerbatim("HGCalGeom") << "neighborz2:: ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
0500                                 << id.iCell1 << ":" << id.iCell2 << " New Layer " << lay << " Range "
0501                                 << m_topology.dddConstants().firstLayer() << ":"
0502                                 << m_topology.dddConstants().lastLayer(true) << " pz " << momentum.z();
0503 #endif
0504   if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
0505       (momentum.z() != 0.0)) {
0506     GlobalPoint v = getPosition(idin, false);
0507     double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
0508     FreeTrajectoryState fts(v, momentum, charge, bField);
0509     Plane::PlanePointer nPlane = Plane::build(Plane::PositionType(0, 0, z), Plane::RotationType());
0510     AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
0511     TrajectoryStateOnSurface tsos = myAP.propagate(fts, *nPlane);
0512     GlobalPoint p;
0513     auto rlimit = topology().dddConstants().rangeR(z, true);
0514     if (tsos.isValid()) {
0515       p = tsos.globalPosition();
0516       double r = p.perp();
0517       if (r >= rlimit.first && r <= rlimit.second)
0518         idnew = getClosestCell(p);
0519     }
0520 #ifdef EDM_ML_DEBUG
0521     edm::LogVerbatim("HGCalGeom") << "neighborz2:: Position " << v << " New Z " << z << ":" << charge << ":"
0522                                   << tsos.isValid() << " new position " << p << " r limits " << rlimit.first << ":"
0523                                   << rlimit.second;
0524 #endif
0525   }
0526   return idnew;
0527 }
0528 
0529 DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const {
0530   unsigned int cellIndex = getClosestCellIndex(r);
0531   if ((cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) ||
0532       (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc)) {
0533     HGCalTopology::DecodedDetId id = m_topology.decode(m_validGeomIds[cellIndex]);
0534     if (id.det == 0)
0535       id.det = static_cast<int>(m_topology.detector());
0536     HepGeom::Point3D<float> local;
0537     if (r.z() > 0) {
0538       local = HepGeom::Point3D<float>(r.x(), r.y(), 0);
0539       id.zSide = 1;
0540     } else {
0541       local = HepGeom::Point3D<float>(-r.x(), r.y(), 0);
0542       id.zSide = -1;
0543     }
0544     if (m_topology.waferHexagon6()) {
0545       const auto& kxy = m_topology.dddConstants().assignCell(local.x(), local.y(), id.iLay, id.iType, true);
0546       id.iCell1 = kxy.second;
0547       id.iSec1 = kxy.first;
0548       id.iType = m_topology.dddConstants().waferTypeT(kxy.first);
0549       if (id.iType != 1)
0550         id.iType = -1;
0551     } else if (m_topology.tileTrapezoid()) {
0552       id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
0553       const auto& kxy = m_topology.dddConstants().assignCellTrap(r.x(), r.y(), r.z(), id.iLay, true);
0554       id.iSec1 = kxy[0];
0555       id.iCell1 = kxy[1];
0556       id.iType = kxy[2];
0557     } else {
0558       id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
0559       int zside = (r.z() > 0) ? 1 : -1;
0560 #ifdef EDM_ML_DEBUG
0561       edm::LogVerbatim("HGCalGeom") << "ZZ " << r.z() << ":" << zside << " Layer " << id.iLay << " Global " << r
0562                                     << "  Local " << local;
0563 #endif
0564       const auto& kxy =
0565           m_topology.dddConstants().assignCellHex(local.x(), local.y(), zside, id.iLay, true, false, true);
0566       id.iSec1 = kxy[0];
0567       id.iSec2 = kxy[1];
0568       id.iType = kxy[2];
0569       id.iCell1 = kxy[3];
0570       id.iCell2 = kxy[4];
0571     }
0572 #ifdef EDM_ML_DEBUG
0573     edm::LogVerbatim("HGCalGeom") << "getClosestCell: local " << local << " Id " << id.det << ":" << id.zSide << ":"
0574                                   << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":" << id.iType << ":"
0575                                   << id.iCell1 << ":" << id.iCell2;
0576 #endif
0577 
0578     //check if returned cell is valid
0579     if (id.iCell1 >= 0)
0580       return m_topology.encode(id);
0581   }
0582 
0583   //if not valid or out of bounds return a null DetId
0584   return DetId();
0585 }
0586 
0587 DetId HGCalGeometry::getClosestCellHex(const GlobalPoint& r, bool extend) const {
0588   unsigned int cellIndex = getClosestCellIndex(r);
0589   if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
0590     HGCalTopology::DecodedDetId id = m_topology.decode(m_validGeomIds[cellIndex]);
0591     if (id.det == 0)
0592       id.det = static_cast<int>(m_topology.detector());
0593     HepGeom::Point3D<float> local;
0594     if (r.z() > 0) {
0595       local = HepGeom::Point3D<float>(r.x(), r.y(), 0);
0596       id.zSide = 1;
0597     } else {
0598       local = HepGeom::Point3D<float>(-r.x(), r.y(), 0);
0599       id.zSide = -1;
0600     }
0601     if (m_topology.waferHexagon8()) {
0602       id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
0603       int zside = (r.z() > 0) ? 1 : -1;
0604 #ifdef EDM_ML_DEBUG
0605       edm::LogVerbatim("HGCalGeom") << "ZZ " << r.z() << ":" << zside << " Layer " << id.iLay << " Global " << r
0606                                     << "  Local " << local;
0607 #endif
0608       const auto& kxy =
0609           m_topology.dddConstants().assignCellHex(local.x(), local.y(), zside, id.iLay, true, extend, true);
0610       id.iSec1 = kxy[0];
0611       id.iSec2 = kxy[1];
0612       id.iType = kxy[2];
0613       id.iCell1 = kxy[3];
0614       id.iCell2 = kxy[4];
0615     }
0616 #ifdef EDM_ML_DEBUG
0617     edm::LogVerbatim("HGCalGeom") << "getClosestCell: local " << local << " Id " << id.det << ":" << id.zSide << ":"
0618                                   << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":" << id.iType << ":"
0619                                   << id.iCell1 << ":" << id.iCell2;
0620 #endif
0621 
0622     //check if returned cell is valid
0623     if (id.iCell1 >= 0)
0624       return m_topology.encode(id);
0625   }
0626 
0627   //if not valid or out of bounds return a null DetId
0628   return DetId();
0629 }
0630 
0631 HGCalGeometry::DetIdSet HGCalGeometry::getCells(const GlobalPoint& r, double dR) const {
0632   HGCalGeometry::DetIdSet dss;
0633   return dss;
0634 }
0635 
0636 std::string HGCalGeometry::cellElement() const {
0637   if (m_subdet == HGCEE || m_det == DetId::HGCalEE)
0638     return "HGCalEE";
0639   else if (m_subdet == HGCHEF || m_det == DetId::HGCalHSi)
0640     return "HGCalHEFront";
0641   else if (m_subdet == HGCHEB || m_det == DetId::HGCalHSc)
0642     return "HGCalHEBack";
0643   else
0644     return "Unknown";
0645 }
0646 
0647 unsigned int HGCalGeometry::indexFor(const DetId& detId) const {
0648   unsigned int cellIndex = ((m_det == DetId::HGCalHSc) ? m_cellVec2.size() : m_cellVec.size());
0649   if (detId != DetId()) {
0650     DetId geomId = getGeometryDetId(detId);
0651     cellIndex = m_topology.detId2denseGeomId(geomId);
0652 #ifdef EDM_ML_DEBUG
0653     edm::LogVerbatim("HGCalGeom") << "indexFor " << std::hex << detId.rawId() << ":" << geomId.rawId() << std::dec
0654                                   << " index " << cellIndex;
0655 #endif
0656   }
0657   return cellIndex;
0658 }
0659 
0660 unsigned int HGCalGeometry::sizeForDenseIndex() const { return m_topology.totalGeomModules(); }
0661 
0662 const CaloCellGeometry* HGCalGeometry::getGeometryRawPtr(uint32_t index) const {
0663   // Modify the RawPtr class
0664   if (m_det == DetId::HGCalHSc) {
0665     if (m_cellVec2.size() < index)
0666       return nullptr;
0667     const CaloCellGeometry* cell(&m_cellVec2[index]);
0668     return (nullptr == cell->param() ? nullptr : cell);
0669   } else {
0670     if (m_cellVec.size() < index)
0671       return nullptr;
0672     const CaloCellGeometry* cell(&m_cellVec[index]);
0673     return (nullptr == cell->param() ? nullptr : cell);
0674   }
0675 }
0676 
0677 std::shared_ptr<const CaloCellGeometry> HGCalGeometry::cellGeomPtr(uint32_t index) const {
0678   if ((index >= m_cellVec.size() && m_det != DetId::HGCalHSc) ||
0679       (index >= m_cellVec2.size() && m_det == DetId::HGCalHSc) || (m_validGeomIds[index].rawId() == 0))
0680     return nullptr;
0681   static const auto do_not_delete = [](const void*) {};
0682   if (m_det == DetId::HGCalHSc) {
0683     auto cell = std::shared_ptr<const CaloCellGeometry>(&m_cellVec2[index], do_not_delete);
0684     if (nullptr == cell->param())
0685       return nullptr;
0686     return cell;
0687   } else {
0688     auto cell = std::shared_ptr<const CaloCellGeometry>(&m_cellVec[index], do_not_delete);
0689     if (nullptr == cell->param())
0690       return nullptr;
0691     return cell;
0692   }
0693 }
0694 
0695 std::shared_ptr<const CaloCellGeometry> HGCalGeometry::cellGeomPtr(uint32_t index, const GlobalPoint& pos) const {
0696   if ((index >= m_cellVec.size() && m_det != DetId::HGCalHSc) ||
0697       (index >= m_cellVec2.size() && m_det == DetId::HGCalHSc) || (m_validGeomIds[index].rawId() == 0))
0698     return nullptr;
0699   if (pos == GlobalPoint())
0700     return cellGeomPtr(index);
0701   if (m_det == DetId::HGCalHSc) {
0702     auto cell = std::make_shared<FlatTrd>(m_cellVec2[index]);
0703     cell->setPosition(pos);
0704 #ifdef EDM_ML_DEBUG
0705     edm::LogVerbatim("HGCalGeom") << "cellGeomPtr " << index << ":" << cell;
0706 #endif
0707     if (nullptr == cell->param())
0708       return nullptr;
0709     return cell;
0710   } else {
0711     auto cell = std::make_shared<FlatHexagon>(m_cellVec[index]);
0712     cell->setPosition(pos);
0713 #ifdef EDM_ML_DEBUG
0714     edm::LogVerbatim("HGCalGeom") << "cellGeomPtr " << index << ":" << cell;
0715 #endif
0716     if (nullptr == cell->param())
0717       return nullptr;
0718     return cell;
0719   }
0720 }
0721 
0722 void HGCalGeometry::addValidID(const DetId& id) {
0723   edm::LogError("HGCalGeom") << "HGCalGeometry::addValidID is not implemented";
0724 }
0725 
0726 unsigned int HGCalGeometry::getClosestCellIndex(const GlobalPoint& r) const {
0727   return ((m_det == DetId::HGCalHSc) ? getClosestCellIndex(r, m_cellVec2) : getClosestCellIndex(r, m_cellVec));
0728 }
0729 
0730 template <class T>
0731 unsigned int HGCalGeometry::getClosestCellIndex(const GlobalPoint& r, const std::vector<T>& vec) const {
0732   float phip = r.phi();
0733   float zp = r.z();
0734   float dzmin(9999), dphimin(9999), dphi10(0.175);
0735   unsigned int cellIndex = vec.size();
0736   for (unsigned int k = 0; k < vec.size(); ++k) {
0737     float dphi = phip - vec[k].phiPos();
0738     while (dphi > M_PI)
0739       dphi -= 2 * M_PI;
0740     while (dphi <= -M_PI)
0741       dphi += 2 * M_PI;
0742     if (std::abs(dphi) < dphi10) {
0743       float dz = std::abs(zp - vec[k].getPosition().z());
0744       if (dz < (dzmin + 0.001)) {
0745         dzmin = dz;
0746         if (std::abs(dphi) < (dphimin + 0.01)) {
0747           cellIndex = k;
0748           dphimin = std::abs(dphi);
0749         } else {
0750           if (cellIndex >= vec.size())
0751             cellIndex = k;
0752         }
0753       }
0754     }
0755   }
0756 #ifdef EDM_ML_DEBUG
0757   edm::LogVerbatim("HGCalGeom") << "getClosestCellIndex::Input " << zp << ":" << phip << " Index " << cellIndex;
0758   if (cellIndex < vec.size())
0759     edm::LogVerbatim("HGCalGeom") << " Cell z " << vec[cellIndex].getPosition().z() << ":" << dzmin << " phi "
0760                                   << vec[cellIndex].phiPos() << ":" << dphimin;
0761 #endif
0762   return cellIndex;
0763 }
0764 
0765 // FIXME: Change sorting algorithm if needed
0766 namespace {
0767   struct rawIdSort {
0768     bool operator()(const DetId& a, const DetId& b) { return (a.rawId() < b.rawId()); }
0769   };
0770 }  // namespace
0771 
0772 void HGCalGeometry::sortDetIds(void) {
0773   m_validIds.shrink_to_fit();
0774   std::sort(m_validIds.begin(), m_validIds.end(), rawIdSort());
0775 }
0776 
0777 void HGCalGeometry::getSummary(CaloSubdetectorGeometry::TrVec& trVector,
0778                                CaloSubdetectorGeometry::IVec& iVector,
0779                                CaloSubdetectorGeometry::DimVec& dimVector,
0780                                CaloSubdetectorGeometry::IVec& dinsVector) const {
0781   unsigned int numberOfCells = m_topology.totalGeomModules();  // total Geom Modules both sides
0782   unsigned int numberOfShapes = k_NumberOfShapes;
0783   unsigned int numberOfParametersPerShape = ((m_det == DetId::HGCalHSc) ? (unsigned int)(k_NumberOfParametersPerTrd)
0784                                                                         : (unsigned int)(k_NumberOfParametersPerHex));
0785 
0786   trVector.reserve(numberOfCells * numberOfTransformParms());
0787   iVector.reserve(numberOfCells);
0788   dimVector.reserve(numberOfShapes * numberOfParametersPerShape);
0789   dinsVector.reserve(numberOfCells);
0790 
0791   for (unsigned itr = 0; itr < m_topology.dddConstants().getTrFormN(); ++itr) {
0792     HGCalParameters::hgtrform mytr = m_topology.dddConstants().getTrForm(itr);
0793     int layer = mytr.lay;
0794 
0795     if (m_topology.waferHexagon6()) {
0796       for (int wafer = 0; wafer < m_topology.dddConstants().sectors(); ++wafer) {
0797         if (m_topology.dddConstants().waferInLayer(wafer, layer, true)) {
0798           HGCalParameters::hgtrap vol = m_topology.dddConstants().getModule(wafer, true, true);
0799           ParmVec params(numberOfParametersPerShape, 0);
0800           params[FlatHexagon::k_dZ] = vol.dz;
0801           params[FlatHexagon::k_r] = vol.cellSize;
0802           params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
0803           dimVector.insert(dimVector.end(), params.begin(), params.end());
0804         }
0805       }
0806     } else if (m_topology.tileTrapezoid()) {
0807       int indx = m_topology.dddConstants().layerIndex(layer, true);
0808       for (int md = m_topology.dddConstants().getParameter()->firstModule_[indx];
0809            md <= m_topology.dddConstants().getParameter()->lastModule_[indx];
0810            ++md) {
0811         HGCalParameters::hgtrap vol = m_topology.dddConstants().getModule(md, true, true);
0812         ParmVec params(numberOfParametersPerShape, 0);
0813         params[FlatTrd::k_dZ] = vol.dz;
0814         params[FlatTrd::k_Theta] = params[FlatTrd::k_Phi] = 0;
0815         params[FlatTrd::k_dY1] = params[FlatTrd::k_dY2] = vol.h;
0816         params[FlatTrd::k_dX1] = params[FlatTrd::k_dX3] = vol.bl;
0817         params[FlatTrd::k_dX2] = params[FlatTrd::k_dX4] = vol.tl;
0818         params[FlatTrd::k_Alp1] = params[FlatTrd::k_Alp2] = vol.alpha;
0819         params[FlatTrd::k_Cell] = vol.cellSize;
0820         dimVector.insert(dimVector.end(), params.begin(), params.end());
0821       }
0822     } else {
0823       for (int wafer = 0; wafer < m_topology.dddConstants().sectors(); ++wafer) {
0824         if (m_topology.dddConstants().waferInLayer(wafer, layer, true)) {
0825           HGCalParameters::hgtrap vol = m_topology.dddConstants().getModule(wafer, true, true);
0826           ParmVec params(numberOfParametersPerShape, 0);
0827           params[FlatHexagon::k_dZ] = vol.dz;
0828           params[FlatHexagon::k_r] = vol.cellSize;
0829           params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
0830           dimVector.insert(dimVector.end(), params.begin(), params.end());
0831         }
0832       }
0833     }
0834   }
0835 
0836   for (unsigned int i(0); i < numberOfCells; ++i) {
0837     DetId detId = m_validGeomIds[i];
0838     int layer(0);
0839     if (m_topology.waferHexagon6()) {
0840       layer = HGCalDetId(detId).layer();
0841     } else if (m_topology.tileTrapezoid()) {
0842       layer = HGCScintillatorDetId(detId).layer();
0843     } else if (m_topology.isHFNose()) {
0844       layer = HFNoseDetId(detId).layer();
0845     } else {
0846       layer = HGCSiliconDetId(detId).layer();
0847     }
0848     dinsVector.emplace_back(m_topology.detId2denseGeomId(detId));
0849     iVector.emplace_back(layer);
0850 
0851     Tr3D tr;
0852     auto ptr = cellGeomPtr(i);
0853     if (nullptr != ptr) {
0854       ptr->getTransform(tr, (Pt3DVec*)nullptr);
0855 
0856       if (Tr3D() == tr) {  // there is no rotation
0857         const GlobalPoint& gp(ptr->getPosition());
0858         tr = HepGeom::Translate3D(gp.x(), gp.y(), gp.z());
0859       }
0860 
0861       const CLHEP::Hep3Vector tt(tr.getTranslation());
0862       trVector.emplace_back(tt.x());
0863       trVector.emplace_back(tt.y());
0864       trVector.emplace_back(tt.z());
0865       if (6 == numberOfTransformParms()) {
0866         const CLHEP::HepRotation rr(tr.getRotation());
0867         const ROOT::Math::Transform3D rtr(
0868             rr.xx(), rr.xy(), rr.xz(), tt.x(), rr.yx(), rr.yy(), rr.yz(), tt.y(), rr.zx(), rr.zy(), rr.zz(), tt.z());
0869         ROOT::Math::EulerAngles ea;
0870         rtr.GetRotation(ea);
0871         trVector.emplace_back(ea.Phi());
0872         trVector.emplace_back(ea.Theta());
0873         trVector.emplace_back(ea.Psi());
0874       }
0875     }
0876   }
0877 }
0878 
0879 DetId HGCalGeometry::getGeometryDetId(DetId detId) const {
0880   DetId geomId;
0881   if (m_topology.waferHexagon6()) {
0882     geomId = static_cast<DetId>(HGCalDetId(detId).geometryCell());
0883   } else if (m_topology.tileTrapezoid()) {
0884     geomId = static_cast<DetId>(HGCScintillatorDetId(detId).geometryCell());
0885   } else if (m_topology.isHFNose()) {
0886     geomId = static_cast<DetId>(HFNoseDetId(detId).geometryCell());
0887   } else {
0888     geomId = static_cast<DetId>(HGCSiliconDetId(detId).geometryCell());
0889   }
0890   return geomId;
0891 }
0892 
0893 #include "FWCore/Utilities/interface/typelookup.h"
0894 
0895 TYPELOOKUP_DATA_REG(HGCalGeometry);