Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-08 23:50:39

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