Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-05-29 03:17:44

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