CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/Geometry/CaloGeometry/src/CaloCellGeometry.cc

Go to the documentation of this file.
00001 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00002 #include <CLHEP/Geometry/Plane3D.h>
00003 
00004 const float CaloCellGeometry::k_ScaleFromDDDtoGeant ( 0.1 ) ;
00005 
00006 const CaloCellGeometry::CornersVec&
00007 CaloCellGeometry::getCorners() const
00008 {
00009    return m_corners ;
00010 }
00011 
00012 std::ostream& operator<<( std::ostream& s, const CaloCellGeometry& cell ) 
00013 {
00014    s << ", Center: " <<  cell.getPosition() << std::endl;
00015 
00016    if( cell.emptyCorners() )
00017    {
00018       s << "Corners vector is empty." << std::endl ;
00019    }
00020    else
00021    {
00022       const CaloCellGeometry::CornersVec& corners ( cell.getCorners() ) ;
00023       for ( unsigned int i ( 0 ) ; i != corners.size() ; ++i ) 
00024       {
00025          s << "Corner: " << corners[ i ] << std::endl;
00026       }
00027    }
00028    return s ;
00029 }
00030 
00031 HepGeom::Transform3D
00032 CaloCellGeometry::getTransform( std::vector<HepGeom::Point3D<double> >* lptr ) const 
00033 {
00034    const GlobalPoint& p ( CaloCellGeometry::getPosition() ) ;
00035    const HepGeom::Point3D<double>    gFront ( p.x(), p.y(), p.z() ) ;
00036 
00037    HepGeom::Point3D<double>  lFront ;
00038    assert(                               0 != param() ) ;
00039    std::vector<HepGeom::Point3D<double> > lc ( vocalCorners( param(), lFront ) ) ;
00040 
00041    HepGeom::Point3D<double>  lBack  ( 0.25*(lc[4]+lc[5]+lc[6]+lc[7]) ) ;
00042 
00043    const HepGeom::Point3D<double>  lOne  ( lc[0] ) ;
00044 
00045    const CornersVec& cor ( getCorners() ) ;
00046    std::vector<HepGeom::Point3D<double> > kor ( 8, HepGeom::Point3D<double> (0,0,0) ) ;
00047    for( unsigned int i ( 0 ) ; i != 8 ; ++i )
00048    {
00049       kor[i] = HepGeom::Point3D<double> ( cor[i].x(), cor[i].y(), cor[i].z() ) ;
00050    }
00051 
00052    HepGeom::Point3D<double>  gBack ( 0.25*( kor[4]+kor[5]+kor[6]+kor[7] ) ) ;
00053 
00054    const HepGeom::Vector3D<double>  gAxis ( (gBack-gFront).unit() ) ;
00055 
00056    gBack = ( gFront + (lBack-lFront).mag()*gAxis ) ;
00057    const HepGeom::Point3D<double>  gOneT ( gFront + ( lOne - lFront ).mag()*( kor[0] - gFront ).unit() ) ;
00058 
00059    const double langle ( ( lBack - lFront).angle( lOne - lFront ) ) ;
00060    const double gangle ( ( gBack - gFront).angle( gOneT- gFront ) ) ;
00061    const double dangle ( langle - gangle ) ;
00062 
00063    const HepGeom::Plane3D<double>  gPl (  gFront, gOneT, gBack ) ;
00064    const HepGeom::Point3D<double>  p2  ( gFront + gPl.normal().unit() ) ;
00065 
00066    const HepGeom::Point3D<double>  gOne ( gFront + HepGeom::Rotate3D( -dangle, gFront, p2 )*
00067                            HepGeom::Vector3D<double> ( gOneT - gFront ) ) ;
00068 
00069    const HepGeom::Transform3D tr ( lFront , lBack , lOne ,
00070                              gFront , gBack , gOne    ) ;
00071 
00072    if( 0 != lptr ) (*lptr) = lc ;
00073 
00074    return tr ;
00075 }
00076 
00077 const double* 
00078 CaloCellGeometry::checkParmPtr(
00079    const std::vector<double>&   vv  ,
00080    CaloCellGeometry::ParVecVec& pvv  )
00081 {
00082    const double* pP ( 0 ) ;
00083 
00084    for( unsigned int ii ( 0 ) ; ii != pvv.size() ; ++ii )
00085    {
00086       const ParVec& v ( pvv[ii] ) ;
00087       assert( v.size() == vv.size() ) ;
00088 
00089       bool same ( true ) ;
00090       for( unsigned int j ( 0 ) ; j != vv.size() ; ++j )
00091       {
00092          same = same && ( fabs( vv[j] - v[j] )<1.e-6 ) ;
00093          if( !same ) break ;
00094       }
00095       if( same )
00096       {
00097          pP = &(*v.begin()) ;
00098          break ;
00099       }
00100    }
00101    return pP ;
00102 }
00103 
00104 const double* 
00105 CaloCellGeometry::getParmPtr(
00106    const std::vector<double>&   vv  ,
00107    CaloCellGeometry::ParMgr*    mgr ,
00108    CaloCellGeometry::ParVecVec& pvv  )
00109 {
00110    const double* pP ( checkParmPtr( vv, pvv ) ) ;
00111 
00112    if( 0 == pP )
00113    {
00114       pvv.push_back( ParVec( mgr ) ) ;
00115       ParVec& back ( pvv.back() ) ;
00116       for( unsigned int i ( 0 ) ; i != vv.size() ; ++i )
00117       {
00118          back[i] = vv[i] ;
00119       }
00120       pP = &(*back.begin()) ;
00121    }
00122    return pP ;
00123 }
00124 
00125 bool 
00126 CaloCellGeometry::inside( const GlobalPoint& point ) const
00127 {
00128    bool ans ( false ) ;
00129    const HepGeom::Point3D<double>  p ( point.x(), point.y(), point.z() ) ;
00130    const CornersVec& cog ( getCorners() ) ;
00131    HepGeom::Point3D<double>  co[8] ;
00132    for( unsigned int i ( 0 ) ; i != 8 ; ++i )
00133    {
00134       co[i] = HepGeom::Point3D<double> ( cog[i].x(), cog[i].y(), cog[i].z() ) ;
00135    }
00136 
00137    const HepGeom::Plane3D<double>  AA ( co[0], co[1], co[2] ) ; // z<0
00138    const HepGeom::Plane3D<double>  BB ( co[6], co[5], co[4] ) ; // z>0
00139 
00140    if( AA.distance(p)*BB.distance(p) >= 0 )
00141    {
00142      const HepGeom::Plane3D<double>  CC ( co[0], co[4], co[5] ) ; // x<0
00143      const HepGeom::Plane3D<double>  DD ( co[2], co[6], co[7] ) ; // x>0
00144      if( CC.distance(p)*DD.distance(p) >= 0 )
00145      {
00146        const HepGeom::Plane3D<double>  EE ( co[3], co[7], co[4] ) ; // y<0
00147        const HepGeom::Plane3D<double>  FF ( co[1], co[5], co[6] ) ; // y>0
00148        if( EE.distance(p)*FF.distance(p) >= 0 )
00149        {
00150            ans = true ;
00151        }
00152       }
00153    }
00154    return ans ;
00155 }