From 3fd5fc3d333b7600b5cba1e2b74f680ab028b258 Mon Sep 17 00:00:00 2001 From: Frank Gaede <frank.gaede@desy.de> Date: Tue, 27 May 2014 09:28:11 +0000 Subject: [PATCH] - implemented u,v,normal at point for cylindircal surfaces - changed ISurface interface for u(),v(),normal() to return Vector3D rather than const Vector3d& - added some documentation for Surface methods --- DDRec/CMakeLists.txt | 2 + DDRec/include/DDRec/Surface.h | 119 +++++++++++++++++++---- DDRec/src/Surface.cpp | 95 +++++++++++++++++- DDSurfaces/include/DDSurfaces/ISurface.h | 6 +- DDSurfaces/include/DDSurfaces/Vector3D.h | 17 +++- DDTest/src/test_surface.cc | 21 +++- examples/ILDExDet/src/ILDExTPC_geo.cpp | 6 +- examples/ILDExDet/src/ILDExVXD_geo.cpp | 2 +- 8 files changed, 232 insertions(+), 36 deletions(-) diff --git a/DDRec/CMakeLists.txt b/DDRec/CMakeLists.txt index efbad2154..062dd8590 100644 --- a/DDRec/CMakeLists.txt +++ b/DDRec/CMakeLists.txt @@ -19,6 +19,8 @@ if(DD4HEP_WITH_GEAR) set(sources ${sources} src/gear/DDGear.cpp) + add_definitions("-D DD4HEP_WITH_GEAR") + endif() #------------------------------------------------------------------------------------- diff --git a/DDRec/include/DDRec/Surface.h b/DDRec/include/DDRec/Surface.h index c105ef82c..d41216f18 100644 --- a/DDRec/include/DDRec/Surface.h +++ b/DDRec/include/DDRec/Surface.h @@ -75,16 +75,17 @@ namespace DD4hep { SurfaceMaterial _innerMat ; SurfaceMaterial _outerMat ; + /// default c'tor. SurfaceData(); + /// Standard c'tor for initialization. SurfaceData( SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D u ,Vector3D v ,Vector3D n ,Vector3D o ) ; /// Default destructor virtual ~SurfaceData() {} - /// Copy the object - + /// Copy the from object void copy(const SurfaceData& c) { _type = c._type ; _u = c._u ; @@ -111,17 +112,29 @@ namespace DD4hep { Geometry::Volume _vol ; + /// setter for daughter classes + virtual void setU(const Vector3D& u) { object<SurfaceData>()._u = u ; } + + /// setter for daughter classes + virtual void setV(const Vector3D& v) { object<SurfaceData>()._v = v ; } + + /// setter for daughter classes + virtual void setNormal(const Vector3D& n) { object<SurfaceData>()._n = n ; } + + public: virtual ~VolSurface() {} + ///default c'tor VolSurface() { } - + + /// Standrad c'tor for initialization. VolSurface( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D u ,Vector3D v ,Vector3D n , Vector3D o = Vector3D(0.,0.,0.) ) ; + /// the volume to which this surface is attached. Geometry::Volume volume() const { return _vol ; } - /// The id of this surface - always 0 for VolSurfaces virtual long64 id() const { return 0 ; } @@ -134,13 +147,13 @@ namespace DD4hep { //==== geometry ==== /** First direction of measurement U */ - virtual const Vector3D& u( const Vector3D& point = Vector3D() ) const { point.x() ; return object<SurfaceData>()._u ; } + virtual Vector3D u( const Vector3D& point = Vector3D() ) const { point.x() ; return object<SurfaceData>()._u ; } /** Second direction of measurement V */ - virtual const Vector3D& v(const Vector3D& point = Vector3D() ) const { point.x() ; return object<SurfaceData>()._v ; } + virtual Vector3D v(const Vector3D& point = Vector3D() ) const { point.x() ; return object<SurfaceData>()._v ; } /// Access to the normal direction at the given point - virtual const Vector3D& normal(const Vector3D& point = Vector3D() ) const { point.x() ; return object<SurfaceData>()._n ; } + virtual Vector3D normal(const Vector3D& point = Vector3D() ) const { point.x() ; return object<SurfaceData>()._n ; } /** Get Origin of local coordinate system on surface */ virtual const Vector3D& origin() const { return object<SurfaceData>()._o ;} @@ -189,11 +202,13 @@ namespace DD4hep { VolSurfaceList() {} - // required c'tors for extension mechanism + // required c'tor for extension mechanism VolSurfaceList(Geometry::DetElement& det){ // det.addExtension<VolSurfaceList>( this ) ; } + + // required c'tor for extension mechanism VolSurfaceList(const VolSurfaceList& vsl, Geometry::DetElement& det ){ //fixme: this causes a seg fault ... @@ -215,6 +230,9 @@ namespace DD4hep { } } ; + /** Helper function for accessing the list assigned to a DetElement - attaches + * empty list if needed. + */ VolSurfaceList* volSurfaceList( Geometry::DetElement& det ) ; //====================================================================================================== @@ -228,10 +246,13 @@ namespace DD4hep { public: + ///default c'tor VolPlane() : VolSurface() { } + /// copy c'tor VolPlane(const VolSurface& vs ) : VolSurface( vs ) { } + /// standard c'tor with all necessary arguments - origin is (0,0,0) if not given. VolPlane( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D u ,Vector3D v ,Vector3D n , Vector3D o = Vector3D(0.,0.,0.) ) : @@ -239,6 +260,8 @@ namespace DD4hep { object<SurfaceData>()._type.setProperty( SurfaceType::Plane , true ) ; object<SurfaceData>()._type.setProperty( SurfaceType::Cylinder , false ) ; + object<SurfaceData>()._type.checkParallelToZ( *this ) ; + object<SurfaceData>()._type.checkOrthogonalToZ( *this ) ; } @@ -262,20 +285,43 @@ namespace DD4hep { public: + /// default c'tor VolCylinder() : VolSurface() { } + /// copy c'tor VolCylinder(const VolSurface& vs ) : VolSurface( vs ) { } - VolCylinder( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, - Vector3D u ,Vector3D v ,Vector3D n , Vector3D o = Vector3D(0.,0.,0.) ) : - - VolSurface( vol, type, thickness_inner, thickness_outer, u,v,n,o ) { + /** The standard constructor. The origin vector points to the origin of the coordinate system on the cylinder, + * its rho defining the radius of the cylinder. The measurement direction u is set to be (0,0,1), the normal is + * chosen to be parallel to the origin vector and v = n X u. + */ + VolCylinder( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D origin ) ; - object<SurfaceData>()._type.setProperty( SurfaceType::Plane , false ) ; - object<SurfaceData>()._type.setProperty( SurfaceType::Cylinder , true ) ; - } + //fg: for now we don't allow to set u and v for cylinders + // VolCylinder( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, + // Vector3D u ,Vector3D v ,Vector3D n , Vector3D o = Vector3D(0.,0.,0.) ) : + // VolSurface( vol, type, thickness_inner, thickness_outer, u,v,n,o ) { + // object<SurfaceData>()._type.setProperty( SurfaceType::Plane , false ) ; + // object<SurfaceData>()._type.setProperty( SurfaceType::Cylinder , true ) ; + // } + + /** First direction of measurement U - rotated to point projected onto the cylinder. + * No check is done whether the point actually is on the cylinder surface + */ + virtual Vector3D u( const Vector3D& point = Vector3D() ) const ; + + /** Second direction of measurement V - rotated to point projected onto the cylinder. + * No check is done whether the point actually is on the cylinder surface + */ + virtual Vector3D v(const Vector3D& point = Vector3D() ) const ; + + /** The normal direction at the given point, projected onto the cylinder. + * No check is done whether the point actually is on the cylinder surface + */ + virtual Vector3D normal(const Vector3D& point = Vector3D() ) const ; + /** Distance to surface */ virtual double distance(const Vector3D& point ) const ; @@ -309,15 +355,19 @@ namespace DD4hep { Vector3D _n ; Vector3D _o ; + /// default c'tor Surface() :_det( Geometry::DetElement() ), _volSurf( VolSurface() ), _wtM( 0 ) , _id( 0) { } public: virtual ~Surface() {} + /** Standard c'tor initializes the surface from the parameters of the VolSurface and the + * transform (placement) of the corresponding volume, if found in DetElement + */ Surface( Geometry::DetElement det, VolSurface volSurf ) ; - /// The id of this surface - corresponds to DetElement id ( or'ed with the placement ids ) + /// The id of this surface - corresponds to DetElement id. virtual long64 id() const { return _id ; } /** properties of the surface encoded in Type. @@ -325,21 +375,23 @@ namespace DD4hep { */ virtual const SurfaceType& type() const { return _type ; } + /// The volume that has the surface attached. Geometry::Volume volume() const { return _volSurf.volume() ; } + /// The VolSurface attched to the volume. VolSurface volSurface() const { return _volSurf ; } //==== geometry ==== /** First direction of measurement U */ - virtual const Vector3D& u( const Vector3D& point = Vector3D() ) const { point.x() ; return _u ; } + virtual Vector3D u( const Vector3D& point = Vector3D() ) const { point.x() ; return _u ; } /** Second direction of measurement V */ - virtual const Vector3D& v(const Vector3D& point = Vector3D() ) const { point.x() ; return _v ; } + virtual Vector3D v(const Vector3D& point = Vector3D() ) const { point.x() ; return _v ; } /// Access to the normal direction at the given point - virtual const Vector3D& normal(const Vector3D& point = Vector3D() ) const { point.x() ; return _n ; } + virtual Vector3D normal(const Vector3D& point = Vector3D() ) const { point.x() ; return _n ; } /** Get Origin of local coordinate system on surface */ virtual const Vector3D& origin() const { return _o ;} @@ -380,14 +432,35 @@ namespace DD4hep { //====================================================================================================== + /** Specialization of Surface for cylinders. Provides acces to the cylinder radius and + * implements the access to the rotated surface vectors for points on the cylinder. + * @author F.Gaede, DESY + * @date May, 10 2014 + */ class CylinderSurface: public Surface, public ICylinder { public: + ///Standard c'tor. CylinderSurface( Geometry::DetElement det, VolSurface volSurf ) : Surface( det, volSurf ) { } + /** First direction of measurement U - rotated to point projected onto the cylinder. + * No check is done whether the point actually is on the cylinder surface + */ + virtual Vector3D u( const Vector3D& point = Vector3D() ) const ; + + /** Second direction of measurement V - rotated to point projected onto the cylinder. + * No check is done whether the point actually is on the cylinder surface + */ + virtual Vector3D v(const Vector3D& point = Vector3D() ) const ; + + /** The normal direction at the given point - rotated to point projected onto the cylinder. + * No check is done whether the point actually is on the cylinder surface + */ + virtual Vector3D normal(const Vector3D& point = Vector3D() ) const ; + + /// the radius of the cylinder (rho of the origin vector) virtual double radius() const { - return _volSurf.origin().rho() ; } } ; @@ -404,18 +477,22 @@ namespace DD4hep { bool _isOwner ; public: + /// defaul c'tor - allow to set ownership for surfaces SurfaceList(bool isOwner=false ) : _isOwner( isOwner ) {} + /// copy c'tor SurfaceList(const SurfaceList& other ) : std::list< Surface* >( other ), _isOwner( false ){} - // required c'tors for extension mechanism + /// required c'tor for extension mechanism SurfaceList(const Geometry::DetElement& ){ // anything to do here ? } + /// required c'tor for extension mechanism SurfaceList(const SurfaceList& ,const Geometry::DetElement& ){ // anything to do here ? } + /// d'tor deletes all owned surfaces virtual ~SurfaceList(){ if( _isOwner ) { diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp index 5b02db720..bdebf7b5c 100644 --- a/DDRec/src/Surface.cpp +++ b/DDRec/src/Surface.cpp @@ -55,12 +55,13 @@ namespace DD4hep { } - + + VolSurface::VolSurface( Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, - Vector3D u ,Vector3D v ,Vector3D n ,Vector3D o ) : + Vector3D u ,Vector3D v ,Vector3D n ,Vector3D o ) : Geometry::Handle< SurfaceData >( new SurfaceData( type, thickness_inner ,thickness_outer, u,v,n,o) ) , - + _vol( vol ) { } @@ -94,6 +95,58 @@ namespace DD4hep { } + + //============================================================================================================= + + Vector3D VolCylinder::u( const Vector3D& point ) const { + // for now we just have u const as (0,0,1) + point.x() ; return VolSurface::u() ; + } + + Vector3D VolCylinder::v(const Vector3D& point ) const { + + Vector3D n( 1. , point.phi() , 0. , Vector3D::cylindrical ) ; + + // std::cout << " u : " << u() + // << " n : " << n + // << " u X n :" << u().cross( n ) ; + return u().cross( n ) ; + } + + Vector3D VolCylinder::normal(const Vector3D& point ) const { + + // normal is just given by phi of the point + return Vector3D( 1. , point.phi() , 0. , Vector3D::cylindrical ) ; + } + + + + VolCylinder::VolCylinder( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D o ) : + + VolSurface( vol, type, thickness_inner, thickness_outer, Vector3D() , Vector3D() , Vector3D() , o ) { + + Vector3D u( 0., 0., 1. ) ; + + Vector3D o_rphi( o.x() , o.y() , 0. ) ; + + Vector3D n = o_rphi.unit() ; + + Vector3D v = u.cross( n ) ; + + setU( u ) ; + setV( v ) ; + setNormal( n ) ; + + object<SurfaceData>()._type.setProperty( SurfaceType::Plane , false ) ; + + object<SurfaceData>()._type.setProperty( SurfaceType::Cylinder , true ) ; + + object<SurfaceData>()._type.checkParallelToZ( *this ) ; + + object<SurfaceData>()._type.checkOrthogonalToZ( *this ) ; + } + + /** Distance to surface */ double VolCylinder::distance(const Vector3D& point ) const { @@ -122,7 +175,8 @@ namespace DD4hep { - //==================== + + //================================================================================================================ VolSurfaceList* volSurfaceList( DetElement& det ) { @@ -210,12 +264,14 @@ namespace DD4hep { } - Surface::Surface( Geometry::DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ), _wtM(0) , _id( 0) , _type( _volSurf.type() ) { + Surface::Surface( Geometry::DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ), + _wtM(0) , _id( 0) , _type( _volSurf.type() ) { initialize() ; } + const IMaterial& Surface::innerMaterial() const { SurfaceMaterial& mat = _volSurf->_innerMat ; @@ -646,6 +702,35 @@ namespace DD4hep { } + //================================================================================================================ + Vector3D CylinderSurface::u( const Vector3D& point ) const { + + Vector3D lp , u ; + _wtM->MasterToLocal( point , lp.array() ) ; + const DDSurfaces::Vector3D& lu = _volSurf.u( lp ) ; + _wtM->LocalToMasterVect( lu , u.array() ) ; + return u ; + } + + Vector3D CylinderSurface::v(const Vector3D& point ) const { + Vector3D lp , v ; + _wtM->MasterToLocal( point , lp.array() ) ; + const DDSurfaces::Vector3D& lv = _volSurf.v( lp ) ; + _wtM->LocalToMasterVect( lv , v.array() ) ; + return v ; + } + + Vector3D CylinderSurface::normal(const Vector3D& point ) const { + Vector3D lp , n ; + _wtM->MasterToLocal( point , lp.array() ) ; + const DDSurfaces::Vector3D& ln = _volSurf.normal( lp ) ; + _wtM->LocalToMasterVect( ln , n.array() ) ; + return n ; + } + //================================================================================================================ + + + } // namespace } // namespace diff --git a/DDSurfaces/include/DDSurfaces/ISurface.h b/DDSurfaces/include/DDSurfaces/ISurface.h index fb5070402..1efd20a6d 100644 --- a/DDSurfaces/include/DDSurfaces/ISurface.h +++ b/DDSurfaces/include/DDSurfaces/ISurface.h @@ -37,13 +37,13 @@ namespace DDSurfaces { virtual bool insideBounds(const Vector3D& point, double epsilon=1.e-4) const =0 ; /** First direction of measurement U */ - virtual const Vector3D& u( const Vector3D& point = Vector3D() ) const =0 ; + virtual Vector3D u( const Vector3D& point = Vector3D() ) const =0 ; /** Second direction of measurement V */ - virtual const Vector3D& v(const Vector3D& point = Vector3D() ) const =0 ; + virtual Vector3D v(const Vector3D& point = Vector3D() ) const =0 ; /// Access to the normal direction at the given point - virtual const Vector3D& normal(const Vector3D& point = Vector3D() ) const =0 ; + virtual Vector3D normal(const Vector3D& point = Vector3D() ) const =0 ; /** Get Origin of local coordinate system on surface */ virtual const Vector3D& origin() const =0 ; diff --git a/DDSurfaces/include/DDSurfaces/Vector3D.h b/DDSurfaces/include/DDSurfaces/Vector3D.h index ea273e5d9..89c20bb11 100644 --- a/DDSurfaces/include/DDSurfaces/Vector3D.h +++ b/DDSurfaces/include/DDSurfaces/Vector3D.h @@ -206,6 +206,19 @@ namespace DDSurfaces { } + /** Component wise comparison of two vectors - true if all components differ less than epsilon */ + inline bool isEqual( const Vector3D& b , double epsilon=1e-6) { + + if( std::fabs( x() - b.x() ) < epsilon && + std::fabs( y() - b.y() ) < epsilon && + std::fabs( z() - b.z() ) < epsilon ) + return true; + else + return false; + } + + + // this causes template lookup errors on some machines : // -> use explicit conversion with to<T>() // /** Implicit templated conversion to anything that has a c'tor T(x,y,z) @@ -265,7 +278,7 @@ namespace DDSurfaces { return Vector3D( a.x() - b.x() , a.y() - b.y(), a.z() - b.z() ) ; } - /** Comparison of two vectors */ + /** Exact comparison of two vectors */ inline bool operator==( const Vector3D& a, const Vector3D& b ) { if( a.x() == b.x() && a.y() == b.y() && a.z() == b.z() ) @@ -273,7 +286,7 @@ namespace DDSurfaces { else return false; } - + /** Multiplication with scalar */ inline Vector3D operator*( double s , const Vector3D& v ) { diff --git a/DDTest/src/test_surface.cc b/DDTest/src/test_surface.cc index 69e9271c3..4fb1ddd55 100644 --- a/DDTest/src/test_surface.cc +++ b/DDTest/src/test_surface.cc @@ -136,7 +136,9 @@ int main(int argc, char** argv ){ Vector3D o_radius = Vector3D( radius , 0. , 0. ) ; - VolCylinder surfT( vol , SurfaceType( SurfaceType::Sensitive ), thick/2, thick/2 , u,v,n, o_radius ) ; + VolCylinder surfT( vol , SurfaceType( SurfaceType::Sensitive ), thick/2, thick/2 , o_radius ) ; + + std::cout << " *** cylinder surface : " << surfT << std::endl ; test( surfT.insideBounds( Vector3D( radius * sin(0.75) , radius * cos( 0.75 ) , 49. )) , true , " insideBounds Vector3D( radius * sin(0.75) , radius * cos( 0.75 ) , 49. ) " ) ; @@ -152,6 +154,23 @@ int main(int argc, char** argv ){ , " insideBounds Vector3D( (radius+0.00005) * sin(0.75) , (radius+0.00005) * cos( 0.75 ) , 49. ) " ) ; + Vector3D y( 0. , radius , 42 ) ; + + Vector3D yn = surfT.normal( y ) ; + + bool dummy = yn.isEqual( Vector3D( 0. , 1. , 0 ) ) ; + test( dummy , true , " normal at (0.,radius,42) is Vector3D( 0. , 1. , 0 ) " ) ; + + if( ! dummy ) + std::cout << " ** yn = " << yn << std::endl ; + + + Vector3D yv = surfT.v( y ) ; + dummy = yv.isEqual( Vector3D( -1. , 0. , 0 ) ) ; + test( dummy , true , " v at (0.,radius,42) is Vector3D( -1. , 0. , 0 ) " ) ; + if( ! dummy ) + std::cout << " ** yv = " << yv << std::endl ; + // test surface type: diff --git a/examples/ILDExDet/src/ILDExTPC_geo.cpp b/examples/ILDExDet/src/ILDExTPC_geo.cpp index db0d0f555..7fdb30b8e 100644 --- a/examples/ILDExDet/src/ILDExTPC_geo.cpp +++ b/examples/ILDExDet/src/ILDExTPC_geo.cpp @@ -109,7 +109,7 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { case 2: { tpcData->innerWall=part_det; // add a surface to the det element - VolCylinder surf( part_vol , SurfaceType( SurfaceType::Helper ) , drcyl*.5 , drcyl*.5 , u,v,n , ocyl ) ; + VolCylinder surf( part_vol , SurfaceType( SurfaceType::Helper ) , drcyl*.5 , drcyl*.5 , ocyl ) ; volSurfaceList( part_det )->push_back( surf ) ; } break; @@ -117,7 +117,7 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { case 3: { tpcData->outerWall=part_det; // add a surface to the det element - VolCylinder surf( part_vol , SurfaceType( SurfaceType::Helper ) , drcyl*.5 , drcyl*.5 , u,v,n , ocyl ) ; + VolCylinder surf( part_vol , SurfaceType( SurfaceType::Helper ) , drcyl*.5 , drcyl*.5 , ocyl ) ; volSurfaceList( part_det )->push_back( surf ) ; } break; @@ -165,7 +165,7 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { Vector3D o( r0 + (2*i+1) * dR , 0. , 0. ) ; - VolCylinder surf( gas_volU , SurfaceType(SurfaceType::Sensitive, SurfaceType::Invisible ) , dR , dR , u,v,n ,o ) ; + VolCylinder surf( gas_volU , SurfaceType(SurfaceType::Sensitive, SurfaceType::Invisible ) , dR , dR , o ) ; volSurfaceList( layerDE )->push_back( surf ) ; } diff --git a/examples/ILDExDet/src/ILDExVXD_geo.cpp b/examples/ILDExDet/src/ILDExVXD_geo.cpp index 8ad79b2bd..b0f8e131e 100644 --- a/examples/ILDExDet/src/ILDExVXD_geo.cpp +++ b/examples/ILDExDet/src/ILDExVXD_geo.cpp @@ -69,7 +69,7 @@ static Ref_t create_element(LCDD& lcdd, xml_h e, SensitiveDetector sens) { /// ======== test layer assembly -> results in assembly in assembly and -#define use_layer_assembly 0 // crashes ild_exsimu +#define use_layer_assembly 1 // crashes ild_exsimu #if use_layer_assembly // --- create an assembly and DetElement for the layer -- GitLab