00001 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00002 #include <CLHEP/Geometry/Plane3D.h>
00003
00004 typedef CaloCellGeometry::CCGFloat CCGFloat ;
00005 typedef CaloCellGeometry::Pt3D Pt3D ;
00006 typedef CaloCellGeometry::Pt3DVec Pt3DVec ;
00007 typedef CaloCellGeometry::Tr3D Tr3D ;
00008
00009 typedef HepGeom::Vector3D<CCGFloat> Vec3D ;
00010 typedef HepGeom::Plane3D<CCGFloat> Plane3D ;
00011
00012 typedef HepGeom::Vector3D<double> DVec3D ;
00013 typedef HepGeom::Plane3D<double> DPlane3D ;
00014 typedef HepGeom::Point3D<double> DPt3D ;
00015 typedef std::vector<DPt3D> DPt3DVec ;
00016
00017 const float CaloCellGeometry::k_ScaleFromDDDtoGeant ( 0.1 ) ;
00018
00019 CaloCellGeometry::CaloCellGeometry() :
00020 m_refPoint ( 0., 0., 0. ),
00021 m_corners ( ) ,
00022 m_parms ( (CCGFloat*) 0 ), m_eta(0), m_phi(0)
00023 {}
00024
00025
00026 CaloCellGeometry::~CaloCellGeometry()
00027 {}
00028
00029
00030
00031
00032
00033 CaloCellGeometry::CaloCellGeometry( CornersVec::const_reference gp ,
00034 const CornersMgr* mgr,
00035 const CCGFloat* par ) :
00036 m_refPoint ( gp ),
00037 m_corners ( mgr ),
00038 m_parms ( par ),
00039 m_eta(gp.eta()), m_phi(gp.phi())
00040 {}
00041
00042 CaloCellGeometry::CaloCellGeometry( const CornersVec& cv,
00043 const CCGFloat* par ) :
00044 m_refPoint ( 0.25*( cv[0].x() + cv[1].x() + cv[2].x() + cv[3].x() ),
00045 0.25*( cv[0].y() + cv[1].y() + cv[2].y() + cv[3].y() ),
00046 0.25*( cv[0].z() + cv[1].z() + cv[2].z() + cv[3].z() ) ),
00047 m_corners ( cv ),
00048 m_parms ( par ),
00049 m_eta( m_refPoint.eta()), m_phi(m_refPoint.phi())
00050 {}
00051
00052 CaloCellGeometry::CornersVec&
00053 CaloCellGeometry::setCorners() const
00054 {
00055 return m_corners ;
00056 }
00057
00058 const CaloCellGeometry::CornersVec&
00059 CaloCellGeometry::getCorners() const
00060 {
00061 return m_corners ;
00062 }
00063
00064 std::ostream& operator<<( std::ostream& s, const CaloCellGeometry& cell )
00065 {
00066 s << ", Center: " << cell.getPosition() << std::endl;
00067
00068 if( cell.emptyCorners() )
00069 {
00070 s << "Corners vector is empty." << std::endl ;
00071 }
00072 else
00073 {
00074 const CaloCellGeometry::CornersVec& corners ( cell.getCorners() ) ;
00075 for ( unsigned int i ( 0 ) ; i != corners.size() ; ++i )
00076 {
00077 s << "Corner: " << corners[ i ] << std::endl;
00078 }
00079 }
00080 return s ;
00081 }
00082
00083 void
00084 CaloCellGeometry::getTransform( Tr3D& tr , Pt3DVec* lptr ) const
00085 {
00086 const GlobalPoint& p ( CaloCellGeometry::getPosition() ) ;
00087 const Pt3D gFront ( p.x(), p.y(), p.z() ) ;
00088 const DPt3D dgFront ( p.x(), p.y(), p.z() ) ;
00089
00090 Pt3D lFront ;
00091 assert( 0 != param() ) ;
00092 Pt3DVec lc ( 8, Pt3D(0,0,0) ) ;
00093 vocalCorners( lc, param(), lFront ) ;
00094
00095 DPt3D dlFront ( lFront.x(), lFront.y(), lFront.z() ) ;
00096
00097 const Pt3D lBack ( 0.25*(lc[4]+lc[5]+lc[6]+lc[7]) ) ;
00098 const DPt3D dlBack ( lBack.x(), lBack.y(), lBack.z() ) ;
00099
00100 const Pt3D dlOne ( lc[0].x(), lc[0].y(), lc[0].z() ) ;
00101
00102 const CornersVec& cor ( getCorners() ) ;
00103 DPt3DVec dkor ( 8, DPt3D(0,0,0) ) ;
00104 for( unsigned int i ( 0 ) ; i != 8 ; ++i )
00105 {
00106 dkor[i] = DPt3D ( cor[i].x(), cor[i].y(), cor[i].z() ) ;
00107 }
00108
00109 DPt3D dgBack ( 0.25*( dkor[4]+dkor[5]+dkor[6]+dkor[7] ) ) ;
00110
00111 const DVec3D dgAxis ( (dgBack-dgFront).unit() ) ;
00112
00113 dgBack = ( dgFront + (dlBack-dlFront).mag()*dgAxis ) ;
00114 const DPt3D dgOneT ( dgFront + ( dlOne - dlFront ).mag()*( dkor[0] - dgFront ).unit() ) ;
00115
00116 const double dlangle ( ( dlBack - dlFront).angle( dlOne - dlFront ) ) ;
00117 const double dgangle ( ( dgBack - dgFront).angle( dgOneT- dgFront ) ) ;
00118 const double ddangle ( dlangle - dgangle ) ;
00119
00120 const DPlane3D dgPl ( dgFront, dgOneT, dgBack ) ;
00121 const DPt3D dp2 ( dgFront + dgPl.normal().unit() ) ;
00122
00123 const DPt3D dgOne ( dgFront + HepGeom::Rotate3D( -ddangle, dgFront, dp2 )*
00124 DVec3D ( dgOneT - dgFront ) ) ;
00125
00126 tr = Tr3D( dlFront , dlBack , dlOne ,
00127 dgFront , dgBack , dgOne ) ;
00128
00129 if( 0 != lptr ) (*lptr) = lc ;
00130 }
00131
00132 const float*
00133 CaloCellGeometry::checkParmPtr(
00134 const std::vector<float>& vv ,
00135 CaloCellGeometry::ParVecVec& pvv )
00136 {
00137 const float* pP ( 0 ) ;
00138
00139 for( unsigned int ii ( 0 ) ; ii != pvv.size() ; ++ii )
00140 {
00141 const ParVec& v ( pvv[ii] ) ;
00142 assert( v.size() == vv.size() ) ;
00143
00144 bool same ( true ) ;
00145 for( unsigned int j ( 0 ) ; j != vv.size() ; ++j )
00146 {
00147 same = same && ( fabs( vv[j] - v[j] )<1.e-6 ) ;
00148 if( !same ) break ;
00149 }
00150 if( same )
00151 {
00152 pP = &(*v.begin()) ;
00153 break ;
00154 }
00155 }
00156 return pP ;
00157 }
00158
00159 const float*
00160 CaloCellGeometry::getParmPtr(
00161 const std::vector<float>& vv ,
00162 CaloCellGeometry::ParMgr* mgr ,
00163 CaloCellGeometry::ParVecVec& pvv )
00164 {
00165 const float* pP ( checkParmPtr( vv, pvv ) ) ;
00166
00167 if( 0 == pP )
00168 {
00169 pvv.push_back( ParVec( mgr ) ) ;
00170 ParVec& back ( pvv.back() ) ;
00171 for( unsigned int i ( 0 ) ; i != vv.size() ; ++i )
00172 {
00173 back[i] = vv[i] ;
00174 }
00175 pP = &(*back.begin()) ;
00176 }
00177 return pP ;
00178 }
00179
00180 bool
00181 CaloCellGeometry::inside( const GlobalPoint& point ) const
00182 {
00183 bool ans ( false ) ;
00184 const Pt3D p ( point.x(), point.y(), point.z() ) ;
00185 const CornersVec& cog ( getCorners() ) ;
00186 Pt3D co[8] ;
00187 for( unsigned int i ( 0 ) ; i != 8 ; ++i )
00188 {
00189 co[i] = Pt3D ( cog[i].x(), cog[i].y(), cog[i].z() ) ;
00190 }
00191
00192 const Plane3D AA ( co[0], co[1], co[2] ) ;
00193 const Plane3D BB ( co[6], co[5], co[4] ) ;
00194
00195 if( AA.distance(p)*BB.distance(p) >= 0 )
00196 {
00197 const Plane3D CC ( co[0], co[4], co[5] ) ;
00198 const Plane3D DD ( co[2], co[6], co[7] ) ;
00199 if( CC.distance(p)*DD.distance(p) >= 0 )
00200 {
00201 const Plane3D EE ( co[3], co[7], co[4] ) ;
00202 const Plane3D FF ( co[1], co[5], co[6] ) ;
00203 if( EE.distance(p)*FF.distance(p) >= 0 )
00204 {
00205 ans = true ;
00206 }
00207 }
00208 }
00209 return ans ;
00210 }