File indexing completed on 2023-03-17 13:02:29
0001 #include "Geometry/CaloGeometry/interface/Line3D.h"
0002 #include "Geometry/CaloGeometry/interface/CaloCellCrossing.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 CaloCellCrossing::CaloCellCrossing(const GlobalPoint& gp,
0006 const GlobalVector& gv,
0007 const CaloCellCrossing::DetIds* di,
0008 const CaloSubdetectorGeometry* sg,
0009 DetId::Detector det,
0010 int subdet,
0011 double small,
0012 bool onewayonly)
0013 : m_gp(gp), m_gv(gv) {
0014 const double eps(fabs(small) > 1.e-15 ? fabs(small) : 1.e-10);
0015 static unsigned int k_rlen(30);
0016 m_detId.reserve(k_rlen);
0017 m_ctr.reserve(k_rlen);
0018 m_entr.reserve(k_rlen);
0019 m_exit.reserve(k_rlen);
0020
0021 const HepLine3D line(
0022 HepGeom::Point3D<double>(gp.x(), gp.y(), gp.z()), HepGeom::Vector3D<double>(gv.x(), gv.y(), gv.z()), eps);
0023
0024 LogDebug("CaloCellCrossing") << "*** Line: pt=" << line.pt() << ", unitvec=" << line.uv();
0025
0026 const DetIds& ids(nullptr == di ? sg->getValidDetIds(det, subdet) : *di);
0027
0028 for (auto dId : ids) {
0029 unsigned int found(0);
0030 auto cg = sg->getGeometry(dId);
0031 const CaloCellGeometry::CornersVec& gc(cg->getCorners());
0032 const HepGeom::Point3D<double> fr(cg->getPosition().x(), cg->getPosition().y(), cg->getPosition().z());
0033 const double bCut2((gc[0] - gc[6]).mag2());
0034
0035 if ((!onewayonly || eps < HepGeom::Vector3D<double>(fr - line.pt()).dot(line.uv())) &&
0036 bCut2 > line.dist2(fr))
0037 {
0038 LogDebug("CaloCellCrossing") << "*** fr=" << fr << ", bCut =" << sqrt(bCut2) << ", dis=" << line.dist(fr);
0039 const HepGeom::Point3D<double> cv[8] = {HepGeom::Point3D<double>(gc[0].x(), gc[0].y(), gc[0].z()),
0040 HepGeom::Point3D<double>(gc[1].x(), gc[1].y(), gc[1].z()),
0041 HepGeom::Point3D<double>(gc[2].x(), gc[2].y(), gc[2].z()),
0042 HepGeom::Point3D<double>(gc[3].x(), gc[3].y(), gc[3].z()),
0043 HepGeom::Point3D<double>(gc[4].x(), gc[4].y(), gc[4].z()),
0044 HepGeom::Point3D<double>(gc[5].x(), gc[5].y(), gc[5].z()),
0045 HepGeom::Point3D<double>(gc[6].x(), gc[6].y(), gc[6].z()),
0046 HepGeom::Point3D<double>(gc[7].x(), gc[7].y(), gc[7].z())};
0047 const HepGeom::Point3D<double> ctr(0.125 * (cv[0] + cv[1] + cv[2] + cv[3] + cv[4] + cv[5] + cv[6] + cv[7]));
0048 const double dCut2(bCut2 / 4.);
0049 if (dCut2 > line.dist2(ctr))
0050 {
0051 LogDebug("CaloCellCrossing") << "** 2nd cut: ctr=" << ctr << ", dist=" << line.dist(ctr);
0052 static const unsigned int nc[6][4] = {
0053 {0, 1, 2, 3}, {0, 4, 5, 1}, {0, 4, 7, 3}, {6, 7, 4, 5}, {6, 2, 3, 7}, {6, 2, 1, 5}};
0054 for (unsigned int face(0); face != 6; ++face) {
0055 const unsigned int* ic(&nc[face][0]);
0056 const HepGeom::Plane3D<double> pl(cv[ic[0]], cv[ic[1]], cv[ic[2]]);
0057 bool parallel;
0058 const HepGeom::Point3D<double> pt(line.point(pl, parallel));
0059 LogDebug("CaloCellCrossing") << "***Face: " << face << ", pt=" << pt;
0060 if (!parallel) {
0061 LogDebug("CaloCellCrossing") << "Not parallel";
0062 const HepLine3D la(cv[ic[0]], cv[ic[1]], eps);
0063 const HepLine3D lb(cv[ic[2]], cv[ic[3]], eps);
0064 LogDebug("CaloCellCrossing") << "la.point=" << la.point(pt);
0065
0066
0067
0068 if (eps > (la.point(pt) - pt).dot((lb.point(pt) - pt))) {
0069 const HepLine3D lc(cv[ic[0]], cv[ic[3]], eps);
0070 const HepLine3D ld(cv[ic[1]], cv[ic[2]], eps);
0071
0072
0073 if (eps > (lc.point(pt) - pt).dot((ld.point(pt) - pt))) {
0074 if (0 == found) {
0075 ++found;
0076 m_detId.emplace_back(dId);
0077 m_ctr.emplace_back(GlobalPoint(ctr.x(), ctr.y(), ctr.z()));
0078 m_entr.emplace_back(GlobalPoint(pt.x(), pt.y(), pt.z()));
0079 } else {
0080 if (1 == found) {
0081 ++found;
0082 m_exit.emplace_back(GlobalPoint(pt.x(), pt.y(), pt.z()));
0083 } else {
0084 const double dist1(
0085 (pt - HepGeom::Point3D<double>(m_entr.back().x(), m_entr.back().y(), m_entr.back().z())).mag());
0086 const double dist2(
0087 (pt - HepGeom::Point3D<double>(m_exit.back().x(), m_exit.back().y(), m_exit.back().z())).mag());
0088 if (eps < dist1 && eps < dist2) {
0089 edm::LogWarning("CaloCellCrossing")
0090 << "********For DetId = " << dId << " distances too big: " << dist1 << ", " << dist2;
0091 }
0092 }
0093 }
0094 }
0095 }
0096 }
0097 }
0098 }
0099 }
0100 assert(2 >= found);
0101 if (1 == found)
0102 m_exit.emplace_back(m_entr.back());
0103 }
0104
0105 assert(m_detId.size() == m_entr.size() && m_detId.size() == m_ctr.size() && m_detId.size() == m_exit.size());
0106
0107 m_len.reserve(m_entr.size());
0108 for (unsigned int i(0); i != m_entr.size(); ++i) {
0109 m_len.emplace_back((m_exit[i] - m_entr[i]).mag());
0110 }
0111 }