From 0263ae480120e8a062d4e41a6aa5f03a69a9a4d9 Mon Sep 17 00:00:00 2001
From: Frank Gaede <frank.gaede@desy.de>
Date: Thu, 15 Oct 2015 13:04:18 +0000
Subject: [PATCH] - added implemation of conical surfaces - ICone,
VolConeImpl, VolCone, ConeSurface
---
DDRec/include/DDRec/Surface.h | 96 ++++++++++--
DDRec/src/DetectorSurfaces.cpp | 11 +-
DDRec/src/Surface.cpp | 185 +++++++++++++++++++++--
DDSurfaces/include/DDSurfaces/ISurface.h | 45 +++++-
DDSurfaces/include/DDSurfaces/Vector3D.h | 17 +--
5 files changed, 318 insertions(+), 36 deletions(-)
diff --git a/DDRec/include/DDRec/Surface.h b/DDRec/include/DDRec/Surface.h
index cfb96df1f..2a3481c97 100644
--- a/DDRec/include/DDRec/Surface.h
+++ b/DDRec/include/DDRec/Surface.h
@@ -53,6 +53,8 @@ namespace DD4hep {
virtual void setV(const Vector3D& v) ;
/// setter for daughter classes
virtual void setNormal(const Vector3D& n) ;
+ /// setter for daughter classes
+ virtual void setOrigin(const Vector3D& o) ;
public:
@@ -354,6 +356,7 @@ namespace DD4hep {
_type.setProperty( SurfaceType::Plane , true ) ;
_type.setProperty( SurfaceType::Cylinder , false ) ;
+ _type.setProperty( SurfaceType::Cone , false ) ;
_type.checkParallelToZ( *this ) ;
_type.checkOrthogonalToZ( *this ) ;
}
@@ -382,7 +385,6 @@ namespace DD4hep {
*/
VolCylinderImpl( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D origin ) ;
-
/** 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
*/
@@ -396,13 +398,62 @@ namespace DD4hep {
/** Distance to surface */
virtual double distance(const Vector3D& point ) const ;
- /** Convert the global position to the local position (u,v) on the surface - u runs along the axis of the cylinder, v is r*phi */
+ /** Convert the global position to the local position (u,v) on the surface - v runs along the axis of the cylinder, u is r*phi */
virtual Vector2D globalToLocal( const Vector3D& point) const ;
- /** Convert the local position (u,v) on the surface to the global position - u runs along the axis of the cylinder, v is r*phi*/
+ /** Convert the local position (u,v) on the surface to the global position - v runs along the axis of the cylinder, u is r*phi*/
virtual Vector3D localToGlobal( const Vector2D& point) const ;
} ;
+ //======================================================================================================
+ /** Implementation of conical surface attached to a volume
+ * @author F.Gaede, DESY
+ * @date Nov, 6 2015
+ * @version $Id$
+ */
+ class VolConeImpl : public VolSurfaceBase {
+
+ public:
+
+ /// default c'tor
+ VolConeImpl() : VolSurfaceBase() { }
+
+
+ /** The standard constructor. The origin vector points to the origin of the coordinate system on the cone,
+ * its rho defining the mean radius of the cone (z-component of the origin is ignored !).
+ * The measurement direction v defines the opening angle of the cone,
+ * the normal is chosen to be orthogonal to v. NB: the cone is always parallel to the local z axis.
+ */
+ VolConeImpl( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer,
+ Vector3D v, Vector3D origin ) ;
+
+ /** First direction of measurement U - rotated to point projected onto the cone.
+ * No check is done whether the point actually is on the cone surface
+ */
+ virtual Vector3D u( const Vector3D& point = Vector3D() ) const ;
+
+ /** Second direction of measurement V - rotated to point projected onto the cone.
+ * No check is done whether the point actually is on the cone surface
+ */
+ virtual Vector3D v( const Vector3D& point = Vector3D() ) const ;
+
+ /** The normal direction at the given point, projected onto the cone.
+ * No check is done whether the point actually is on the cone surface
+ */
+ virtual Vector3D normal(const Vector3D& point = Vector3D() ) const ;
+
+ /** Distance to surface */
+ virtual double distance(const Vector3D& point ) const ;
+
+ /** Convert the global position to the local position (u,v) on the surface - v runs along the axis of the cone, u is r*phi */
+ virtual Vector2D globalToLocal( const Vector3D& point) const ;
+
+ /** Convert the local position (u,v) on the surface to the global position - v runs along the axis of the cone, u is r*phi*/
+ virtual Vector3D localToGlobal( const Vector2D& point) const ;
+
+ virtual std::vector< std::pair<DDSurfaces::Vector3D, DDSurfaces::Vector3D> > getLines(unsigned nMax=100) ;
+ } ;
+
//======================================================================================================
@@ -436,6 +487,12 @@ namespace DD4hep {
VolSurface( new VolCylinderImpl( vol, type, thickness_inner , thickness_outer, origin ) ) {}
} ;
+ class VolCone : public VolSurface{
+ public:
+ VolCone( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D v, Vector3D origin ) :
+ VolSurface( new VolConeImpl( vol, type, thickness_inner , thickness_outer, v, origin ) ) {}
+ } ;
+
//======================================================================================================
/** Implementation of Surface class holding a local surface attached to a volume and the DetElement
@@ -539,12 +596,6 @@ namespace DD4hep {
*/
virtual double length_along_v() const ;
- //---------------------------------------------------
- /** Get vertices constraining the surface for drawing ( might not be exact boundaries) -
- * at most nMax points are returned.
- */
- // std::vector< Vector3D > getVertices( unsigned nMax=360 ) ;
-
/** Get lines constraining the surface for drawing ( might not be exact boundaries) -
* at most nMax lines are returned.
*/
@@ -596,10 +647,33 @@ namespace DD4hep {
/// the center of the cylinder
virtual Vector3D center() const ;
-
-
} ;
+ //======================================================================================================
+
+ /** Specialization of Surface for cones - specializes CyliderSurface (for lazyness ...)
+ * @author F.Gaede, DESY
+ * @date May, 10 2014
+ */
+ class ConeSurface : public CylinderSurface, public ICone {
+ public:
+ ConeSurface( Geometry::DetElement det, VolSurface volSurf ) : CylinderSurface( det, volSurf ) { }
+
+ /// the start radius of the cone
+ virtual double radius0() const ;
+
+ /// the end radius of the cone
+ virtual double radius1() const ;
+ /// the start z of the cone
+ virtual double z0() const ;
+
+ /// the end z of the cone
+ virtual double z1() const ;
+
+ /// the center of the cone
+ virtual Vector3D center() const ;
+
+ };
//======================================================================================================
/** std::list of Surfaces that optionally takes ownership.
* @author F.Gaede, DESY
diff --git a/DDRec/src/DetectorSurfaces.cpp b/DDRec/src/DetectorSurfaces.cpp
index 632021ab4..1885133b7 100644
--- a/DDRec/src/DetectorSurfaces.cpp
+++ b/DDRec/src/DetectorSurfaces.cpp
@@ -43,8 +43,17 @@ namespace DD4hep {
VolSurface volSurf = *it ;
- Surface* surf = ( volSurf.type().isCylinder() ? new CylinderSurface( det, volSurf ) : new Surface( det, volSurf ) ) ;
+ Surface* surf = 0 ;
+ if( volSurf.type().isCylinder() )
+ surf = new CylinderSurface( det, volSurf ) ;
+
+ else if( volSurf.type().isCone() )
+ surf = new ConeSurface( det, volSurf ) ;
+
+ else
+ surf = new Surface( det, volSurf ) ;
+
// std::cout << " ------------------------- "
// << " surface: " << *surf << std::endl
// << " ------------------------- " << std::endl ;
diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp
index 65d146bf7..7f22f0a71 100644
--- a/DDRec/src/Surface.cpp
+++ b/DDRec/src/Surface.cpp
@@ -25,6 +25,7 @@ namespace DD4hep {
void VolSurfaceBase::setU(const Vector3D& u_val) { _u = u_val ; }
void VolSurfaceBase::setV(const Vector3D& v_val) { _v = v_val ; }
void VolSurfaceBase::setNormal(const Vector3D& n) { _n = n ; }
+ void VolSurfaceBase::setOrigin(const Vector3D& o) { _o = o ; }
long64 VolSurfaceBase::id() const { return _id ; }
@@ -268,6 +269,7 @@ namespace DD4hep {
_type.setProperty( SurfaceType::Plane , false ) ;
_type.setProperty( SurfaceType::Cylinder , true ) ;
+ _type.setProperty( SurfaceType::Cone , false ) ;
_type.checkParallelToZ( *this ) ;
_type.checkOrthogonalToZ( *this ) ;
}
@@ -316,7 +318,140 @@ namespace DD4hep {
}
//================================================================================================================
+ VolConeImpl::VolConeImpl( Geometry::Volume vol, SurfaceType typ,
+ double thickness_inner ,double thickness_outer, Vector3D v, Vector3D o ) :
+
+ VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , v , Vector3D() , Vector3D() , vol, 0) {
+
+ Vector3D o_rphi( o.x() , o.y() , 0. ) ;
+
+ // sanity check: v and o have to have a common phi
+ double dphi = v.phi() - o_rphi.phi() ;
+ while( dphi < -M_PI ) dphi += 2.*M_PI ;
+ while( dphi > M_PI ) dphi -= 2.*M_PI ;
+
+ if( std::fabs( dphi ) > 1e-6 ){
+ std::stringstream sst ; sst << "VolConeImpl::VolConeImpl() - incompatibel vector v and o given "
+ << v << " - " << o ;
+ throw std::runtime_error( sst.str() ) ;
+ }
+
+ Vector3D n( 1. , v.phi() , ( v.theta() + M_PI/2. ) , Vector3D::spherical ) ;
+ Vector3D u = v.cross( n ) ;
+
+ setU( u ) ;
+ setOrigin( o_rphi ) ;
+ setNormal( n ) ;
+
+ _type.setProperty( SurfaceType::Plane , false ) ;
+ _type.setProperty( SurfaceType::Cylinder, false ) ;
+ _type.setProperty( SurfaceType::Cone , true ) ;
+ _type.setProperty( SurfaceType::ParallelToZ, true ) ;
+ _type.setProperty( SurfaceType::OrthogonalToZ, false ) ;
+ }
+
+ Vector3D VolConeImpl::v(const Vector3D& point ) const {
+ // just take phi from point
+ Vector3D av( 1. , point.phi() , _v.theta() , Vector3D::spherical ) ;
+ return av ;
+ }
+
+ Vector3D VolConeImpl::u(const Vector3D& point ) const {
+ // compute from v X n
+ const Vector3D& av = this->v( point ) ;
+ const Vector3D& n = normal( point ) ;
+ return av.cross( n ) ;
+ }
+
+ Vector3D VolConeImpl::normal(const Vector3D& point ) const {
+ // just take phi from point
+ Vector3D n( 1. , point.phi() , _n.theta() , Vector3D::spherical ) ;
+ return n ;
+ }
+
+ Vector2D VolConeImpl::globalToLocal( const Vector3D& point) const {
+
+ // cone is parallel to z here, so u is r *Phi and v is "along" z
+ double phi = point.phi() - origin().phi() ;
+
+ while( phi < -M_PI ) phi += 2.*M_PI ;
+ while( phi > M_PI ) phi -= 2.*M_PI ;
+
+ return Vector2D( origin().rho()*phi, ( point.z() - origin().z() ) / cos( _v.theta() ) ) ;
+ }
+
+
+ Vector3D VolConeImpl::localToGlobal( const Vector2D& point) const {
+
+ double z = point.v() * cos( _v.theta() ) + origin().z() ;
+
+ double phi = point.u() / origin().rho() + origin().phi() ;
+
+ while( phi < -M_PI ) phi += 2.*M_PI ;
+ while( phi > M_PI ) phi -= 2.*M_PI ;
+
+ return Vector3D( origin().rho() , phi, z , Vector3D::cylindrical ) ;
+ }
+
+
+ /** Distance to surface */
+ double VolConeImpl::distance(const Vector3D& point ) const {
+
+ //fixme: there are probably faster ways to compute this
+ // e.g by using the intercept theorem - tbd. ...
+ const Vector2D& lp = globalToLocal( point ) ;
+ const Vector3D& gp = localToGlobal( lp ) ;
+
+ Vector3D dz = point - gp ;
+
+ return dz * normal( point ) ;
+ }
+
+ /// create outer bounding lines for the given symmetry of the polyhedron
+ std::vector< std::pair<DDSurfaces::Vector3D, DDSurfaces::Vector3D> > VolConeImpl::getLines(unsigned nMax){
+
+ std::vector< std::pair<DDSurfaces::Vector3D, DDSurfaces::Vector3D> > lines ;
+
+ lines.reserve( nMax ) ;
+
+ double theta = v().theta() ;
+ double half_length = 0.5 * length_along_v() * cos( theta ) ;
+
+ Vector3D zv( 0. , 0. , half_length ) ;
+
+ double dr = half_length * tan( theta ) ;
+
+ double r0 = origin().rho() - dr ;
+ double r1 = origin().rho() + dr ;
+
+
+ unsigned n = nMax / 4 ;
+ double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
+
+ for( unsigned i = 0 ; i < n ; ++i ) {
+
+ Vector3D r0v0( r0*sin( i *dPhi ) , r0*cos( i *dPhi ) , 0. ) ;
+ Vector3D r0v1( r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi ) , 0. ) ;
+
+ Vector3D r1v0( r1*sin( i *dPhi ) , r1*cos( i *dPhi ) , 0. ) ;
+ Vector3D r1v1( r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi ) , 0. ) ;
+
+ Vector3D pl0 = zv + r1v0 ;
+ Vector3D pl1 = zv + r1v1 ;
+ Vector3D pl2 = -zv + r0v1 ;
+ Vector3D pl3 = -zv + r0v0 ;
+
+ lines.push_back( std::make_pair( pl0, pl1 ) ) ;
+ lines.push_back( std::make_pair( pl1, pl2 ) ) ;
+ lines.push_back( std::make_pair( pl2, pl3 ) ) ;
+ lines.push_back( std::make_pair( pl3, pl0 ) ) ;
+ }
+ return lines;
+ }
+
+ //================================================================================================================
+
VolSurfaceList::~VolSurfaceList(){
// // delete all surfaces attached to this volume
@@ -527,7 +662,6 @@ namespace DD4hep {
Vector3D localPoint( pa ) ;
return _volSurf.distance( localPoint ) ;
- //FG return ( _volSurf.type().isPlane() ? VolPlane(_volSurf).distance( localPoint ) : VolCylinder(_volSurf).distance( localPoint ) ) ;
}
bool Surface::insideBounds(const Vector3D& point, double epsilon) const {
@@ -537,8 +671,6 @@ namespace DD4hep {
Vector3D localPoint( pa ) ;
return _volSurf.insideBounds( localPoint , epsilon) ;
-
- //FG return ( _volSurf.type().isPlane() ? VolPlane(_volSurf).insideBounds( localPoint, epsilon ) : VolCylinder(_volSurf).insideBounds( localPoint , epsilon) ) ;
}
void Surface::initialize() {
@@ -554,11 +686,6 @@ namespace DD4hep {
throw std::runtime_error( sst.str() ) ;
}
- // std::cout << " **** Surface::initialize() # placements for surface = " << pVList.size()
- // << " worldTransform : "
- // << std::endl ;
-
-
//=========== compute and cache world transform for surface ==========
const TGeoHMatrix& wm = _det.worldTransformation() ;
@@ -674,8 +801,8 @@ namespace DD4hep {
for( unsigned i=0;i<n;++i){
DDSurfaces::Vector3D av,bv;
- _wtM->LocalToMasterVect( local_lines[i].first , av.array() ) ;
- _wtM->LocalToMasterVect( local_lines[i].second , bv.array() ) ;
+ _wtM->LocalToMaster( local_lines[i].first , av.array() ) ;
+ _wtM->LocalToMaster( local_lines[i].second , bv.array() ) ;
lines.push_back( std::make_pair( av, bv ) );
}
@@ -1075,6 +1202,44 @@ namespace DD4hep {
//================================================================================================================
+
+ double ConeSurface::radius0() const {
+
+ double theta = _volSurf.v().theta() ;
+ double l = length_along_v() * cos( theta ) ;
+
+ return origin().rho() - 0.5 * l * tan( theta ) ;
+ }
+
+ double ConeSurface::radius1() const {
+
+ double theta = _volSurf.v().theta() ;
+ double l = length_along_v() * cos( theta ) ;
+
+ return origin().rho() + 0.5 * l * tan( theta ) ;
+ }
+
+ double ConeSurface::z0() const {
+
+ double theta = _volSurf.v().theta() ;
+ double l = length_along_v() * cos( theta ) ;
+
+ return origin().z() - 0.5 * l ;
+ }
+
+ double ConeSurface::z1() const {
+
+ double theta = _volSurf.v().theta() ;
+ double l = length_along_v() * cos( theta ) ;
+
+ return origin().z() + 0.5 * l ;
+ }
+
+ Vector3D ConeSurface::center() const { return volumeOrigin() ; }
+
+
+ //================================================================================================================
+
} // namespace
} // namespace
diff --git a/DDSurfaces/include/DDSurfaces/ISurface.h b/DDSurfaces/include/DDSurfaces/ISurface.h
index 5c37a7cc2..38d7bfdc5 100644
--- a/DDSurfaces/include/DDSurfaces/ISurface.h
+++ b/DDSurfaces/include/DDSurfaces/ISurface.h
@@ -98,6 +98,23 @@ namespace DDSurfaces {
virtual double radius() const=0 ;
virtual Vector3D center() const=0 ;
};
+ //==============================================================================================
+ /** Minimal interface to provide acces to radii of conical surfaces.
+ * @author F. Gaede, DESY
+ * @version $Id$
+ * @date Nov 6 2015
+ */
+ class ICone {
+
+ public:
+ /// Destructor
+ virtual ~ICone() {}
+ virtual double radius0() const=0 ;
+ virtual double radius1() const=0 ;
+ virtual double z0() const=0 ;
+ virtual double z1() const=0 ;
+ virtual Vector3D center() const=0 ;
+ };
//==============================================================================================
@@ -120,7 +137,8 @@ namespace DDSurfaces {
ParallelToZ,
OrthogonalToZ,
Invisible,
- Measurement1D
+ Measurement1D,
+ Cone
} ;
///default c'tor
@@ -176,6 +194,9 @@ namespace DDSurfaces {
/// true if this a cylindrical surface
bool isCylinder() const { return _bits[ SurfaceType::Cylinder ] ; }
+ /// true if this a conical surface
+ bool isCone() const { return _bits[ SurfaceType::Cone ] ; }
+
/// true if surface is parallel to Z
bool isParallelToZ() const { return _bits[ SurfaceType::ParallelToZ ] ; }
@@ -188,6 +209,9 @@ namespace DDSurfaces {
/// true if this is a cylinder parallel to Z
bool isZCylinder() const { return ( _bits[ SurfaceType::Cylinder ] && _bits[ SurfaceType::ParallelToZ ] ) ; }
+ /// true if this is a cone parallel to Z
+ bool isZCone() const { return ( _bits[ SurfaceType::Cone ] && _bits[ SurfaceType::ParallelToZ ] ) ; }
+
/// true if this is a plane parallel to Z
bool isZPlane() const { return ( _bits[ SurfaceType::Plane ] && _bits[ SurfaceType::ParallelToZ ] ) ;
}
@@ -246,9 +270,17 @@ namespace DDSurfaces {
/// dump SurfaceType operator
inline std::ostream& operator<<( std::ostream& os , const SurfaceType& t ) {
- os << "sensitive[" << t.isSensitive() << "] helper[" << t.isHelper() << "] plane[" << t.isPlane() << "] cylinder[" << t.isCylinder()
- << "] parallelToZ[" << t.isParallelToZ() << "] orthogonalToZ[" << t. isOrthogonalToZ() << "] zCylinder[" << t.isZCylinder()
- << "] zPlane[" << t.isZPlane() << "] zDisk[" << t.isZDisk() << "]" ;
+ os << "sensitive[" << t.isSensitive()
+ << "] helper[" << t.isHelper()
+ << "] plane[" << t.isPlane()
+ << "] cylinder[" << t.isCylinder()
+ << "] cone[" << t.isCone()
+ << "] parallelToZ[" << t.isParallelToZ()
+ << "] orthogonalToZ[" << t.isOrthogonalToZ()
+ << "] zCylinder[" << t.isZCylinder()
+ << "] zCone[" << t.isZCone()
+ << "] zPlane[" << t.isZPlane()
+ << "] zDisk[" << t.isZDisk() << "]" ;
return os ;
}
@@ -264,10 +296,13 @@ namespace DDSurfaces {
<< " outerMaterial : " << s.outerMaterial() << " thickness: " << s.outerThickness() << std::endl ;
const ICylinder* cyl = dynamic_cast< const ICylinder* > ( &s ) ;
-
if( cyl )
os << " cylinder radius : " << cyl->radius() << std::endl ;
+ const ICone* cone = dynamic_cast< const ICone* > ( &s ) ;
+ if( cone )
+ os << " cone radius0: " << cone->radius0() << " cone radius1: " << cone->radius1() << std::endl ;
+
return os ;
}
diff --git a/DDSurfaces/include/DDSurfaces/Vector3D.h b/DDSurfaces/include/DDSurfaces/Vector3D.h
index 819de8b88..23a826e30 100644
--- a/DDSurfaces/include/DDSurfaces/Vector3D.h
+++ b/DDSurfaces/include/DDSurfaces/Vector3D.h
@@ -344,15 +344,14 @@ namespace DDSurfaces {
/** Output operator */
inline std::ostream & operator << (std::ostream & os, const Vector3D &v) {
- os << "( " << v[0] << ", " << v[1] << ", " << v[2] << " )" ;
-
- // os << " ( " << v[0]
- // << ", " << v[1]
- // << ", " << v[2]
- // << " ) - [ phi: " << v.phi()
- // << " , rho: " << v.rho() << " ] "
- // << " [ theta: " << v.theta()
- // << " , r: " << v.r() << " ] " ;
+ // os << "( " << v[0] << ", " << v[1] << ", " << v[2] << " )" ;
+ os << " ( " << v[0]
+ << ", " << v[1]
+ << ", " << v[2]
+ << " ) - [ phi: " << v.phi()
+ << " , rho: " << v.rho() << " ] "
+ << " [ theta: " << v.theta()
+ << " , r: " << v.r() << " ] " ;
return os ;
}
--
GitLab