CaloGeometry/src/CaloCellCrossing.cc

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))  // first loose cut
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))  // tighter cut
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             //        const double dot (  ( la.point( pt ) - pt ).dot( ( lb.point( pt ) - pt ) ) ) ;
0067             //        LogDebug("CaloCellCrossing")<<"***Dot1="<<dot;
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               //             const double dot (  ( lc.point( pt ) - pt ).dot( ( ld.point( pt ) - pt ) ) ) ;
0072               //             LogDebug("CaloCellCrossing")<<"***Dot2="<<dot;
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 }