HGCalGeometry/src/FastTimeGeometry.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "Geometry/HGCalGeometry/interface/FastTimeGeometry.h"
0003 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007
0008 #include <cmath>
0009
0010 #include <Math/Transform3D.h>
0011 #include <Math/EulerAngles.h>
0012
0013 typedef CaloCellGeometry::Tr3D Tr3D;
0014 typedef std::vector<float> ParmVec;
0015
0016 //#define EDM_ML_DEBUG
0017
0018 FastTimeGeometry::FastTimeGeometry(const FastTimeTopology& topology_)
0019     : m_topology(topology_),
0020       m_cellVec(topology_.totalGeomModules()),
0021       m_validGeomIds(topology_.totalGeomModules()),
0022       m_Type(topology_.detectorType()),
0023       m_subdet(topology_.subDetector()) {
0024   m_validIds.reserve(topology().totalModules());
0025 #ifdef EDM_ML_DEBUG
0026   edm::LogVerbatim("FastTimeGeom") << "Expected total # of Geometry Modules " << topology().totalGeomModules();
0027 #endif
0028 }
0029
0030 FastTimeGeometry::~FastTimeGeometry() {}
0031
0032 void FastTimeGeometry::fillNamedParams(DDFilteredView fv) {}
0033
0034 void FastTimeGeometry::initializeParms() {}
0035
0036 void FastTimeGeometry::localCorners(Pt3DVec& lc, const CCGFloat* pv, unsigned int i, Pt3D& ref) {
0037   FlatTrd::localCorners(lc, pv, ref);
0038 }
0039
0040 void FastTimeGeometry::newCell(
0041     const GlobalPoint& f1, const GlobalPoint& f2, const GlobalPoint& f3, const CCGFloat* parm, const DetId& detId) {
0042   FastTimeTopology::DecodedDetId id = topology().decode(detId);
0043   DetId geomId = (DetId)(FastTimeDetId(detId).geometryCell());
0044   int nEtaZ = topology().dddConstants().numberEtaZ(m_Type);
0045   int nPhi = topology().dddConstants().numberPhi(m_Type);
0046
0047   const uint32_t cellIndex(topology().detId2denseGeomId(detId));
0048
0049   m_cellVec.at(cellIndex) = FlatTrd(cornersMgr(), f1, f2, f3, parm);
0050   m_validGeomIds.at(cellIndex) = geomId;
0051
0052 #ifdef EDM_ML_DEBUG
0053   unsigned int nOld = m_validIds.size();
0054 #endif
0055   for (int etaZ = 1; etaZ <= nEtaZ; ++etaZ) {
0056     id.iEtaZ = etaZ;
0057     for (int phi = 1; phi <= nPhi; ++phi) {
0058       id.iPhi = phi;
0059       DetId idc = topology().encode(id);
0060       if (topology().valid(idc)) {
0061         m_validIds.emplace_back(idc);
0062       }
0063     }
0064   }
0065
0066 #ifdef EDM_ML_DEBUG
0067   edm::LogVerbatim("FastTimeGeom") << "FastTimeGeometry::newCell-> [" << cellIndex << "] front:" << f1.x() << '/'
0068                                    << f1.y() << '/' << f1.z() << " back:" << f2.x() << '/' << f2.y() << '/' << f2.z()
0069                                    << " eta|phi " << m_cellVec[cellIndex].etaPos() << ":"
0070                                    << m_cellVec[cellIndex].phiPos() << " id:" << FastTimeDetId(detId)
0071                                    << " with valid DetId from " << nOld << " to " << m_validIds.size();
0072   edm::LogVerbatim("FastTimeGeom") << "Cell[" << cellIndex << "] " << std::hex << geomId.rawId() << ":"
0073                                    << m_validGeomIds[cellIndex].rawId() << std::dec;
0074 #endif
0075 }
0076
0077 std::shared_ptr<const CaloCellGeometry> FastTimeGeometry::getGeometry(const DetId& id) const {
0078   if (id == DetId())
0079     return nullptr;  // nothing to get
0080   DetId geoId = (DetId)(FastTimeDetId(id).geometryCell());
0081   const uint32_t cellIndex(topology().detId2denseGeomId(geoId));
0082   const GlobalPoint pos = (id != geoId) ? getPosition(id) : GlobalPoint();
0083   return cellGeomPtr(cellIndex, pos);
0084 }
0085
0086 bool FastTimeGeometry::present(const DetId& id) const {
0087   if (id == DetId())
0088     return false;
0089   DetId geoId = (DetId)(FastTimeDetId(id).geometryCell());
0090   const uint32_t index(topology().detId2denseGeomId(geoId));
0091   return (nullptr != getGeometryRawPtr(index));
0092 }
0093
0094 GlobalPoint FastTimeGeometry::getPosition(const DetId& id) const {
0095   FastTimeDetId id_ = FastTimeDetId(id);
0096   auto pos = topology().dddConstants().getPosition(m_Type, id_.ieta(), id_.iphi(), id_.zside());
0097   return GlobalPoint(0.1 * pos.x(), 0.1 * pos.y(), 0.1 * pos.z());
0098 }
0099
0100 FastTimeGeometry::CornersVec FastTimeGeometry::getCorners(const DetId& id) const {
0101   FastTimeDetId id_ = FastTimeDetId(id);
0102   auto corners = topology().dddConstants().getCorners(m_Type, id_.ieta(), id_.iphi(), id_.zside());
0103   FastTimeGeometry::CornersVec out;
0104   for (const auto& corner : corners) {
0105     out.emplace_back(0.1 * corner.x(), 0.1 * corner.y(), 0.1 * corner.z());
0106   }
0107   return out;
0108 }
0109
0110 DetId FastTimeGeometry::getClosestCell(const GlobalPoint& r) const {
0111   int zside = (r.z() > 0) ? 1 : -1;
0112   std::pair<int, int> etaZPhi;
0113   if (m_Type == 1) {
0114     double zz = (zside > 0) ? r.z() : -r.z();
0115     etaZPhi = topology().dddConstants().getZPhi(zz, r.phi());
0116   } else {
0117     double phi = (zside > 0) ? static_cast<double>(r.phi()) : atan2(r.y(), -r.x());
0118     // Cast needed to resolve compile-time ambiguity of ? operator between
0119     // convertible Phi class and atan2 template function.
0120
0121     etaZPhi = topology().dddConstants().getEtaPhi(r.perp(), phi);
0122   }
0123   FastTimeDetId id = FastTimeDetId(m_Type, etaZPhi.first, etaZPhi.second, zside);
0124 #ifdef EDM_ML_DEBUG
0125   edm::LogVerbatim("FastTimeGeom") << "getClosestCell: for (" << r.x() << ", " << r.y() << ", " << r.z() << ")  Id "
0126                                    << id.type() << ":" << id.zside() << ":" << id.ieta() << ":" << id.iphi();
0127 #endif
0128
0129   return (topology().valid(id) ? DetId(id) : DetId());
0130 }
0131
0132 FastTimeGeometry::DetIdSet FastTimeGeometry::getCells(const GlobalPoint& r, double dR) const {
0133   FastTimeGeometry::DetIdSet dss;
0134   return dss;
0135 }
0136
0137 std::string FastTimeGeometry::cellElement() const {
0138   if (m_Type == 1)
0139     return "FastTimeBarrel";
0140   else if (m_Type == 2)
0141     return "FastTimeEndcap";
0142   else
0143     return "Unknown";
0144 }
0145
0146 unsigned int FastTimeGeometry::indexFor(const DetId& id) const {
0147   unsigned int cellIndex = m_cellVec.size();
0148   if (id != DetId()) {
0149     DetId geoId = (DetId)(FastTimeDetId(id).geometryCell());
0150     cellIndex = topology().detId2denseGeomId(geoId);
0151 #ifdef EDM_ML_DEBUG
0152     edm::LogVerbatim("FastTimeGeom") << "indexFor " << std::hex << id.rawId() << ":" << geoId.rawId() << std::dec
0153                                      << " index " << cellIndex;
0154 #endif
0155   }
0156   return cellIndex;
0157 }
0158
0159 unsigned int FastTimeGeometry::sizeForDenseIndex() const { return topology().totalGeomModules(); }
0160
0161 const CaloCellGeometry* FastTimeGeometry::getGeometryRawPtr(uint32_t index) const {
0162   // Modify the RawPtr class
0163   const CaloCellGeometry* cell(&m_cellVec[index]);
0164   return (m_cellVec.size() < index || nullptr == cell->param() ? nullptr : cell);
0165 }
0166
0167 std::shared_ptr<const CaloCellGeometry> FastTimeGeometry::cellGeomPtr(uint32_t index) const {
0168   if ((index >= m_cellVec.size()) || (m_validGeomIds[index].rawId() == 0))
0169     return nullptr;
0170   static const auto do_not_delete = [](const void*) {};
0171   auto cell = std::shared_ptr<const CaloCellGeometry>(&m_cellVec[index], do_not_delete);
0172   if (nullptr == cell->param())
0173     return nullptr;
0174   return cell;
0175 }
0176
0177 std::shared_ptr<const CaloCellGeometry> FastTimeGeometry::cellGeomPtr(uint32_t index, const GlobalPoint& pos) const {
0178   if ((index >= m_cellVec.size()) || (m_validGeomIds[index].rawId() == 0))
0179     return nullptr;
0180   if (pos == GlobalPoint())
0181     return cellGeomPtr(index);
0182   auto cell = std::make_shared<FlatTrd>(m_cellVec[index]);
0183   cell->setPosition(pos);
0184 #ifdef EDM_ML_DEBUG
0185   edm::LogVerbatim("FastTimeGeomX") << "cellGeomPtr " << pos << ":" << cell;
0186 #endif
0187   if (nullptr == cell->param())
0188     return nullptr;
0189   return cell;
0190 }
0191
0192 void FastTimeGeometry::addValidID(const DetId& id) {
0193   edm::LogError("FastTimeGeom") << "FastTimeGeometry::addValidID is not implemented";
0194 }
0195
0196 // FIXME: Change sorting algorithm if needed
0197 namespace {
0198   struct rawIdSort {
0199     bool operator()(const DetId& a, const DetId& b) { return (a.rawId() < b.rawId()); }
0200   };
0201 }  // namespace
0202
0203 void FastTimeGeometry::sortDetIds(void) {
0204   m_validIds.shrink_to_fit();
0205   std::sort(m_validIds.begin(), m_validIds.end(), rawIdSort());
0206 }
0207
0208 void FastTimeGeometry::getSummary(CaloSubdetectorGeometry::TrVec& trVector,
0209                                   CaloSubdetectorGeometry::IVec& iVector,
0210                                   CaloSubdetectorGeometry::DimVec& dimVector,
0211                                   CaloSubdetectorGeometry::IVec& dinsVector) const {
0212   unsigned int numberOfCells = topology().totalGeomModules();  // total Geom Modules both sides
0213   unsigned int numberOfShapes = FastTimeGeometry::k_NumberOfShapes;
0214   unsigned int numberOfParametersPerShape = FastTimeGeometry::k_NumberOfParametersPerShape;
0215
0216   trVector.reserve(numberOfCells * numberOfTransformParms());
0217   iVector.reserve(numberOfCells);
0218   dimVector.reserve(numberOfShapes * numberOfParametersPerShape);
0219   dinsVector.reserve(numberOfCells);
0220
0221   for (unsigned int k = 0; k < topology().totalGeomModules(); ++k) {
0222     ParmVec params(FastTimeGeometry::k_NumberOfParametersPerShape, 0);
0223     params[0] = topology().dddConstants().getZHalf(m_Type);
0224     params[1] = params[2] = 0;
0225     params[3] = params[7] = topology().dddConstants().getRin(m_Type);
0226     params[4] = params[8] = topology().dddConstants().getRout(m_Type);
0227     params[5] = params[9] = topology().dddConstants().getRout(m_Type);
0228     params[6] = params[10] = 0;
0229     params[11] = (k == 0) ? 1.0 : -1.0;
0230     dimVector.insert(dimVector.end(), params.begin(), params.end());
0231   }
0232
0233   for (unsigned int i(0); i < numberOfCells; ++i) {
0234     DetId detId = m_validGeomIds[i];
0235     dinsVector.emplace_back(topology().detId2denseGeomId(detId));
0236     iVector.emplace_back(1);
0237
0238     Tr3D tr;
0239     auto ptr(cellGeomPtr(i));
0240     if (nullptr != ptr) {
0241       ptr->getTransform(tr, (Pt3DVec*)nullptr);
0242
0243       if (Tr3D() == tr) {  // there is no rotation
0244         const GlobalPoint& gp(ptr->getPosition());
0245         tr = HepGeom::Translate3D(gp.x(), gp.y(), gp.z());
0246       }
0247
0248       const CLHEP::Hep3Vector tt(tr.getTranslation());
0249       trVector.emplace_back(tt.x());
0250       trVector.emplace_back(tt.y());
0251       trVector.emplace_back(tt.z());
0252       if (6 == numberOfTransformParms()) {
0253         const CLHEP::HepRotation rr(tr.getRotation());
0254         const ROOT::Math::Transform3D rtr(
0255             rr.xx(), rr.xy(), rr.xz(), tt.x(), rr.yx(), rr.yy(), rr.yz(), tt.y(), rr.zx(), rr.zy(), rr.zz(), tt.z());
0256         ROOT::Math::EulerAngles ea;
0257         rtr.GetRotation(ea);
0258         trVector.emplace_back(ea.Phi());
0259         trVector.emplace_back(ea.Theta());
0260         trVector.emplace_back(ea.Psi());
0261       }
0262     }
0263   }
0264 }
0265
0266 #include "FWCore/Utilities/interface/typelookup.h"
0267
0268 TYPELOOKUP_DATA_REG(FastTimeGeometry);