HGCalGeometry.cc

Go to the documentation of this file.

/* for High Granularity Calorimeter
 * This geometry is essentially driven by topology, 
 * which is thus encapsulated in this class. 
 * This makes this geometry not suitable to be loaded
 * by regular CaloGeometryLoader<T>
 */
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"

#include <cmath>

#include <Math/Transform3D.h>
#include <Math/EulerAngles.h>

typedef CaloCellGeometry::Tr3D Tr3D;
typedef std::vector<float> ParmVec;

//#define EDM_ML_DEBUG

HGCalGeometry::HGCalGeometry(const HGCalTopology& topology_)
    : m_topology(topology_),
      m_validGeomIds(topology_.totalGeomModules()),
      mode_(topology_.geomMode()),
      m_det(topology_.detector()),
      m_subdet(topology_.subDetector()),
      twoBysqrt3_(2.0 / std::sqrt(3.0)) {
  if (m_det == DetId::HGCalHSc) {
    m_cellVec2 = CellVec2(topology_.totalGeomModules());
  } else {
    m_cellVec = CellVec(topology_.totalGeomModules());
  }
  m_validIds.reserve(m_topology.totalModules());
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "Expected total # of Geometry Modules " << m_topology.totalGeomModules();
#endif
}

HGCalGeometry::~HGCalGeometry() {}

void HGCalGeometry::fillNamedParams(DDFilteredView fv) {}

void HGCalGeometry::initializeParms() {}

void HGCalGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int i, Pt3D& ref) {
  if (m_det == DetId::HGCalHSc) {
    FlatTrd::localCorners(lc, pv, ref);
  } else {
    FlatHexagon::localCorners(lc, pv, ref);
  }
}

void HGCalGeometry::newCell(
    const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
  DetId geomId = getGeometryDetId(detId);
  int cells(0);
  HGCalTopology::DecodedDetId id = m_topology.decode(detId);
  if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
    cells = m_topology.dddConstants().numberCellsHexagon(id.iSec1);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "NewCell " << HGCalDetId(detId) << " GEOM " << HGCalDetId(geomId);
#endif
  } else if (mode_ == HGCalGeometryMode::Trapezoid) {
    cells = 1;
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "NewCell " << HGCScintillatorDetId(detId) << " GEOM "
                                  << HGCScintillatorDetId(geomId);
#endif
  } else {
    cells = m_topology.dddConstants().numberCellsHexagon(id.iLay, id.iSec1, id.iSec2, false);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "NewCell " << HGCSiliconDetId(detId) << " GEOM " << HGCSiliconDetId(geomId);
#endif
  }
  const uint32_t cellIndex(m_topology.detId2denseGeomId(geomId));

  if (m_det == DetId::HGCalHSc) {
    m_cellVec2.at(cellIndex) = FlatTrd(cornersMgr(), f1, f2, f3, parm);
  } else {
    m_cellVec.at(cellIndex) = FlatHexagon(cornersMgr(), f1, f2, f3, parm);
  }
  m_validGeomIds.at(cellIndex) = geomId;

#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "Store for DetId " << std::hex << detId.rawId() << " GeomId " << geomId.rawId()
                                << std::dec << " Index " << cellIndex << " cells " << cells;
  unsigned int nOld = m_validIds.size();
#endif
  if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
    for (int cell = 0; cell < cells; ++cell) {
      id.iCell1 = cell;
      DetId idc = m_topology.encode(id);
      if (m_topology.valid(idc)) {
        m_validIds.emplace_back(idc);
#ifdef EDM_ML_DEBUG
        edm::LogVerbatim("HGCalGeom") << "Valid Id [" << cell << "] " << HGCalDetId(idc);
#endif
      }
    }
  } else if (mode_ == HGCalGeometryMode::Trapezoid) {
    DetId idc = m_topology.encode(id);
    if (m_topology.valid(idc)) {
      m_validIds.emplace_back(idc);
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HGCalGeom") << "Valid Id [0] " << HGCScintillatorDetId(idc);
#endif
    } else {
      edm::LogWarning("HGCalGeom") << "Check " << HGCScintillatorDetId(idc) << " from " << HGCScintillatorDetId(detId)
                                   << " ERROR ???";
    }
  } else {
#ifdef EDM_ML_DEBUG
    unsigned int cellAll(0), cellSelect(0);
#endif
    for (int u = 0; u < 2 * cells; ++u) {
      for (int v = 0; v < 2 * cells; ++v) {
        if (((v - u) < cells) && (u - v) <= cells) {
          id.iCell1 = u;
          id.iCell2 = v;
          DetId idc = m_topology.encode(id);
#ifdef EDM_ML_DEBUG
          ++cellAll;
#endif
          if (m_topology.dddConstants().cellInLayer(id.iSec1, id.iSec2, u, v, id.iLay, true)) {
            m_validIds.emplace_back(idc);
#ifdef EDM_ML_DEBUG
            ++cellSelect;
            edm::LogVerbatim("HGCalGeom") << "Valid Id [" << u << ", " << v << "] " << HGCSiliconDetId(idc);
#endif
          }
        }
      }
    }
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "HGCalGeometry keeps " << cellSelect << " out of " << cellAll << " for wafer "
                                  << id.iSec1 << ":" << id.iSec2 << " in "
                                  << " layer " << id.iLay;
#endif
  }
#ifdef EDM_ML_DEBUG
  if (m_det == DetId::HGCalHSc) {
    edm::LogVerbatim("HGCalGeom") << "HGCalGeometry::newCell-> [" << cellIndex << "]"
                                  << " front:" << f1.x() << '/' << f1.y() << '/' << f1.z() << " back:" << f2.x() << '/'
                                  << f2.y() << '/' << f2.z() << " eta|phi " << m_cellVec2[cellIndex].etaPos() << ":"
                                  << m_cellVec2[cellIndex].phiPos();
  } else {
    edm::LogVerbatim("HGCalGeom") << "HGCalGeometry::newCell-> [" << cellIndex << "]"
                                  << " front:" << f1.x() << '/' << f1.y() << '/' << f1.z() << " back:" << f2.x() << '/'
                                  << f2.y() << '/' << f2.z() << " eta|phi " << m_cellVec[cellIndex].etaPos() << ":"
                                  << m_cellVec[cellIndex].phiPos();
  }
  unsigned int nNew = m_validIds.size();
  if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
    edm::LogVerbatim("HGCalGeom") << "ID: " << HGCalDetId(detId) << " with valid DetId from " << nOld << " to " << nNew;
  } else if (mode_ == HGCalGeometryMode::Trapezoid) {
    edm::LogVerbatim("HGCalGeom") << "ID: " << HGCScintillatorDetId(detId) << " with valid DetId from " << nOld
                                  << " to " << nNew;
  } else if (m_topology.isHFNose()) {
    edm::LogVerbatim("HGCalGeom") << "ID: " << HFNoseDetId(detId) << " with valid DetId from " << nOld << " to "
                                  << nNew;
  } else {
    edm::LogVerbatim("HGCalGeom") << "ID: " << HGCSiliconDetId(detId) << " with valid DetId from " << nOld << " to "
                                  << nNew;
  }
  edm::LogVerbatim("HGCalGeom") << "Cell[" << cellIndex << "] " << std::hex << geomId.rawId() << ":"
                                << m_validGeomIds[cellIndex].rawId() << std::dec;
#endif
}

std::shared_ptr<const CaloCellGeometry> HGCalGeometry::getGeometry(const DetId& detId) const {
  if (detId == DetId())
    return nullptr;  // nothing to get
  DetId geomId = getGeometryDetId(detId);
  const uint32_t cellIndex(m_topology.detId2denseGeomId(geomId));
  const GlobalPoint pos = (detId != geomId) ? getPosition(detId) : GlobalPoint();
  return cellGeomPtr(cellIndex, pos);
}

bool HGCalGeometry::present(const DetId& detId) const {
  if (detId == DetId())
    return false;
  DetId geomId = getGeometryDetId(detId);
  const uint32_t index(m_topology.detId2denseGeomId(geomId));
  return (nullptr != getGeometryRawPtr(index));
}

GlobalPoint HGCalGeometry::getPosition(const DetId& detid) const {
  unsigned int cellIndex = indexFor(detid);
  GlobalPoint glob;
  unsigned int maxSize = ((mode_ == HGCalGeometryMode::Trapezoid) ? m_cellVec2.size() : m_cellVec.size());
  if (cellIndex < maxSize) {
    HGCalTopology::DecodedDetId id = m_topology.decode(detid);
    std::pair<float, float> xy;
    if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
      xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
      const HepGeom::Point3D<float> lcoord(xy.first, xy.second, 0);
      glob = m_cellVec[cellIndex].getPosition(lcoord);
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HGCalGeom") << "getPosition:: index " << cellIndex << " Local " << lcoord.x() << ":"
                                    << lcoord.y() << " ID " << id.iCell1 << ":" << id.iSec1 << " Global " << glob;
#endif
    } else if (mode_ == HGCalGeometryMode::Trapezoid) {
      const HepGeom::Point3D<float> lcoord(0, 0, 0);
      glob = m_cellVec2[cellIndex].getPosition(lcoord);
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HGCalGeom") << "getPositionTrap:: index " << cellIndex << " Local " << lcoord.x() << ":"
                                    << lcoord.y() << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iCell1
                                    << " Global " << glob;
#endif
    } else {
      xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false, true);
      const HepGeom::Point3D<float> lcoord(xy.first, xy.second, 0);
      glob = m_cellVec[cellIndex].getPosition(lcoord);
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HGCalGeom") << "getPositionWafer:: index " << cellIndex << " Local " << lcoord.x() << ":"
                                    << lcoord.y() << " ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
                                    << id.iCell1 << ":" << id.iCell2 << " Global " << glob;
#endif
    }
  }
  return glob;
}

double HGCalGeometry::getArea(const DetId& detid) const {
  HGCalGeometry::CornersVec corners = getNewCorners(detid);
  double area(0);
  if (corners.size() > 1) {
    int n = corners.size() - 1;
    int j = n - 1;
    for (int i = 0; i < n; ++i) {
      area += ((corners[j].x() + corners[i].x()) * (corners[i].y() - corners[j].y()));
      j = i;
    }
  }
  return (0.5 * area);
}

HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& detid) const {
  unsigned int ncorner = ((m_det == DetId::HGCalHSc) ? FlatTrd::ncorner_ : FlatHexagon::ncorner_);
  HGCalGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
  unsigned int cellIndex = indexFor(detid);
  HGCalTopology::DecodedDetId id = m_topology.decode(detid);
  if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
    GlobalPoint v = getPosition(detid);
    std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(id.iType, id.iSec1);
    float dr = k_half * (rr.second - rr.first);
    float dfi = m_cellVec2[cellIndex].param()[FlatTrd::k_Cell];
    float dz = id.zSide * m_cellVec2[cellIndex].param()[FlatTrd::k_dZ];
    float r = v.perp();
    float fi = v.phi();
    static const int signr[] = {1, 1, -1, -1, 1, 1, -1, -1};
    static const int signf[] = {-1, 1, 1, -1, -1, 1, 1, -1};
    static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
    for (unsigned int i = 0; i < ncorner; ++i) {
      co[i] = GlobalPoint((r + signr[i] * dr) * cos(fi + signf[i] * dfi),
                          (r + signr[i] * dr) * sin(fi + signf[i] * dfi),
                          (v.z() + signz[i] * dz));
    }
  } else if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
    std::pair<float, float> xy;
    if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
      xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
      float dx = m_cellVec[cellIndex].param()[FlatHexagon::k_r];
      float dy = k_half * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
      float dz = m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
      static const int signx[] = {0, -1, -1, 0, 1, 1, 0, -1, -1, 0, 1, 1};
      static const int signy[] = {-2, -1, 1, 2, 1, -1, -2, -1, 1, 2, 1, -1};
      static const int signz[] = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
      for (unsigned int i = 0; i < ncorner; ++i) {
        const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, signz[i] * dz);
        co[i] = m_cellVec[cellIndex].getPosition(lcoord);
      }
    } else {
      xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false);
      float dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
      float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
      float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
      static const int signx[] = {1, -1, -2, -1, 1, 2, 1, -1, -2, -1, 1, 2};
      static const int signy[] = {1, 1, 0, -1, -1, 0, 1, 1, 0, -1, -1, 0};
      static const int signz[] = {-1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1};
      for (unsigned int i = 0; i < ncorner; ++i) {
        const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, signz[i] * dz);
        co[i] = m_cellVec[cellIndex].getPosition(lcoord);
      }
    }
  }
  return co;
}

HGCalGeometry::CornersVec HGCalGeometry::get8Corners(const DetId& detid) const {
  unsigned int ncorner = FlatTrd::ncorner_;
  HGCalGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
  unsigned int cellIndex = indexFor(detid);
  HGCalTopology::DecodedDetId id = m_topology.decode(detid);
  if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
    GlobalPoint v = getPosition(detid);
    std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(id.iType, id.iSec1);
    float dr = k_half * (rr.second - rr.first);
    float dfi = m_cellVec2[cellIndex].param()[FlatTrd::k_Cell];
    float dz = id.zSide * m_cellVec2[cellIndex].param()[FlatTrd::k_dZ];
    float r = v.perp();
    float fi = v.phi();
    static const int signr[] = {1, 1, -1, -1, 1, 1, -1, -1};
    static const int signf[] = {-1, 1, 1, -1, -1, 1, 1, -1};
    static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
    for (unsigned int i = 0; i < ncorner; ++i) {
      co[i] = GlobalPoint((r + signr[i] * dr) * cos(fi + signf[i] * dfi),
                          (r + signr[i] * dr) * sin(fi + signf[i] * dfi),
                          (v.z() + signz[i] * dz));
    }
  } else if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
    std::pair<float, float> xy;
    float dx(0);
    if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
      xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
      dx = m_cellVec[cellIndex].param()[FlatHexagon::k_r];
    } else {
      xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false);
      dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
    }
    static const int signx[] = {-1, -1, 1, 1, -1, -1, 1, 1};
    static const int signy[] = {-1, 1, 1, -1, -1, 1, 1, -1};
    static const int signz[] = {-1, -1, -1, -1, 1, 1, 1, 1};
    float dz = m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
    for (unsigned int i = 0; i < ncorner; ++i) {
      const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dx, signz[i] * dz);
      co[i] = m_cellVec[cellIndex].getPosition(lcoord);
    }
  }
  return co;
}

HGCalGeometry::CornersVec HGCalGeometry::getNewCorners(const DetId& detid) const {
  unsigned int ncorner = (m_det == DetId::HGCalHSc) ? 5 : 7;
  HGCalGeometry::CornersVec co(ncorner, GlobalPoint(0, 0, 0));
  unsigned int cellIndex = indexFor(detid);
  HGCalTopology::DecodedDetId id = m_topology.decode(detid);
  if (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc) {
    GlobalPoint v = getPosition(detid);
    std::pair<double, double> rr = m_topology.dddConstants().cellSizeTrap(id.iType, id.iSec1);
    float dr = k_half * (rr.second - rr.first);
    float dfi = m_cellVec2[cellIndex].param()[FlatTrd::k_Cell];
    float dz = -id.zSide * m_cellVec2[cellIndex].param()[FlatTrd::k_dZ];
    float r = v.perp();
    float fi = v.phi();
    static const int signr[] = {1, 1, -1, -1};
    static const int signf[] = {-1, 1, 1, -1};
    for (unsigned int i = 0; i < ncorner - 1; ++i) {
      co[i] = GlobalPoint(
          (r + signr[i] * dr) * cos(fi + signf[i] * dfi), (r + signr[i] * dr) * sin(fi + signf[i] * dfi), (v.z() + dz));
    }
    co[ncorner - 1] = GlobalPoint(0, 0, -2 * dz);
  } else if (cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) {
    std::pair<float, float> xy;
    if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
      xy = m_topology.dddConstants().locateCellHex(id.iCell1, id.iSec1, true);
    } else {
      xy = m_topology.dddConstants().locateCell(id.iLay, id.iSec1, id.iSec2, id.iCell1, id.iCell2, true, false);
    }
    float dx = k_fac2 * m_cellVec[cellIndex].param()[FlatHexagon::k_r];
    float dy = k_fac1 * m_cellVec[cellIndex].param()[FlatHexagon::k_R];
    float dz = -id.zSide * m_cellVec[cellIndex].param()[FlatHexagon::k_dZ];
    static const int signx[] = {1, -1, -2, -1, 1, 2};
    static const int signy[] = {1, 1, 0, -1, -1, 0};
    for (unsigned int i = 0; i < ncorner - 1; ++i) {
      const HepGeom::Point3D<float> lcoord(xy.first + signx[i] * dx, xy.second + signy[i] * dy, dz);
      co[i] = m_cellVec[cellIndex].getPosition(lcoord);
    }
    co[ncorner - 1] = GlobalPoint(0, 0, -2 * dz);
  }
  return co;
}

DetId HGCalGeometry::neighborZ(const DetId& idin, const GlobalVector& momentum) const {
  DetId idnew;
  HGCalTopology::DecodedDetId id = m_topology.decode(idin);
  int lay = ((momentum.z() * id.zSide > 0) ? (id.iLay + 1) : (id.iLay - 1));
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "neighborz1:: ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
                                << id.iCell1 << ":" << id.iCell2 << " New Layer " << lay << " Range "
                                << m_topology.dddConstants().firstLayer() << ":"
                                << m_topology.dddConstants().lastLayer(true) << " pz " << momentum.z();
#endif
  if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
      (momentum.z() != 0.0)) {
    GlobalPoint v = getPosition(idin);
    double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
    double grad = (z - v.z()) / momentum.z();
    GlobalPoint p(v.x() + grad * momentum.x(), v.y() + grad * momentum.y(), z);
    double r = p.perp();
    auto rlimit = topology().dddConstants().rangeR(z, true);
    if (r >= rlimit.first && r <= rlimit.second)
      idnew = getClosestCell(p);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "neighborz1:: Position " << v << " New Z " << z << ":" << grad << " new position "
                                  << p << " r-limit " << rlimit.first << ":" << rlimit.second;
#endif
  }
  return idnew;
}

DetId HGCalGeometry::neighborZ(const DetId& idin,
                               const MagneticField* bField,
                               int charge,
                               const GlobalVector& momentum) const {
  DetId idnew;
  HGCalTopology::DecodedDetId id = m_topology.decode(idin);
  int lay = ((momentum.z() * id.zSide > 0) ? (id.iLay + 1) : (id.iLay - 1));
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "neighborz2:: ID " << id.iLay << ":" << id.iSec1 << ":" << id.iSec2 << ":"
                                << id.iCell1 << ":" << id.iCell2 << " New Layer " << lay << " Range "
                                << m_topology.dddConstants().firstLayer() << ":"
                                << m_topology.dddConstants().lastLayer(true) << " pz " << momentum.z();
#endif
  if ((lay >= m_topology.dddConstants().firstLayer()) && (lay <= m_topology.dddConstants().lastLayer(true)) &&
      (momentum.z() != 0.0)) {
    GlobalPoint v = getPosition(idin);
    double z = id.zSide * m_topology.dddConstants().waferZ(lay, true);
    FreeTrajectoryState fts(v, momentum, charge, bField);
    Plane::PlanePointer nPlane = Plane::build(Plane::PositionType(0, 0, z), Plane::RotationType());
    AnalyticalPropagator myAP(bField, alongMomentum, 2 * M_PI);
    TrajectoryStateOnSurface tsos = myAP.propagate(fts, *nPlane);
    GlobalPoint p;
    auto rlimit = topology().dddConstants().rangeR(z, true);
    if (tsos.isValid()) {
      p = tsos.globalPosition();
      double r = p.perp();
      if (r >= rlimit.first && r <= rlimit.second)
        idnew = getClosestCell(p);
    }
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "neighborz2:: Position " << v << " New Z " << z << ":" << charge << ":"
                                  << tsos.isValid() << " new position " << p << " r limits " << rlimit.first << ":"
                                  << rlimit.second;
#endif
  }
  return idnew;
}

DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const {
  unsigned int cellIndex = getClosestCellIndex(r);
  if ((cellIndex < m_cellVec.size() && m_det != DetId::HGCalHSc) ||
      (cellIndex < m_cellVec2.size() && m_det == DetId::HGCalHSc)) {
    HGCalTopology::DecodedDetId id = m_topology.decode(m_validGeomIds[cellIndex]);
    HepGeom::Point3D<float> local;
    if (r.z() > 0) {
      local = HepGeom::Point3D<float>(r.x(), r.y(), 0);
      id.zSide = 1;
    } else {
      local = HepGeom::Point3D<float>(-r.x(), r.y(), 0);
      id.zSide = -1;
    }
    if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
      const auto& kxy = m_topology.dddConstants().assignCell(local.x(), local.y(), id.iLay, id.iType, true);
      id.iCell1 = kxy.second;
      id.iSec1 = kxy.first;
      id.iType = m_topology.dddConstants().waferTypeT(kxy.first);
      if (id.iType != 1)
        id.iType = -1;
    } else if (mode_ == HGCalGeometryMode::Trapezoid) {
      id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
      const auto& kxy = m_topology.dddConstants().assignCellTrap(r.x(), r.y(), r.z(), id.iLay, true);
      id.iSec1 = kxy[0];
      id.iCell1 = kxy[1];
      id.iType = kxy[2];
    } else {
      id.iLay = m_topology.dddConstants().getLayer(r.z(), true);
      const auto& kxy = m_topology.dddConstants().assignCellHex(local.x(), local.y(), id.iLay, true);
      id.iSec1 = kxy[0];
      id.iSec2 = kxy[1];
      id.iType = kxy[2];
      id.iCell1 = kxy[3];
      id.iCell2 = kxy[4];
    }
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "getClosestCell: local " << local << " Id " << id.zSide << ":" << id.iLay << ":"
                                  << id.iSec1 << ":" << id.iSec2 << ":" << id.iType << ":" << id.iCell1 << ":"
                                  << id.iCell2;
#endif

    //check if returned cell is valid
    if (id.iCell1 >= 0)
      return m_topology.encode(id);
  }

  //if not valid or out of bounds return a null DetId
  return DetId();
}

HGCalGeometry::DetIdSet HGCalGeometry::getCells(const GlobalPoint& r, double dR) const {
  HGCalGeometry::DetIdSet dss;
  return dss;
}

std::string HGCalGeometry::cellElement() const {
  if (m_subdet == HGCEE || m_det == DetId::HGCalEE)
    return "HGCalEE";
  else if (m_subdet == HGCHEF || m_det == DetId::HGCalHSi)
    return "HGCalHEFront";
  else if (m_subdet == HGCHEB || m_det == DetId::HGCalHSc)
    return "HGCalHEBack";
  else
    return "Unknown";
}

unsigned int HGCalGeometry::indexFor(const DetId& detId) const {
  unsigned int cellIndex = ((m_det == DetId::HGCalHSc) ? m_cellVec2.size() : m_cellVec.size());
  if (detId != DetId()) {
    DetId geomId = getGeometryDetId(detId);
    cellIndex = m_topology.detId2denseGeomId(geomId);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "indexFor " << std::hex << detId.rawId() << ":" << geomId.rawId() << std::dec
                                  << " index " << cellIndex;
#endif
  }
  return cellIndex;
}

unsigned int HGCalGeometry::sizeForDenseIndex() const { return m_topology.totalGeomModules(); }

const CaloCellGeometry* HGCalGeometry::getGeometryRawPtr(uint32_t index) const {
  // Modify the RawPtr class
  if (m_det == DetId::HGCalHSc) {
    if (m_cellVec2.size() < index)
      return nullptr;
    const CaloCellGeometry* cell(&m_cellVec2[index]);
    return (nullptr == cell->param() ? nullptr : cell);
  } else {
    if (m_cellVec2.size() < index)
      return nullptr;
    const CaloCellGeometry* cell(&m_cellVec[index]);
    return (nullptr == cell->param() ? nullptr : cell);
  }
}

std::shared_ptr<const CaloCellGeometry> HGCalGeometry::cellGeomPtr(uint32_t index) const {
  if ((index >= m_cellVec.size() && m_det != DetId::HGCalHSc) ||
      (index >= m_cellVec2.size() && m_det == DetId::HGCalHSc) || (m_validGeomIds[index].rawId() == 0))
    return nullptr;
  static const auto do_not_delete = [](const void*) {};
  if (m_det == DetId::HGCalHSc) {
    auto cell = std::shared_ptr<const CaloCellGeometry>(&m_cellVec2[index], do_not_delete);
    if (nullptr == cell->param())
      return nullptr;
    return cell;
  } else {
    auto cell = std::shared_ptr<const CaloCellGeometry>(&m_cellVec[index], do_not_delete);
    if (nullptr == cell->param())
      return nullptr;
    return cell;
  }
}

std::shared_ptr<const CaloCellGeometry> HGCalGeometry::cellGeomPtr(uint32_t index, const GlobalPoint& pos) const {
  if ((index >= m_cellVec.size() && m_det != DetId::HGCalHSc) ||
      (index >= m_cellVec2.size() && m_det == DetId::HGCalHSc) || (m_validGeomIds[index].rawId() == 0))
    return nullptr;
  if (pos == GlobalPoint())
    return cellGeomPtr(index);
  if (m_det == DetId::HGCalHSc) {
    auto cell = std::make_shared<FlatTrd>(m_cellVec2[index]);
    cell->setPosition(pos);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "cellGeomPtr " << index << ":" << cell;
#endif
    if (nullptr == cell->param())
      return nullptr;
    return cell;
  } else {
    auto cell = std::make_shared<FlatHexagon>(m_cellVec[index]);
    cell->setPosition(pos);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HGCalGeom") << "cellGeomPtr " << index << ":" << cell;
#endif
    if (nullptr == cell->param())
      return nullptr;
    return cell;
  }
}

void HGCalGeometry::addValidID(const DetId& id) {
  edm::LogError("HGCalGeom") << "HGCalGeometry::addValidID is not implemented";
}

unsigned int HGCalGeometry::getClosestCellIndex(const GlobalPoint& r) const {
  return ((m_det == DetId::HGCalHSc) ? getClosestCellIndex(r, m_cellVec2) : getClosestCellIndex(r, m_cellVec));
}

template <class T>
unsigned int HGCalGeometry::getClosestCellIndex(const GlobalPoint& r, const std::vector<T>& vec) const {
  float phip = r.phi();
  float zp = r.z();
  float dzmin(9999), dphimin(9999), dphi10(0.175);
  unsigned int cellIndex = vec.size();
  for (unsigned int k = 0; k < vec.size(); ++k) {
    float dphi = phip - vec[k].phiPos();
    while (dphi > M_PI)
      dphi -= 2 * M_PI;
    while (dphi <= -M_PI)
      dphi += 2 * M_PI;
    if (std::abs(dphi) < dphi10) {
      float dz = std::abs(zp - vec[k].getPosition().z());
      if (dz < (dzmin + 0.001)) {
        dzmin = dz;
        if (std::abs(dphi) < (dphimin + 0.01)) {
          cellIndex = k;
          dphimin = std::abs(dphi);
        } else {
          if (cellIndex >= vec.size())
            cellIndex = k;
        }
      }
    }
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "getClosestCellIndex::Input " << zp << ":" << phip << " Index " << cellIndex;
  if (cellIndex < vec.size())
    edm::LogVerbatim("HGCalGeom") << " Cell z " << vec[cellIndex].getPosition().z() << ":" << dzmin << " phi "
                                  << vec[cellIndex].phiPos() << ":" << dphimin;
#endif
  return cellIndex;
}

// FIXME: Change sorting algorithm if needed
namespace {
  struct rawIdSort {
    bool operator()(const DetId& a, const DetId& b) { return (a.rawId() < b.rawId()); }
  };
}  // namespace

void HGCalGeometry::sortDetIds(void) {
  m_validIds.shrink_to_fit();
  std::sort(m_validIds.begin(), m_validIds.end(), rawIdSort());
}

void HGCalGeometry::getSummary(CaloSubdetectorGeometry::TrVec& trVector,
                               CaloSubdetectorGeometry::IVec& iVector,
                               CaloSubdetectorGeometry::DimVec& dimVector,
                               CaloSubdetectorGeometry::IVec& dinsVector) const {
  unsigned int numberOfCells = m_topology.totalGeomModules();  // total Geom Modules both sides
  unsigned int numberOfShapes = k_NumberOfShapes;
  unsigned int numberOfParametersPerShape = ((m_det == DetId::HGCalHSc) ? (unsigned int)(k_NumberOfParametersPerTrd)
                                                                        : (unsigned int)(k_NumberOfParametersPerHex));

  trVector.reserve(numberOfCells * numberOfTransformParms());
  iVector.reserve(numberOfCells);
  dimVector.reserve(numberOfShapes * numberOfParametersPerShape);
  dinsVector.reserve(numberOfCells);

  for (unsigned itr = 0; itr < m_topology.dddConstants().getTrFormN(); ++itr) {
    HGCalParameters::hgtrform mytr = m_topology.dddConstants().getTrForm(itr);
    int layer = mytr.lay;

    if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
      for (int wafer = 0; wafer < m_topology.dddConstants().sectors(); ++wafer) {
        if (m_topology.dddConstants().waferInLayer(wafer, layer, true)) {
          HGCalParameters::hgtrap vol = m_topology.dddConstants().getModule(wafer, true, true);
          ParmVec params(numberOfParametersPerShape, 0);
          params[FlatHexagon::k_dZ] = vol.dz;
          params[FlatHexagon::k_r] = vol.cellSize;
          params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
          dimVector.insert(dimVector.end(), params.begin(), params.end());
        }
      }
    } else if (mode_ == HGCalGeometryMode::Trapezoid) {
      int indx = m_topology.dddConstants().layerIndex(layer, true);
      for (int md = m_topology.dddConstants().getParameter()->firstModule_[indx];
           md <= m_topology.dddConstants().getParameter()->lastModule_[indx];
           ++md) {
        HGCalParameters::hgtrap vol = m_topology.dddConstants().getModule(md, true, true);
        ParmVec params(numberOfParametersPerShape, 0);
        params[FlatTrd::k_dZ] = vol.dz;
        params[FlatTrd::k_Theta] = params[FlatTrd::k_Phi] = 0;
        params[FlatTrd::k_dY1] = params[FlatTrd::k_dY2] = vol.h;
        params[FlatTrd::k_dX1] = params[FlatTrd::k_dX3] = vol.bl;
        params[FlatTrd::k_dX2] = params[FlatTrd::k_dX4] = vol.tl;
        params[FlatTrd::k_Alp1] = params[FlatTrd::k_Alp2] = vol.alpha;
        params[FlatTrd::k_Cell] = vol.cellSize;
        dimVector.insert(dimVector.end(), params.begin(), params.end());
      }
    } else {
      for (int wafer = 0; wafer < m_topology.dddConstants().sectors(); ++wafer) {
        if (m_topology.dddConstants().waferInLayer(wafer, layer, true)) {
          HGCalParameters::hgtrap vol = m_topology.dddConstants().getModule(wafer, true, true);
          ParmVec params(numberOfParametersPerShape, 0);
          params[FlatHexagon::k_dZ] = vol.dz;
          params[FlatHexagon::k_r] = vol.cellSize;
          params[FlatHexagon::k_R] = twoBysqrt3_ * params[FlatHexagon::k_r];
          dimVector.insert(dimVector.end(), params.begin(), params.end());
        }
      }
    }
  }

  for (unsigned int i(0); i < numberOfCells; ++i) {
    DetId detId = m_validGeomIds[i];
    int layer(0);
    if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
      layer = HGCalDetId(detId).layer();
    } else if (mode_ == HGCalGeometryMode::Trapezoid) {
      layer = HGCScintillatorDetId(detId).layer();
    } else if (m_topology.isHFNose()) {
      layer = HFNoseDetId(detId).layer();
    } else {
      layer = HGCSiliconDetId(detId).layer();
    }
    dinsVector.emplace_back(m_topology.detId2denseGeomId(detId));
    iVector.emplace_back(layer);

    Tr3D tr;
    auto ptr = cellGeomPtr(i);
    if (nullptr != ptr) {
      ptr->getTransform(tr, (Pt3DVec*)nullptr);

      if (Tr3D() == tr) {  // there is no rotation
        const GlobalPoint& gp(ptr->getPosition());
        tr = HepGeom::Translate3D(gp.x(), gp.y(), gp.z());
      }

      const CLHEP::Hep3Vector tt(tr.getTranslation());
      trVector.emplace_back(tt.x());
      trVector.emplace_back(tt.y());
      trVector.emplace_back(tt.z());
      if (6 == numberOfTransformParms()) {
        const CLHEP::HepRotation rr(tr.getRotation());
        const ROOT::Math::Transform3D rtr(
            rr.xx(), rr.xy(), rr.xz(), tt.x(), rr.yx(), rr.yy(), rr.yz(), tt.y(), rr.zx(), rr.zy(), rr.zz(), tt.z());
        ROOT::Math::EulerAngles ea;
        rtr.GetRotation(ea);
        trVector.emplace_back(ea.Phi());
        trVector.emplace_back(ea.Theta());
        trVector.emplace_back(ea.Psi());
      }
    }
  }
}

DetId HGCalGeometry::getGeometryDetId(DetId detId) const {
  DetId geomId;
  if ((mode_ == HGCalGeometryMode::Hexagon) || (mode_ == HGCalGeometryMode::HexagonFull)) {
    geomId = static_cast<DetId>(HGCalDetId(detId).geometryCell());
  } else if (mode_ == HGCalGeometryMode::Trapezoid) {
    geomId = static_cast<DetId>(HGCScintillatorDetId(detId).geometryCell());
  } else if (m_topology.isHFNose()) {
    geomId = static_cast<DetId>(HFNoseDetId(detId).geometryCell());
  } else {
    geomId = static_cast<DetId>(HGCSiliconDetId(detId).geometryCell());
  }
  return geomId;
}

#include "FWCore/Utilities/interface/typelookup.h"

TYPELOOKUP_DATA_REG(HGCalGeometry);