CaloGeometry/src/CaloCellGeometry.cc

0001 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0002 #include <CLHEP/Geometry/Plane3D.h>
0003
0004 typedef CaloCellGeometry::CCGFloat CCGFloat;
0005 typedef CaloCellGeometry::Pt3D Pt3D;
0006 typedef CaloCellGeometry::Pt3DVec Pt3DVec;
0007 typedef CaloCellGeometry::Tr3D Tr3D;
0008
0009 typedef HepGeom::Vector3D<CCGFloat> FVec3D;
0010 typedef HepGeom::Plane3D<CCGFloat> Plane3D;
0011
0012 typedef HepGeom::Vector3D<double> DVec3D;
0013 typedef HepGeom::Plane3D<double> DPlane3D;
0014 typedef HepGeom::Point3D<double> DPt3D;
0015 typedef std::vector<DPt3D> DPt3DVec;
0016
0017 const float CaloCellGeometry::k_ScaleFromDDDtoGeant(0.1);
0018
0019 CaloCellGeometry::CaloCellGeometry()
0020     : m_refPoint(0., 0., 0.), m_corners(), m_parms((CCGFloat*)nullptr), m_dEta(0), m_dPhi(0) {}
0021
0022 CaloCellGeometry::~CaloCellGeometry() {}
0023
0024 CaloCellGeometry::CaloCellGeometry(CornersVec::const_reference gp, CornersMgr* mgr, const CCGFloat* par)
0025     : m_refPoint(gp), m_corners(mgr), m_parms(par), m_rep(gp.perp(), gp.eta(), gp.barePhi()), m_dEta(0), m_dPhi(0) {}
0026
0027 CaloCellGeometry::CaloCellGeometry(const CornersVec& cv, const CCGFloat* par)
0028     : m_refPoint(0.25 * (cv[0].x() + cv[1].x() + cv[2].x() + cv[3].x()),
0029                  0.25 * (cv[0].y() + cv[1].y() + cv[2].y() + cv[3].y()),
0030                  0.25 * (cv[0].z() + cv[1].z() + cv[2].z() + cv[3].z())),
0031       m_corners(cv),
0032       m_parms(par),
0033       m_rep(m_refPoint.perp(), m_refPoint.eta(), m_refPoint.barePhi()),
0034       m_dEta(0),
0035       m_dPhi(0) {}
0036
0037 std::ostream& operator<<(std::ostream& s, const CaloCellGeometry& cell) {
0038   s << ", Center: " << cell.getPosition() << std::endl;
0039
0040   if (cell.emptyCorners()) {
0041     s << "Corners vector is empty." << std::endl;
0042   } else {
0043     const CaloCellGeometry::CornersVec& corners(cell.getCorners());
0044     for (unsigned int i(0); i != corners.size(); ++i) {
0045       s << "Corner: " << corners[i] << std::endl;
0046     }
0047   }
0048   return s;
0049 }
0050
0051 void CaloCellGeometry::getTransform(Tr3D& tr, Pt3DVec* lptr) const {
0052   const GlobalPoint& p(CaloCellGeometry::getPosition());
0053   const Pt3D gFront(p.x(), p.y(), p.z());
0054   const DPt3D dgFront(p.x(), p.y(), p.z());
0055
0056   Pt3D lFront;
0057   assert(nullptr != param());
0058   Pt3DVec lc(8, Pt3D(0, 0, 0));
0059   vocalCorners(lc, param(), lFront);
0060
0061   DPt3D dlFront(lFront.x(), lFront.y(), lFront.z());
0062
0063   const Pt3D lBack(0.25 * (lc[4] + lc[5] + lc[6] + lc[7]));
0064   const DPt3D dlBack(lBack.x(), lBack.y(), lBack.z());
0065
0066   const Pt3D dlOne(lc[0].x(), lc[0].y(), lc[0].z());
0067
0068   const CornersVec& cor(getCorners());
0069   DPt3DVec dkor(8, DPt3D(0, 0, 0));
0070   for (unsigned int i(0); i != 8; ++i) {
0071     dkor[i] = DPt3D(cor[i].x(), cor[i].y(), cor[i].z());
0072   }
0073
0074   DPt3D dgBack(0.25 * (dkor[4] + dkor[5] + dkor[6] + dkor[7]));
0075
0076   const DVec3D dgAxis((dgBack - dgFront).unit());
0077
0078   dgBack = (dgFront + (dlBack - dlFront).mag() * dgAxis);
0079   const DPt3D dgOneT(dgFront + (dlOne - dlFront).mag() * (dkor[0] - dgFront).unit());
0080
0081   const double dlangle((dlBack - dlFront).angle(dlOne - dlFront));
0082   const double dgangle((dgBack - dgFront).angle(dgOneT - dgFront));
0083   const double ddangle(dlangle - dgangle);
0084
0085   const DPlane3D dgPl(dgFront, dgOneT, dgBack);
0086   const DPt3D dp2(dgFront + dgPl.normal().unit());
0087
0088   const DPt3D dgOne(dgFront + HepGeom::Rotate3D(-ddangle, dgFront, dp2) * DVec3D(dgOneT - dgFront));
0089
0090   tr = Tr3D(dlFront, dlBack, dlOne, dgFront, dgBack, dgOne);
0091
0092   if (nullptr != lptr)
0093     (*lptr) = lc;
0094 }
0095
0096 const float* CaloCellGeometry::checkParmPtr(const std::vector<float>& vv, CaloCellGeometry::ParVecVec& pvv) {
0097   const float* pP(nullptr);
0098
0099   for (unsigned int ii(0); ii != pvv.size(); ++ii) {
0100     const ParVec& v(pvv[ii]);
0101     assert(v.size() == vv.size());
0102
0103     bool same(true);
0104     for (unsigned int j(0); j != vv.size(); ++j) {
0105       same = same && (fabs(vv[j] - v[j]) < 1.e-6);
0106       if (!same)
0107         break;
0108     }
0109     if (same) {
0110       pP = &(*v.begin());
0111       break;
0112     }
0113   }
0114   return pP;
0115 }
0116
0117 const float* CaloCellGeometry::getParmPtr(const std::vector<float>& vv,
0118                                           CaloCellGeometry::ParMgr* mgr,
0119                                           CaloCellGeometry::ParVecVec& pvv) {
0120   const float* pP(checkParmPtr(vv, pvv));
0121
0122   if (nullptr == pP) {
0123     pvv.emplace_back(ParVec(mgr));
0124     ParVec& back(pvv.back());
0125     for (unsigned int i(0); i != vv.size(); ++i) {
0126       back[i] = vv[i];
0127     }
0128     pP = &(*back.begin());
0129   }
0130   return pP;
0131 }
0132
0133 bool CaloCellGeometry::inside(const GlobalPoint& point) const {
0134   bool ans(false);
0135   const Pt3D p(point.x(), point.y(), point.z());
0136   const CornersVec& cog(getCorners());
0137   Pt3D co[8];
0138   for (unsigned int i(0); i != 8; ++i) {
0139     co[i] = Pt3D(cog[i].x(), cog[i].y(), cog[i].z());
0140   }
0141
0142   const Plane3D AA(co[0], co[1], co[2]);  // z<0
0143   const Plane3D BB(co[6], co[5], co[4]);  // z>0
0144
0145   if (AA.distance(p) * BB.distance(p) >= 0) {
0146     const Plane3D CC(co[0], co[4], co[5]);  // x<0
0147     const Plane3D DD(co[2], co[6], co[7]);  // x>0
0148     if (CC.distance(p) * DD.distance(p) >= 0) {
0149       const Plane3D EE(co[3], co[7], co[4]);  // y<0
0150       const Plane3D FF(co[1], co[5], co[6]);  // y>0
0151       if (EE.distance(p) * FF.distance(p) >= 0) {
0152         ans = true;
0153       }
0154     }
0155   }
0156   return ans;
0157 }