Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 13:03:39

0001 /* for High Granularity Calorimeter TestBeam
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/HGCalTBGeometry.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 HGCalTBGeometry::HGCalTBGeometry(const HGCalTBTopology& topology_)
0027     : m_topology(topology_),
0028       m_validGeomIds(topology_.totalGeomModules()),
0029       m_det(topology_.detector()),
0030       m_subdet(topology_.subDetector()),
0031       twoBysqrt3_(2.0 / std::sqrt(3.0)) {
0032   m_cellVec = CellVec(topology_.totalGeomModules());
0033   m_validIds.reserve(m_topology.totalModules());
0034 #ifdef EDM_ML_DEBUG
0035   edm::LogVerbatim("HGCalGeom") << "Expected total # of Geometry Modules " << m_topology.totalGeomModules();
0036 #endif
0037 }
0038 
0039 HGCalTBGeometry::~HGCalTBGeometry() {}
0040 
0041 void HGCalTBGeometry::fillNamedParams(DDFilteredView fv) {}
0042 
0043 void HGCalTBGeometry::initializeParms() {}
0044 
0045 void HGCalTBGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int i, Pt3D& ref) {
0046   FlatHexagon::localCorners(lc, pv, ref);
0047 }
0048 
0049 void HGCalTBGeometry::newCell(
0050     const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0051   DetId geomId = getGeometryDetId(detId);
0052   int cells(0);
0053   HGCalTBTopology::DecodedDetId id = m_topology.decode(detId);
0054   cells = m_topology.dddConstants().numberCellsHexagon(id.iSec1);
0055 #ifdef EDM_ML_DEBUG
0056   edm::LogVerbatim("HGCalGeom") << "NewCell " << HGCalDetId(detId) << " GEOM " << HGCalDetId(geomId);
0057 #endif
0058   const uint32_t cellIndex(m_topology.detId2denseGeomId(geomId));
0059 
0060   m_cellVec.at(cellIndex) = FlatHexagon(cornersMgr(), f1, f2, f3, parm);
0061   m_validGeomIds.at(cellIndex) = geomId;
0062 
0063 #ifdef EDM_ML_DEBUG
0064   edm::LogVerbatim("HGCalGeom") << "Store for DetId " << std::hex << detId.rawId() << " GeomId " << geomId.rawId()
0065                                 << std::dec << " Index " << cellIndex << " cells " << cells;
0066   unsigned int nOld = m_validIds.size();
0067 #endif
0068   for (int cell = 0; cell < cells; ++cell) {
0069     id.iCell1 = cell;
0070     DetId idc = m_topology.encode(id);
0071     if (m_topology.valid(idc)) {
0072       m_validIds.emplace_back(idc);
0073 #ifdef EDM_ML_DEBUG
0074       edm::LogVerbatim("HGCalGeom") << "Valid Id [" << cell << "] " << HGCalDetId(idc);
0075 #endif
0076     }
0077   }
0078 #ifdef EDM_ML_DEBUG
0079   edm::LogVerbatim("HGCalGeom") << "HGCalTBGeometry::newCell-> [" << cellIndex << "]"
0080                                 << " front:" << f1.x() << '/' << f1.y() << '/' << f1.z() << " back:" << f2.x() << '/'
0081                                 << f2.y() << '/' << f2.z() << " eta|phi " << m_cellVec[cellIndex].etaPos() << ":"
0082                                 << m_cellVec[cellIndex].phiPos();
0083   unsigned int nNew = m_validIds.size();
0084   edm::LogVerbatim("HGCalGeom") << "ID: " << HGCalDetId(detId) << " with valid DetId from " << nOld << " to " << nNew;
0085 #endif
0086 }
0087 
0088 std::shared_ptr<const CaloCellGeometry> HGCalTBGeometry::getGeometry(const DetId& detId) const {
0089   if (detId == DetId())
0090     return nullptr;  // nothing to get
0091   DetId geomId = getGeometryDetId(detId);
0092   const uint32_t cellIndex(m_topology.detId2denseGeomId(geomId));
0093   const GlobalPoint pos = (detId != geomId) ? getPosition(detId, false) : GlobalPoint();
0094   return cellGeomPtr(cellIndex, pos);
0095 }
0096 
0097 bool HGCalTBGeometry::present(const DetId& detId) const {
0098   if (detId == DetId())
0099     return false;
0100   DetId geomId = getGeometryDetId(detId);
0101   const uint32_t index(m_topology.detId2denseGeomId(geomId));
0102   return (nullptr != getGeometryRawPtr(index));
0103 }
0104 
0105 GlobalPoint HGCalTBGeometry::getPosition(const DetId& detid, bool debug) const {
0106   unsigned int cellIndex = indexFor(detid);
0107   GlobalPoint glob;
0108   unsigned int maxSize = m_cellVec.size();
0109   if (cellIndex < maxSize) {
0110     HGCalTBTopology::DecodedDetId id = m_topology.decode(detid);
0111     std::pair<float, float> xy;
0112     xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0113     const HepGeom::Point3D<float> lcoord(xy.first, xy.second, 0);
0114     glob = m_cellVec[cellIndex].getPosition(lcoord);
0115     if (debug)
0116       edm::LogVerbatim("HGCalGeom") << "getPosition:: index " << cellIndex << " Local " << lcoord.x() << ":"
0117                                     << lcoord.y() << " ID " << id.iCell1 << ":" << id.iSec1 << " Global " << glob;
0118   }
0119   return glob;
0120 }
0121 
0122 GlobalPoint HGCalTBGeometry::getWaferPosition(const DetId& detid) const {
0123   unsigned int cellIndex = indexFor(detid);
0124   GlobalPoint glob;
0125   unsigned int maxSize = m_cellVec.size();
0126   if (cellIndex < maxSize) {
0127     const HepGeom::Point3D<float> lcoord(0, 0, 0);
0128     glob = m_cellVec[cellIndex].getPosition(lcoord);
0129 #ifdef EDM_ML_DEBUG
0130     edm::LogVerbatim("HGCalGeom") << "getWaferPosition:: ID " << std::hex << detid.rawId() << std::dec << " index "
0131                                   << cellIndex << " Global " << glob;
0132 #endif
0133   }
0134   return glob;
0135 }
0136 
0137 double HGCalTBGeometry::getArea(const DetId& detid) const {
0138   HGCalTBGeometry::CornersVec corners = getNewCorners(detid);
0139   double area(0);
0140   if (corners.size() > 1) {
0141     int n = corners.size() - 1;
0142     int j = n - 1;
0143     for (int i = 0; i < n; ++i) {
0144       area += ((corners[j].x() + corners[i].x()) * (corners[i].y() - corners[j].y()));
0145       j = i;
0146     }
0147   }
0148   return std::abs(0.5 * area);
0149 }
0150 
0151 HGCalTBGeometry::CornersVec HGCalTBGeometry::getCorners(const DetId& detid) const {
0152   unsigned int ncorner = ((m_det == DetId::HGCalHSc) ? FlatTrd::ncorner_ : FlatHexagon::ncorner_);
0153   HGCalTBGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
0154   unsigned int cellIndex = indexFor(detid);
0155   HGCalTBTopology::DecodedDetId id = m_topology.decode(detid);
0156   if (cellIndex < m_cellVec.size()) {
0157     std::pair<float, float> xy;
0158     xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0159     float dx = m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0160     float dy = k_half * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
0161     float dz = m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0162     static const int signx[] = {0, -1, -1, 0, 1, 1, 0, -1, -1, 0, 1, 1};
0163     static const int signy[] = {-2, -1, 1, 2, 1, -1, -2, -1, 1, 2, 1, -1};
0164     static const int signz[] = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
0165     for (unsigned int i = 0; i < ncorner; ++i) {
0166       const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, signz[i] * dz);
0167       co[i] = m_cellVec[cellIndex].getPosition(lcoord);
0168     }
0169   }
0170   return co;
0171 }
0172 
0173 HGCalTBGeometry::CornersVec HGCalTBGeometry::get8Corners(const DetId& detid) const {
0174   unsigned int ncorner = FlatTrd::ncorner_;
0175   HGCalTBGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
0176   unsigned int cellIndex = indexFor(detid);
0177   HGCalTBTopology::DecodedDetId id = m_topology.decode(detid);
0178   if (cellIndex < m_cellVec.size()) {
0179     std::pair<float, float> xy;
0180     float dx(0);
0181     static const int signx[] = {-1, -1, 1, 1, -1, -1, 1, 1};
0182     static const int signy[] = {-1, 1, 1, -1, -1, 1, 1, -1};
0183     static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
0184     xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0185     dx = m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0186     float dz = m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0187     for (unsigned int i = 0; i < ncorner; ++i) {
0188       const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dx, signz[i] * dz);
0189       co[i] = m_cellVec[cellIndex].getPosition(lcoord);
0190     }
0191   }
0192   return co;
0193 }
0194 
0195 HGCalTBGeometry::CornersVec HGCalTBGeometry::getNewCorners(const DetId& detid, bool debug) const {
0196   unsigned int ncorner = (m_det == DetId::HGCalHSc) ? 5 : 7;
0197   HGCalTBGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
0198   unsigned int cellIndex = indexFor(detid);
0199   HGCalTBTopology::DecodedDetId id = m_topology.decode(detid);
0200   if (debug)
0201     edm::LogVerbatim("HGCalGeom") << "NewCorners for Layer " << id.iLay << " Wafer " << id.iSec1 << ":" << id.iSec2
0202                                   << " Cell " << id.iCell1 << ":" << id.iCell2;
0203   if (cellIndex < m_cellVec.size()) {
0204     std::pair<float, float> xy;
0205     float dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
0206     float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
0207     float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
0208     static const int signx[] = {1, -1, -2, -1, 1, 2};
0209     static const int signy[] = {1, 1, 0, -1, -1, 0};
0210 #ifdef EDM_ML_DEBUG
0211     if (debug)
0212       edm::LogVerbatim("HGCalGeom") << "kfac " << k_fac1 << ":" << k_fac2 << " dx:dy:dz " << dx << ":" << dy << ":"
0213                                     << dz;
0214 #endif
0215     xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
0216     for (unsigned int i = 0; i < ncorner - 1; ++i) {
0217       const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, dz);
0218       co[i] = m_cellVec[cellIndex].getPosition(lcoord);
0219     }
0220     // Used to pass downstream the thickness of this cell
0221     co[ncorner - 1] = GlobalPoint(0, 0, -2 * dz);
0222   }
0223   return co;
0224 }
0225 
0226 DetId HGCalTBGeometry::neighborZ(const DetId& idin, const GlobalVector& momentum) const {
0227   DetId idnew;
0228   HGCalTBTopology::DecodedDetId id = m_topology.decode(idin);
0229   int lay = ((momentum.z() * id.zSide > 0) ? (id.iLay + 1) : (id.iLay - 1));
0230 #ifdef EDM_ML_DEBUG
0231   edm::LogVerbatim("HGCalGeom") << "neighborz1:: ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
0232                                 << id.iCell1 << ":" << id.iCell2 << " New Layer " << lay << " Range "
0233                                 << m_topology.dddConstants().firstLayer() << ":"
0234                                 << m_topology.dddConstants().lastLayer(true) << " pz " << momentum.z();
0235 #endif
0236   if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
0237       (momentum.z() != 0.0)) {
0238     GlobalPoint v = getPosition(idin, false);
0239     double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
0240     double grad = (z - v.z()) / momentum.z();
0241     GlobalPoint p(v.x() + grad * momentum.x(), v.y() + grad * momentum.y(), z);
0242     double r = p.perp();
0243     auto rlimit = topology().dddConstants().rangeR(z, true);
0244     if (r >= rlimit.first && r <= rlimit.second)
0245       idnew = getClosestCell(p);
0246 #ifdef EDM_ML_DEBUG
0247     edm::LogVerbatim("HGCalGeom") << "neighborz1:: Position " << v << " New Z " << z << ":" << grad << " new position "
0248                                   << p << " r-limit " << rlimit.first << ":" << rlimit.second;
0249 #endif
0250   }
0251   return idnew;
0252 }
0253 
0254 DetId HGCalTBGeometry::neighborZ(const DetId& idin,
0255                                  const MagneticField* bField,
0256                                  int charge,
0257                                  const GlobalVector& momentum) const {
0258   DetId idnew;
0259   HGCalTBTopology::DecodedDetId id = m_topology.decode(idin);
0260   int lay = ((momentum.z() * id.zSide > 0) ? (id.iLay + 1) : (id.iLay - 1));
0261 #ifdef EDM_ML_DEBUG
0262   edm::LogVerbatim("HGCalGeom") << "neighborz2:: ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
0263                                 << id.iCell1 << ":" << id.iCell2 << " New Layer " << lay << " Range "
0264                                 << m_topology.dddConstants().firstLayer() << ":"
0265                                 << m_topology.dddConstants().lastLayer(true) << " pz " << momentum.z();
0266 #endif
0267   if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
0268       (momentum.z() != 0.0)) {
0269     GlobalPoint v = getPosition(idin, false);
0270     double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
0271     FreeTrajectoryState fts(v, momentum, charge, bField);
0272     Plane::PlanePointer nPlane = Plane::build(Plane::PositionType(0, 0, z), Plane::RotationType());
0273     AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
0274     TrajectoryStateOnSurface tsos = myAP.propagate(fts, *nPlane);
0275     GlobalPoint p;
0276     auto rlimit = topology().dddConstants().rangeR(z, true);
0277     if (tsos.isValid()) {
0278       p = tsos.globalPosition();
0279       double r = p.perp();
0280       if (r >= rlimit.first && r <= rlimit.second)
0281         idnew = getClosestCell(p);
0282     }
0283 #ifdef EDM_ML_DEBUG
0284     edm::LogVerbatim("HGCalGeom") << "neighborz2:: Position " << v << " New Z " << z << ":" << charge << ":"
0285                                   << tsos.isValid() << " new position " << p << " r limits " << rlimit.first << ":"
0286                                   << rlimit.second;
0287 #endif
0288   }
0289   return idnew;
0290 }
0291 
0292 DetId HGCalTBGeometry::getClosestCell(const GlobalPoint& r) const {
0293   unsigned int cellIndex = getClosestCellIndex(r);
0294   if (cellIndex < m_cellVec.size()) {
0295     HGCalTBTopology::DecodedDetId id = m_topology.decode(m_validGeomIds[cellIndex]);
0296     if (id.det == 0)
0297       id.det = static_cast<int>(m_topology.detector());
0298     HepGeom::Point3D<float> local;
0299     if (r.z() > 0) {
0300       local = HepGeom::Point3D<float>(r.x(), r.y(), 0);
0301       id.zSide = 1;
0302     } else {
0303       local = HepGeom::Point3D<float>(-r.x(), r.y(), 0);
0304       id.zSide = -1;
0305     }
0306     const auto& kxy = m_topology.dddConstants().assignCell(local.x(), local.y(), id.iLay, id.iType, true);
0307     id.iCell1 = kxy.second;
0308     id.iSec1 = kxy.first;
0309     id.iType = m_topology.dddConstants().waferTypeT(kxy.first);
0310     if (id.iType != 1)
0311       id.iType = -1;
0312 #ifdef EDM_ML_DEBUG
0313     edm::LogVerbatim("HGCalGeom") << "getClosestCell: local " << local << " Id " << id.det << ":" << id.zSide << ":"
0314                                   << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":" << id.iType << ":"
0315                                   << id.iCell1 << ":" << id.iCell2;
0316 #endif
0317     //check if returned cell is valid
0318     if (id.iCell1 >= 0)
0319       return m_topology.encode(id);
0320   }
0321 
0322   //if not valid or out of bounds return a null DetId
0323   return DetId();
0324 }
0325 
0326 HGCalTBGeometry::DetIdSet HGCalTBGeometry::getCells(const GlobalPoint& r, double dR) const {
0327   HGCalTBGeometry::DetIdSet dss;
0328   return dss;
0329 }
0330 
0331 std::string HGCalTBGeometry::cellElement() const {
0332   if (m_subdet == HGCEE)
0333     return "HGCalEE";
0334   else if (m_subdet == HGCHEF)
0335     return "HGCalHEFront";
0336   else
0337     return "Unknown";
0338 }
0339 
0340 unsigned int HGCalTBGeometry::indexFor(const DetId& detId) const {
0341   unsigned int cellIndex = m_cellVec.size();
0342   if (detId != DetId()) {
0343     DetId geomId = getGeometryDetId(detId);
0344     cellIndex = m_topology.detId2denseGeomId(geomId);
0345 #ifdef EDM_ML_DEBUG
0346     edm::LogVerbatim("HGCalGeom") << "indexFor " << std::hex << detId.rawId() << ":" << geomId.rawId() << std::dec
0347                                   << " index " << cellIndex;
0348 #endif
0349   }
0350   return cellIndex;
0351 }
0352 
0353 unsigned int HGCalTBGeometry::sizeForDenseIndex() const { return m_topology.totalGeomModules(); }
0354 
0355 const CaloCellGeometry* HGCalTBGeometry::getGeometryRawPtr(uint32_t index) const {
0356   // Modify the RawPtr class
0357   if (m_cellVec.size() < index)
0358     return nullptr;
0359   const CaloCellGeometry* cell(&m_cellVec[index]);
0360   return (nullptr == cell->param() ? nullptr : cell);
0361 }
0362 
0363 std::shared_ptr<const CaloCellGeometry> HGCalTBGeometry::cellGeomPtr(uint32_t index) const {
0364   if (index >= m_cellVec.size())
0365     return nullptr;
0366   static const auto do_not_delete = [](const void*) {};
0367   auto cell = std::shared_ptr<const CaloCellGeometry>(&m_cellVec[index], do_not_delete);
0368   if (nullptr == cell->param())
0369     return nullptr;
0370   return cell;
0371 }
0372 
0373 std::shared_ptr<const CaloCellGeometry> HGCalTBGeometry::cellGeomPtr(uint32_t index, const GlobalPoint& pos) const {
0374   if ((index >= m_cellVec.size()) || (m_validGeomIds[index].rawId() == 0))
0375     return nullptr;
0376   if (pos == GlobalPoint())
0377     return cellGeomPtr(index);
0378   auto cell = std::make_shared<FlatHexagon>(m_cellVec[index]);
0379   cell->setPosition(pos);
0380 #ifdef EDM_ML_DEBUG
0381   edm::LogVerbatim("HGCalGeom") << "cellGeomPtr " << index << ":" << cell;
0382 #endif
0383   if (nullptr == cell->param())
0384     return nullptr;
0385   return cell;
0386 }
0387 
0388 void HGCalTBGeometry::addValidID(const DetId& id) {
0389   edm::LogError("HGCalGeom") << "HGCalTBGeometry::addValidID is not implemented";
0390 }
0391 
0392 unsigned int HGCalTBGeometry::getClosestCellIndex(const GlobalPoint& r) const {
0393   return (getClosestCellIndex(r, m_cellVec));
0394 }
0395 
0396 template <class T>
0397 unsigned int HGCalTBGeometry::getClosestCellIndex(const GlobalPoint& r, const std::vector<T>& vec) const {
0398   float phip = r.phi();
0399   float zp = r.z();
0400   float dzmin(9999), dphimin(9999), dphi10(0.175);
0401   unsigned int cellIndex = vec.size();
0402   for (unsigned int k = 0; k < vec.size(); ++k) {
0403     float dphi = phip - vec[k].phiPos();
0404     while (dphi > M_PI)
0405       dphi -= 2 * M_PI;
0406     while (dphi <= -M_PI)
0407       dphi += 2 * M_PI;
0408     if (std::abs(dphi) < dphi10) {
0409       float dz = std::abs(zp - vec[k].getPosition().z());
0410       if (dz < (dzmin + 0.001)) {
0411         dzmin = dz;
0412         if (std::abs(dphi) < (dphimin + 0.01)) {
0413           cellIndex = k;
0414           dphimin = std::abs(dphi);
0415         } else {
0416           if (cellIndex >= vec.size())
0417             cellIndex = k;
0418         }
0419       }
0420     }
0421   }
0422 #ifdef EDM_ML_DEBUG
0423   edm::LogVerbatim("HGCalGeom") << "getClosestCellIndex::Input " << zp << ":" << phip << " Index " << cellIndex;
0424   if (cellIndex < vec.size())
0425     edm::LogVerbatim("HGCalGeom") << " Cell z " << vec[cellIndex].getPosition().z() << ":" << dzmin << " phi "
0426                                   << vec[cellIndex].phiPos() << ":" << dphimin;
0427 #endif
0428   return cellIndex;
0429 }
0430 
0431 // FIXME: Change sorting algorithm if needed
0432 namespace {
0433   struct rawIdSort {
0434     bool operator()(const DetId& a, const DetId& b) { return (a.rawId() < b.rawId()); }
0435   };
0436 }  // namespace
0437 
0438 void HGCalTBGeometry::sortDetIds(void) {
0439   m_validIds.shrink_to_fit();
0440   std::sort(m_validIds.begin(), m_validIds.end(), rawIdSort());
0441 }
0442 
0443 void HGCalTBGeometry::getSummary(CaloSubdetectorGeometry::TrVec& trVector,
0444                                  CaloSubdetectorGeometry::IVec& iVector,
0445                                  CaloSubdetectorGeometry::DimVec& dimVector,
0446                                  CaloSubdetectorGeometry::IVec& dinsVector) const {
0447   unsigned int numberOfCells = m_topology.totalGeomModules();  // total Geom Modules both sides
0448   unsigned int numberOfShapes = k_NumberOfShapes;
0449   unsigned int numberOfParametersPerShape = ((m_det == DetId::HGCalHSc) ? (unsigned int)(k_NumberOfParametersPerTrd)
0450                                                                         : (unsigned int)(k_NumberOfParametersPerHex));
0451 
0452   trVector.reserve(numberOfCells * numberOfTransformParms());
0453   iVector.reserve(numberOfCells);
0454   dimVector.reserve(numberOfShapes * numberOfParametersPerShape);
0455   dinsVector.reserve(numberOfCells);
0456 
0457   for (unsigned itr = 0; itr < m_topology.dddConstants().getTrFormN(); ++itr) {
0458     HGCalTBParameters::hgtrform mytr = m_topology.dddConstants().getTrForm(itr);
0459     int layer = mytr.lay;
0460 
0461     for (int wafer = 0; wafer < m_topology.dddConstants().sectors(); ++wafer) {
0462       if (m_topology.dddConstants().waferInLayer(wafer, layer, true)) {
0463         HGCalTBParameters::hgtrap vol = m_topology.dddConstants().getModule(wafer, true, true);
0464         ParmVec params(numberOfParametersPerShape, 0);
0465         params[FlatHexagon::k_dZ] = vol.dz;
0466         params[FlatHexagon::k_r] = vol.cellSize;
0467         params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
0468         dimVector.insert(dimVector.end(), params.begin(), params.end());
0469       }
0470     }
0471   }
0472 
0473   for (unsigned int i(0); i < numberOfCells; ++i) {
0474     DetId detId = m_validGeomIds[i];
0475     int layer = HGCalDetId(detId).layer();
0476     dinsVector.emplace_back(m_topology.detId2denseGeomId(detId));
0477     iVector.emplace_back(layer);
0478 
0479     Tr3D tr;
0480     auto ptr = cellGeomPtr(i);
0481     if (nullptr != ptr) {
0482       ptr->getTransform(tr, (Pt3DVec*)nullptr);
0483 
0484       if (Tr3D() == tr) {  // there is no rotation
0485         const GlobalPoint& gp(ptr->getPosition());
0486         tr = HepGeom::Translate3D(gp.x(), gp.y(), gp.z());
0487       }
0488 
0489       const CLHEP::Hep3Vector tt(tr.getTranslation());
0490       trVector.emplace_back(tt.x());
0491       trVector.emplace_back(tt.y());
0492       trVector.emplace_back(tt.z());
0493       if (6 == numberOfTransformParms()) {
0494         const CLHEP::HepRotation rr(tr.getRotation());
0495         const ROOT::Math::Transform3D rtr(
0496             rr.xx(), rr.xy(), rr.xz(), tt.x(), rr.yx(), rr.yy(), rr.yz(), tt.y(), rr.zx(), rr.zy(), rr.zz(), tt.z());
0497         ROOT::Math::EulerAngles ea;
0498         rtr.GetRotation(ea);
0499         trVector.emplace_back(ea.Phi());
0500         trVector.emplace_back(ea.Theta());
0501         trVector.emplace_back(ea.Psi());
0502       }
0503     }
0504   }
0505 }
0506 
0507 DetId HGCalTBGeometry::getGeometryDetId(DetId detId) const {
0508   return static_cast<DetId>(HGCalDetId(detId).geometryCell());
0509 }
0510 
0511 #include "FWCore/Utilities/interface/typelookup.h"
0512 
0513 TYPELOOKUP_DATA_REG(HGCalTBGeometry);