CaloCellCrossing.cc

Go to the documentation of this file.

#include "Geometry/CaloGeometry/interface/Line3D.h"
#include "Geometry/CaloGeometry/interface/CaloCellCrossing.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

CaloCellCrossing::CaloCellCrossing( const GlobalPoint&              gp ,
                    const GlobalVector&             gv ,
                    const CaloCellCrossing::DetIds* di ,
                    const CaloSubdetectorGeometry*  sg ,
                    DetId::Detector                 det ,
                    int                             subdet,
                    double                          small,
                    bool                            onewayonly ) :
   m_gp ( gp ) ,
   m_gv ( gv ) 
{
   const double eps ( fabs( small ) > 1.e-15 ? fabs( small ) : 1.e-10 ) ;
   static unsigned int k_rlen ( 30 ) ;
   m_detId.reserve( k_rlen ) ;
   m_ctr  .reserve( k_rlen ) ;
   m_entr .reserve( k_rlen ) ;
   m_exit .reserve( k_rlen ) ;
//------------------------------------------------------------
   const HepLine3D line ( HepGeom::Point3D<double> (  gp.x(), gp.y(), gp.z() ), 
              HepGeom::Vector3D<double> ( gv.x(), gv.y(), gv.z() ), eps ) ;

   LogDebug("CaloCellCrossing") << "*** Line: pt=" << line.pt() << ", unitvec=" << line.uv();

   const DetIds& ids ( nullptr == di ? sg->getValidDetIds( det, subdet ) : *di ) ;
//------------------------------------------------------------
   for(auto dId : ids)
   {
      unsigned int found ( 0 ) ;
      auto cg = sg->getGeometry( dId );
      const CaloCellGeometry::CornersVec& gc ( cg->getCorners() ) ;
      const HepGeom::Point3D<double>  fr ( cg->getPosition().x(),
                       cg->getPosition().y(),
                       cg->getPosition().z()  ) ;
      const double bCut2 ( ( gc[0] - gc[6] ).mag2() ) ;

      if( ( !onewayonly ||
        eps < HepGeom::Vector3D<double> ( fr - line.pt() ).dot( line.uv() ) ) &&
      bCut2 > line.dist2( fr ) ) // first loose cut
      {
         LogDebug("CaloCellCrossing") << "*** fr=" << fr << ", bCut =" << sqrt(bCut2) << ", dis=" << line.dist(fr);
     const HepGeom::Point3D<double>  cv[8] = 
        { HepGeom::Point3D<double> ( gc[0].x(), gc[0].y(), gc[0].z() ) ,
          HepGeom::Point3D<double> ( gc[1].x(), gc[1].y(), gc[1].z() ) ,
          HepGeom::Point3D<double> ( gc[2].x(), gc[2].y(), gc[2].z() ) ,
          HepGeom::Point3D<double> ( gc[3].x(), gc[3].y(), gc[3].z() ) ,
          HepGeom::Point3D<double> ( gc[4].x(), gc[4].y(), gc[4].z() ) ,
          HepGeom::Point3D<double> ( gc[5].x(), gc[5].y(), gc[5].z() ) ,
          HepGeom::Point3D<double> ( gc[6].x(), gc[6].y(), gc[6].z() ) ,
          HepGeom::Point3D<double> ( gc[7].x(), gc[7].y(), gc[7].z() ) } ;
     const HepGeom::Point3D<double>  ctr ( 0.125*(cv[0]+cv[1]+cv[2]+cv[3]+
                       cv[4]+cv[5]+cv[6]+cv[7]) ) ;
     const double dCut2 ( bCut2/4. ) ;
     if( dCut2 > line.dist2( ctr ) ) // tighter cut
     {
            LogDebug("CaloCellCrossing") << "** 2nd cut: ctr=" << ctr
                     << ", dist=" << line.dist(ctr);
        static const unsigned int nc[6][4] = 
           { { 0,1,2,3 }, { 0,4,5,1 }, { 0,4,7,3 },
         { 6,7,4,5 }, { 6,2,3,7 }, { 6,2,1,5 } } ;
        for( unsigned int face ( 0 ) ; face != 6 ; ++face )
        {
           const unsigned int* ic ( &nc[face][0] ) ;
           const HepGeom::Plane3D<double>  pl ( cv[ic[0]], cv[ic[1]], cv[ic[2]] ) ;
           bool parallel ;
           const HepGeom::Point3D<double>  pt ( line.point( pl, parallel ) ) ;
               LogDebug("CaloCellCrossing")<<"***Face: "<<face<<", pt="<<pt;
           if( !parallel )
           {
          LogDebug("CaloCellCrossing")<<"Not parallel";
          const HepLine3D la ( cv[ic[0]], cv[ic[1]], eps ) ;
          const HepLine3D lb ( cv[ic[2]], cv[ic[3]], eps ) ;
          LogDebug("CaloCellCrossing")<<"la.point="<<la.point(pt);

//        const double dot (  ( la.point( pt ) - pt ).dot( ( lb.point( pt ) - pt ) ) ) ;
//        LogDebug("CaloCellCrossing")<<"***Dot1="<<dot;
          if( eps > ( la.point( pt ) - pt ).dot( ( lb.point( pt ) - pt ) ) )
          {
             const HepLine3D lc ( cv[ic[0]], cv[ic[3]], eps ) ;
             const HepLine3D ld ( cv[ic[1]], cv[ic[2]], eps ) ;
//           const double dot (  ( lc.point( pt ) - pt ).dot( ( ld.point( pt ) - pt ) ) ) ;
//           LogDebug("CaloCellCrossing")<<"***Dot2="<<dot;
             if( eps > ( lc.point( pt ) - pt ).dot( ( ld.point( pt ) - pt ) ) )
             {
            if( 0 == found )
            {
               ++found ;
               m_detId.emplace_back( dId ) ;
               m_ctr .emplace_back( GlobalPoint( ctr.x(), ctr.y(), ctr.z() ) ) ;
               m_entr.emplace_back( GlobalPoint(  pt.x(),  pt.y(),  pt.z() ) ) ;
            }
            else
            {
               if( 1 == found )
               {
                  ++found ;
                  m_exit.emplace_back( GlobalPoint( pt.x(), pt.y(), pt.z() ) ) ;
               } 
               else
               {
                  const double dist1 (
                 ( pt - HepGeom::Point3D<double> ( m_entr.back().x(),
                            m_entr.back().y(),
                            m_entr.back().z() ) ).mag() ) ;
                  const double dist2 ( 
                 ( pt - HepGeom::Point3D<double> ( m_exit.back().x(),
                            m_exit.back().y(),
                            m_exit.back().z() ) ).mag() ) ;
                  if( eps < dist1 && 
                  eps < dist2    )
                  {
                 edm::LogWarning("CaloCellCrossing") << "********For DetId = " << dId 
                                     << " distances too big: "
                                     << dist1 << ", " << dist2;
                  
                  }
               }
            }
             }
          }
           }
        }
     }
      }
      assert( 2 >= found ) ;
      if( 1 == found ) m_exit.emplace_back( m_entr.back() ) ;
   }
//------------------------------------------------------------
   assert( m_detId.size() == m_entr.size() &&
       m_detId.size() == m_ctr .size() &&
       m_detId.size() == m_exit.size()    ) ;

   m_len.reserve( m_entr.size() ) ;
   for( unsigned int i ( 0 ) ; i != m_entr.size() ; ++i )
   {
      m_len.emplace_back( ( m_exit[i] - m_entr[i] ).mag() ) ;
   }
}