Go to the documentation of this file.00001 #include "Geometry/ForwardGeometry/interface/IdealCastorTrapezoid.h"
00002 #include "CLHEP/Geometry/Plane3D.h"
00003 #include "CLHEP/Geometry/Transform3D.h"
00004 #include <math.h>
00005
00006 namespace calogeom {
00007
00008 std::vector<HepGeom::Point3D<double> >
00009 IdealCastorTrapezoid::localCorners( const double* pv ,
00010 HepGeom::Point3D<double> & ref )
00011 {
00012 assert( 0 != pv ) ;
00013
00014 const double dxl ( pv[0] ) ;
00015 const double dxh ( pv[1] ) ;
00016 const double dh ( pv[2] ) ;
00017 const double dz ( pv[3] ) ;
00018 const double an ( pv[4] ) ;
00019 const double dx ( ( dxl +dxh )/2. ) ;
00020 const double dy ( dh*sin(an) ) ;
00021 const double dhz ( dh*cos(an) ) ;
00022 const double dzb ( dz + dhz ) ;
00023 const double dzs ( dz - dhz ) ;
00024
00025 assert( 0 < (dxl*dxh) ) ;
00026
00027 std::vector<HepGeom::Point3D<double> > lc ( 8, HepGeom::Point3D<double> ( 0,0,0) ) ;
00028
00029 lc[ 0 ] = HepGeom::Point3D<double> ( -dx, -dy , dzb ) ;
00030 lc[ 1 ] = HepGeom::Point3D<double> ( -dx, +dy , dzs ) ;
00031 lc[ 2 ] = HepGeom::Point3D<double> ( +2*dxh -dx, +dy , dzs ) ;
00032 lc[ 3 ] = HepGeom::Point3D<double> ( +2*dxl -dx, -dy , dzb ) ;
00033 lc[ 4 ] = HepGeom::Point3D<double> ( -dx, -dy , -dzs ) ;
00034 lc[ 5 ] = HepGeom::Point3D<double> ( -dx, +dy , -dzb ) ;
00035 lc[ 6 ] = HepGeom::Point3D<double> ( +2*dxh -dx, +dy , -dzb ) ;
00036 lc[ 7 ] = HepGeom::Point3D<double> ( +2*dxl -dx, -dy , -dzs ) ;
00037
00038 ref = 0.25*( lc[0] + lc[1] + lc[2] + lc[3] ) ;
00039 return lc ;
00040 }
00041
00042 const CaloCellGeometry::CornersVec&
00043 IdealCastorTrapezoid::getCorners() const
00044 {
00045 const CornersVec& co ( CaloCellGeometry::getCorners() ) ;
00046 if( co.uninitialized() )
00047 {
00048 CaloCellGeometry::CornersVec& corners ( setCorners() ) ;
00049 const GlobalPoint& p ( getPosition() ) ;
00050 const double zsign ( 0 < p.z() ? 1. : -1. ) ;
00051 const HepGeom::Point3D<double> gf ( p.x(), p.y(), p.z() ) ;
00052
00053 HepGeom::Point3D<double> lf ;
00054 const std::vector<HepGeom::Point3D<double> > lc ( localCorners( param(), lf ) ) ;
00055 const HepGeom::Point3D<double> lb ( lf.x() , lf.y() , lf.z() - 2.*dz() ) ;
00056 const HepGeom::Point3D<double> ls ( lf.x() - dx(), lf.y(), lf.z() ) ;
00057
00058
00059 const double fphi ( atan( dx()/( dR() + dy() ) ) ) ;
00060 const HepGeom::Point3D<double> gb ( gf.x() , gf.y() , gf.z() + 2.*zsign*dz() ) ;
00061
00062 const double rho ( dR() + dy() ) ;
00063 const double phi ( gf.phi() + fphi ) ;
00064 const HepGeom::Point3D<double> gs ( rho*cos(phi) ,
00065 rho*sin(phi) ,
00066 gf.z() ) ;
00067
00068 const HepGeom::Transform3D tr ( lf, lb, ls,
00069 gf, gb, gs ) ;
00070
00071 for( unsigned int i ( 0 ) ; i != 8 ; ++i )
00072 {
00073 const HepGeom::Point3D<double> gl ( tr*lc[i] ) ;
00074 corners[i] = GlobalPoint( gl.x(), gl.y(), gl.z() ) ;
00075 }
00076 }
00077 return co ;
00078 }
00079
00080 std::ostream& operator<<( std::ostream& s, const IdealCastorTrapezoid& cell )
00081 {
00082 s << "Center: " << cell.getPosition() << std::endl ;
00083
00084
00085
00086 return s;
00087 }
00088 }