From 4318bfba4f18d4814337ab508978c964e33e17bf Mon Sep 17 00:00:00 2001
From: Frank Gaede <frank.gaede@desy.de>
Date: Sun, 6 Apr 2014 16:29:09 +0000
Subject: [PATCH] - restructured surfaces: - modified ISurface to be pure
abstract base (interface) - changed method signatures - added
IMaterial as pure abstract base (interface) - moved Surface from DDCore to
DDRec - classes VolSurface (VolPlane, VolCylinder) implement ISurface
- implement insideBounds using volume - changed Vector3D to be copy of
gear::Vector3D (slightly improved)
---
DDRec/CMakeLists.txt | 1 +
DDRec/include/DDRec/Surface.h | 252 ++++++++++++++++
DDRec/src/Surface.cpp | 123 ++++++++
DDSurfaces/include/DDSurfaces/IMaterial.h | 39 +++
DDSurfaces/include/DDSurfaces/ISurface.h | 236 ++++++++++-----
DDSurfaces/include/DDSurfaces/Vector3D.h | 347 ++++++++++++++++++----
6 files changed, 876 insertions(+), 122 deletions(-)
create mode 100644 DDRec/include/DDRec/Surface.h
create mode 100644 DDRec/src/Surface.cpp
create mode 100644 DDSurfaces/include/DDSurfaces/IMaterial.h
diff --git a/DDRec/CMakeLists.txt b/DDRec/CMakeLists.txt
index f4bd6f5d9..efbad2154 100644
--- a/DDRec/CMakeLists.txt
+++ b/DDRec/CMakeLists.txt
@@ -4,6 +4,7 @@ project(DDRec)
include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include
${ROOT_INCLUDE_DIR} ${CMAKE_SOURCE_DIR}/DDSegmentation/include
+ ${ROOT_INCLUDE_DIR} ${CMAKE_SOURCE_DIR}/DDSurfaces/include
${CMAKE_SOURCE_DIR}/DDCore/include)
file(GLOB sources src/*.cpp)
diff --git a/DDRec/include/DDRec/Surface.h b/DDRec/include/DDRec/Surface.h
new file mode 100644
index 000000000..f05ef7948
--- /dev/null
+++ b/DDRec/include/DDRec/Surface.h
@@ -0,0 +1,252 @@
+#ifndef DD4Surfaces_Surface_H
+#define DD4Surfaces_Surface_H
+
+#include "DD4hep/Objects.h"
+#include "DD4hep/Volumes.h"
+
+#include "DDSurfaces/ISurface.h"
+
+#include <list>
+
+
+namespace DD4hep {
+ namespace DDRec {
+
+ using namespace DDSurfaces ;
+
+
+ /** Wrapper class to Geometry::Material that implements the DDSurfaces::IMaterial interface.
+ *
+ * @author F.Gaede, DESY
+ * @date Apr, 6 2014
+ * @version $Id:$
+ */
+ struct SurfaceMaterial : public virtual Geometry::Material , IMaterial{
+
+ /** Copy c'tor - copies handle */
+ SurfaceMaterial( Geometry::Material m ) : Geometry::Material( m ) {}
+
+ /// averaged proton number
+ virtual double Z() const { return Geometry::Material()->GetMaterial()->GetZ() ; }
+
+ /// averaged atomic number
+ virtual double A() const { return Geometry::Material()->GetMaterial()->GetA() ; }
+
+ /// density - units ?
+ virtual double density() const { return Geometry::Material()->GetMaterial()->GetDensity() ; }
+
+ /// radiation length - units ?
+ virtual double radiationLength() const { return Geometry::Material::radLength() ; }
+
+ /// interaction length - units ?
+ virtual double interactionLength() const { return Geometry::Material::intLength() ; }
+
+ };
+
+ /** Helper class for holding surface data.
+ * @author F.Gaede, DESY
+ * @date Apr, 6 2014
+ * @version $Id:$
+ */
+ struct SurfaceData{
+
+ SurfaceType _type ;
+ Vector3D _u ;
+ Vector3D _v ;
+ Vector3D _n ;
+ Vector3D _o ;
+ double _th_i ;
+ double _th_o ;
+ Geometry::Material _innerMat ;
+ Geometry::Material _outerMat ;
+
+ SurfaceData();
+
+ SurfaceData( SurfaceType type, double thickness_inner ,double thickness_outer,
+ Vector3D u ,Vector3D v ,Vector3D n ,Vector3D o ) ;
+
+ /// Default destructor
+ virtual ~SurfaceData() {}
+
+ /// Copy the object
+
+ void copy(const SurfaceData& c) {
+ _type = c._type ;
+ _u = c._u ;
+ _v = c._v ;
+ _n = c._n ;
+ _o = c._o;
+ _th_i = c._th_i ;
+ _th_o = c._th_o ;
+ _innerMat = c._innerMat ;
+ _outerMat = c._innerMat ;
+ }
+ } ;
+
+
+ /** Implementation of ISurface for a local surface attached to a volume.
+ *
+ * @author F.Gaede, DESY
+ * @date Apr, 6 2014
+ * @version $Id:$
+ */
+ class VolSurface : public Geometry::Handle< SurfaceData > , ISurface {
+
+ protected:
+
+ Geometry::Volume _vol ;
+
+ public:
+
+ virtual ~VolSurface() {}
+
+ VolSurface() { }
+
+ VolSurface( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer,
+ Vector3D u ,Vector3D v ,Vector3D n , Vector3D o = Vector3D(0.,0.,0.) ) ;
+
+ Geometry::Volume volume() const { return _vol ; }
+
+
+ /** properties of the surface encoded in Type.
+ * @see SurfaceType
+ */
+ virtual SurfaceType type() const { return object<SurfaceData>()._type ; }
+
+ //==== geometry ====
+
+ /** First direction of measurement U */
+ virtual const Vector3D& u( const Vector3D& point = Vector3D() ) const { return object<SurfaceData>()._u ; }
+
+ /** Second direction of measurement V */
+ virtual const Vector3D& v(const Vector3D& point = Vector3D() ) const { return object<SurfaceData>()._v ; }
+
+ /// Access to the normal direction at the given point
+ virtual const Vector3D& normal(const Vector3D& point = Vector3D() ) const { return object<SurfaceData>()._n ; }
+
+ /** Get Origin of local coordinate system on surface */
+ virtual const Vector3D& origin() const { return object<SurfaceData>()._o ;}
+
+ /// Access to the material in opposite direction of the normal
+ virtual IMaterial innerMaterial() const{ return SurfaceMaterial( object<SurfaceData>()._innerMat ) ; }
+
+ /// Access to the material in direction of the normal
+ virtual IMaterial outerMaterial() const { return SurfaceMaterial( object<SurfaceData>()._outerMat ) ; }
+
+ /** Thickness of inner material */
+ virtual double innerThickness() const { return object<SurfaceData>()._th_i ; }
+
+ /** Thickness of outer material */
+ virtual double outerThickness() const { return object<SurfaceData>()._th_o ; }
+
+
+ // need default implementations for putting it in list....
+
+ /** Distance to surface */
+ virtual double distance(const Vector3D& point ) const { return 0. ; }
+
+ /// Checks if the given point lies within the surface
+ virtual bool insideBounds(const Vector3D& point, double epsilon=1.e-6) const { return false ; }
+ };
+
+
+ //======================================================================================================
+
+
+ /** std::list of VolSurfaces that takes ownership.
+ * @author F.Gaede, DESY
+ * @date Apr, 6 2014
+ * @version $Id:$
+ */
+ struct VolSurfaceList : std::list< VolSurface > {
+
+ VolSurfaceList() {}
+ // required c'tors for extension mechanism
+ VolSurfaceList(const Geometry::DetElement& d){
+ // anything to do here ?
+ }
+ VolSurfaceList(const VolSurfaceList& c,const Geometry::DetElement& det){
+ // anything to do here ?
+ }
+
+ virtual ~VolSurfaceList(){
+
+ // delete all surfaces attached to this volume
+ for( VolSurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
+ delete (*i).ptr() ;
+ }
+
+ }
+ } ;
+
+ VolSurfaceList* volSurfaceList( Geometry::DetElement& det ) ;
+
+ //======================================================================================================
+
+ /** Implementation of planar surface attached to a volume
+ * @author F.Gaede, DESY
+ * @date Apr, 6 2014
+ * @version $Id:$
+ */
+ class VolPlane : public VolSurface {
+
+ public:
+
+ VolPlane() : VolSurface() { }
+
+ VolPlane( 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 ) {
+
+ type._bits.set( SurfaceType::Cylinder , false ) ;
+ type._bits.set( SurfaceType::Plane , true ) ;
+ }
+
+ /** Distance to surface */
+ virtual double distance(const Vector3D& point ) const ;
+
+ /// Checks if the given point lies within the surface
+ virtual bool insideBounds(const Vector3D& point, double epsilon=1.e-4) const ;
+
+
+ } ;
+
+ //======================================================================================================
+
+ /** Implementation of cylindrical surface attached to a volume
+ * @author F.Gaede, DESY
+ * @date Apr, 6 2014
+ * @version $Id:$
+ */
+ class VolCylinder : public VolSurface {
+
+ public:
+
+ VolCylinder() : VolSurface() { }
+
+ 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 ) {
+
+ type._bits.set( SurfaceType::Plane , false ) ;
+ type._bits.set( SurfaceType::Cylinder , true ) ;
+ }
+
+ /** Distance to surface */
+ virtual double distance(const Vector3D& point ) const ;
+
+ /// Checks if the given point lies within the surface
+ virtual bool insideBounds(const Vector3D& point, double epsilon=1.e-4) const ;
+
+ } ;
+
+ //======================================================================================================
+
+
+
+ } /* namespace */
+} /* namespace */
+
+#endif /* Surface */
diff --git a/DDRec/src/Surface.cpp b/DDRec/src/Surface.cpp
new file mode 100644
index 000000000..b1e56758e
--- /dev/null
+++ b/DDRec/src/Surface.cpp
@@ -0,0 +1,123 @@
+#include "DDRec/Surface.h"
+#include "DD4hep/Detector.h"
+
+
+namespace DD4hep {
+ namespace DDRec {
+
+ using namespace Geometry ;
+
+ SurfaceData::SurfaceData() : _type( SurfaceType() ) ,
+ _u( Vector3D() ) ,
+ _v( Vector3D() ) ,
+ _n( Vector3D() ) ,
+ _o( Vector3D() ) ,
+ _th_i( 0. ),
+ _th_o( 0. ),
+ _innerMat( Material() ),
+ _outerMat( Material() ) {
+ }
+
+
+ SurfaceData::SurfaceData( SurfaceType type , double thickness_inner ,double thickness_outer,
+ Vector3D u ,Vector3D v ,Vector3D n ,Vector3D o ) : _type(type ) ,
+ _u( u ) ,
+ _v( v ) ,
+ _n( n ) ,
+ _o( o ),
+ _th_i( thickness_inner ),
+ _th_o( thickness_outer ),
+ _innerMat( Material() ),
+ _outerMat( Material() ) {
+ }
+
+
+
+ VolSurface::VolSurface( Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer,
+ Vector3D u ,Vector3D v ,Vector3D n ,Vector3D o ) :
+
+ Handle( new SurfaceData( type, thickness_inner ,thickness_outer, u,v,n,o) ) ,
+
+ _vol( vol ) {
+ }
+
+
+ /** Distance to surface */
+ double VolPlane::distance(const Vector3D& point ) const {
+
+ return ( point - origin() ) * normal() ;
+ }
+
+ /// Checks if the given point lies within the surface
+ bool VolPlane::insideBounds(const Vector3D& point, double epsilon) const {
+
+
+#if 1
+ double dist = std::abs ( distance( point ) ) ;
+
+ bool inShape = volume()->GetShape()->Contains( point ) ;
+
+ std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << dist
+ << " origin = " << origin() << " normal = " << normal()
+ << " p * n = " << point * normal()
+ << " isInShape : " << inShape << std::endl ;
+
+ return dist < epsilon && inShape ;
+ #else
+
+ return ( std::abs ( distance() ) < epsilon ) && volume()->GetShape()->Contains( point ) ;
+ #endif
+
+ }
+
+ /** Distance to surface */
+ double VolCylinder::distance(const Vector3D& point ) const {
+
+ return point.rho() - origin().rho() ;
+ }
+
+ /// Checks if the given point lies within the surface
+ bool VolCylinder::insideBounds(const Vector3D& point, double epsilon) const {
+
+#if 1
+ double distR = std::abs( distance( point ) ) ;
+
+ bool inShapeT = volume()->GetShape()->Contains( point ) ;
+
+ std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << distR
+ << " origin = " << origin()
+ << " isInShape : " << inShapeT << std::endl ;
+
+ return distR < epsilon && inShapeT ;
+#else
+
+ return ( std::abs ( distance() ) < epsilon ) && volume()->GetShape()->Contains( point ) ;
+
+#endif
+ }
+
+
+
+ //====================
+
+
+ VolSurfaceList* surfaceList( DetElement& det ) {
+
+
+ VolSurfaceList* list = 0 ;
+
+ try {
+
+ list = det.extension< VolSurfaceList >() ;
+
+ } catch( std::runtime_error e){
+
+ list = det.addExtension<VolSurfaceList >( new VolSurfaceList ) ;
+ }
+
+ return list ;
+ }
+
+
+ } // namespace
+} // namespace
diff --git a/DDSurfaces/include/DDSurfaces/IMaterial.h b/DDSurfaces/include/DDSurfaces/IMaterial.h
new file mode 100644
index 000000000..a0b3095a8
--- /dev/null
+++ b/DDSurfaces/include/DDSurfaces/IMaterial.h
@@ -0,0 +1,39 @@
+#ifndef DDSurfaces_IMaterial_H_
+#define DDSurfaces_IMaterial_H_
+
+namespace DDSurfaces {
+
+ /** Interface for material description for tracking.
+ *
+ * @author F.Gaede, DESY
+ * @version $Id$
+ * @date Apr 6 2014
+ */
+
+ class IMaterial {
+
+ public:
+
+ /// Destructor
+ virtual ~IMaterial() {}
+
+ /// averaged proton number
+ virtual double Z() const =0 ;
+
+ /// averaged atomic number
+ virtual double A() const =0 ;
+
+ /// density - units ?
+ virtual double density() const =0 ;
+
+ /// radiation length - units ?
+ virtual double radiationLength() const =0 ;
+
+ /// interaction length - units ?
+ virtual double interactionLength() const =0 ;
+
+ };
+
+} /* namespace DDSurfaces */
+
+#endif /* DDSurfaces_MATERIAL_H_ */
diff --git a/DDSurfaces/include/DDSurfaces/ISurface.h b/DDSurfaces/include/DDSurfaces/ISurface.h
index 52142fab2..b60f6d0c4 100644
--- a/DDSurfaces/include/DDSurfaces/ISurface.h
+++ b/DDSurfaces/include/DDSurfaces/ISurface.h
@@ -1,79 +1,171 @@
-/*
- * ISurface.h
- *
- * Created on: Mar 7, 2014
- * Author: cgrefe
- */
-
-#ifndef DDSurfaces_ISURFACE_H_
-#define DDSurfaces_ISURFACE_H_
-
-#include "DDSurfaces/Material.h"
-#include "DDSurfaces/MeasurementDirections.h"
+#ifndef DDSurfaces_ISurface_H
+#define DDSurfaces_ISurface_H
+
+#include "DDSurfaces/IMaterial.h"
#include "DDSurfaces/Vector3D.h"
+#include <bitset>
+
namespace DDSurfaces {
+
+ struct SurfaceType ;
+
+ /** Interface for tracking surfaces.
+ * The surface provides access to vectors for u,v,n and the orgigin and
+ * the inner and outer materials with corresponding thicknesses.
+ *
+ * @author C.Grefe, CERN, F. Gaede, DESY
+ * @version $Id: $
+ * @date Mar 7 2014
+ */
+ class ISurface {
+
+ public:
+ /// Destructor
+ virtual ~ISurface() {}
+
+ /// properties of the surface encoded in Type.
+ virtual SurfaceType type() const =0 ;
+
+ /// Checks if the given point lies within the surface
+ 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 ;
+
+ /** Second direction of measurement V */
+ virtual const 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 ;
+
+ /** Get Origin of local coordinate system on surface */
+ virtual const Vector3D& origin() const =0 ;
+
+ /// Access to the material in opposite direction of the normal
+ virtual IMaterial innerMaterial() const =0 ;
+
+ /// Access to the material in direction of the normal
+ virtual IMaterial outerMaterial() const =0 ;
+
+ /** Thickness of inner material */
+ virtual double innerThickness() const =0 ;
+
+ /** Thickness of outer material */
+ virtual double outerThickness() const =0 ;
+
+ /** Distance to surface */
+ virtual double distance(const Vector3D& point ) const =0 ;
+
+ } ;
+
+
+ /** Helper class for describing surface properties.
+ * Usage: SurfaceType type( SurfaceType::Plane, SurfaceType::Sensitive ) ;
+ *
+ * @author F. Gaede, DESY
+ * @version $Id: $
+ * @date Apr 6 2014
+ */
+ struct SurfaceType{
+
+ /// enum for defining the bits used to decode the properties
+ enum{
+ Cylinder = 1,
+ Plane,
+ Sensitive,
+ Helper,
+ ParallelToZ,
+ OrthogonalToZ
+ } ;
+
+ ///default c'tor
+ SurfaceType() : _bits(0) {}
+
+ // c'tor that sets one property
+ SurfaceType( unsigned prop0 ) : _bits(0) {
+ _bits.set( prop0 ) ;
+ }
+
+ // c'tor that sets two properties
+ SurfaceType( unsigned prop0 , unsigned prop1 ) : _bits(0) {
+ _bits.set( prop0 ) ;
+ _bits.set( prop1 ) ;
+ }
+
+ // c'tor that sets three properties
+ SurfaceType( unsigned prop0 , unsigned prop1 , unsigned prop2 ) : _bits(0) {
+ _bits.set( prop0 ) ;
+ _bits.set( prop1 ) ;
+ _bits.set( prop2 ) ;
+ }
+
+ // c'tor that sets four properties
+ SurfaceType( unsigned prop0 , unsigned prop1 , unsigned prop2, unsigned prop3 ) : _bits(0) {
+ _bits.set( prop0 ) ;
+ _bits.set( prop1 ) ;
+ _bits.set( prop2 ) ;
+ _bits.set( prop3 ) ;
+ }
+
+ /// set the given peorperty
+ void setProperty( unsigned prop ) { _bits.set( prop ) ; }
+
+ /// true if surface is sensitive
+ bool isSensitive() const { return _bits[ SurfaceType::Sensitive ] ; }
+
+ /// true if surface is helper surface for navigation
+ bool isHelper() const { return _bits[ SurfaceType::Helper ] ; }
+
+ /// true if this a planar surface
+ bool isPlane() const { return _bits[ SurfaceType::Plane ] ; }
+
+ /// true if this a cylindrical surface
+ bool isCylinder() const { return _bits[ SurfaceType::Cylinder ] ; }
+
+ /// true if surface is parallel to Z
+ bool isParallelToZ() const { return _bits[ SurfaceType::ParallelToZ ] ; }
+
+ /// true if surface is orthogonal to Z
+ bool isOrthogonalToZ() const { return _bits[ SurfaceType::OrthogonalToZ ] ; }
+
+ /// true if this is a cylinder parallel to Z
+ bool isZCylinder() const { return ( _bits[ SurfaceType::Cylinder ] && _bits[ SurfaceType::ParallelToZ ] ) ; }
+
+ /// true if this is a plane parallel to Z
+ bool isZPlane() const { return ( _bits[ SurfaceType::Plane ] && _bits[ SurfaceType::ParallelToZ ] ) ; }
+
+ /// true if this is a plane orthogonal to Z
+ bool isZDisk() const { return ( _bits[ SurfaceType::Plane ] && _bits[ SurfaceType::OrthogonalToZ ] ) ; }
+
+
+ /** True if surface is parallel to Z with accuracy epsilon - result is cached in bit SurfaceType::ParallelToZ */
+ bool checkParallelToZ( const ISurface& surf , double epsilon=1.e-6 ) const {
+
+ double proj = std::abs( surf.normal() * Vector3D(0.,0.,1.) ) ;
+
+ _bits.set( SurfaceType::ParallelToZ , ( std::abs( proj - 1. ) < epsilon ) ) ;
+
+ return _bits[ SurfaceType::ParallelToZ ] ;
+ }
+
+ /** True if surface is orthogonal to Z with accuracy epsilon - result is cached in bit SurfaceType::OrthogonalToZ */
+ bool checkOrthogonalToZ( const ISurface& surf , double epsilon=1.e-6 ) const {
+
+ double proj = std::abs( surf.normal() * Vector3D(0.,0.,1.) ) ;
+
+ _bits.set( SurfaceType::OrthogonalToZ , ( proj < epsilon ) ) ;
+
+ return _bits[ SurfaceType::OrthogonalToZ ] ;
+ }
+
+ mutable std::bitset<32> _bits ;
+ } ;
+
+
+
-/** Generic surface class used for tracking */
-class ISurface {
-public:
- /// Destructor
- virtual ~ISurface() {
- if (m_measurement) {
- delete m_measurement;
- }
- }
-
- /// Checks if the surface has measurement directions
- virtual bool hasMeasurement() const {
- return m_measurement != 0;
- }
-
- /// Checks if the given point lies within the surface
- virtual bool isInsideBoundaries(const Vector3D& point) const = 0;
-
- /// Access to the normal direction at the given point
- virtual Vector3D getNormal(const Vector3D& point) const = 0;
-
- /// Access to the measurement directions at the given point
- virtual MeasurementDirections measurement(const Vector3D& point) const = 0;
-
- /// Access to the material in opposite direction of the normal
- const Material& innerMaterial() const {
- return m_innerMaterial;
- }
-
- /// Access to the material in direction of the normal
- const Material& outerMaterial() const {
- return m_outerMaterial;
- }
-
- /// Sets the material in opposite direction of the normal
- void setInnerMaterial(const Material& material) {
- m_innerMaterial = material;
- }
-
- /// Sets the material in direction of the normal
- void setOuterMaterial(const Material& material) {
- m_outerMaterial = material;
- }
-
- /// Sets the nominal measurement directions
- void setMeasurement(const MeasurementDirections& measurement) {
- m_measurement = new MeasurementDirections(measurement);
- }
-
-protected:
- /// Constructor which can be used by derived classes
- ISurface(const Material& innerMaterial = Material(), const Material& outerMaterial = Material()) :
- m_innerMaterial(innerMaterial), m_outerMaterial(outerMaterial), m_measurement(0) {
- }
-
- Material m_innerMaterial; /// the material in opposite direction of the normal
- Material m_outerMaterial; /// the material in direction of the normal
- MeasurementDirections* m_measurement; /// the nominal measurement directions
-};
} /* namespace DDSurfaces */
-#endif /* DDSurfaces_ISURFACE_H_ */
+#endif /* DDSurfaces_ISurface_H */
diff --git a/DDSurfaces/include/DDSurfaces/Vector3D.h b/DDSurfaces/include/DDSurfaces/Vector3D.h
index 3b710a815..7caac864f 100644
--- a/DDSurfaces/include/DDSurfaces/Vector3D.h
+++ b/DDSurfaces/include/DDSurfaces/Vector3D.h
@@ -1,54 +1,301 @@
-/*
- * Vector3D.h
- *
- * Created on: Mar 7, 2014
- * Author: cgrefe
- */
+#ifndef DDSurfaces_Vector3D_h
+#define DDSurfaces_Vector3D_h 1
+
+#include <math.h>
+#include <iostream>
+#include "assert.h"
-#ifndef DDSurfaces_VECTOR3D_H_
-#define DDSurfaces_VECTOR3D_H_
namespace DDSurfaces {
+
+ /** Simple three dimensional vector providing the components for cartesian,
+ * cylindrical and spherical coordinate systems - internal reperesentation is
+ * cartesian. (copy of original version from gear).
+ *
+ * @author F. Gaede, DESY
+ * @version $Id: $
+ * @date Apr 6 2014
+ */
+
+ class Vector3D{
+
+ public:
+
+ /** Default c'tor - zero vector */
+ Vector3D() : _x(0.0),_y(0.0),_z(0.0) {}
+
+
+ /** Constructor for float array.*/
+ Vector3D(const float* v) : _x(v[0]),_y(v[1]),_z(v[2]) {}
+
+ /** Constructor for double array.*/
+ Vector3D(const double* v) : _x(v[0]),_y(v[1]),_z(v[2]) {}
+
+
+ /** Templated c'tor - allows to have overloaded c'tors for different coordinates */
+ template <class T>
+ Vector3D( double x,double y, double z , T(&)() ) ;
+
+
+ /** Default corrdinate system for initialization is cartesian */
+ Vector3D( double x,double y, double z ) :
+ _x(x),
+ _y(y),
+ _z(z) {
+ }
+
+
+ template <class T>
+ /**Copy c'tor for three vectors from other packages - requires T::x(),T::y(), T::z().
+ */
+ Vector3D( const T& t) :
+ _x( t.x() ) ,
+ _y( t.y() ) ,
+ _z( t.z() ){
+ }
+
+ /** Cartesian x coordinate */
+ inline double x() const { return _x ; }
+
+ /** Cartesian y coordinate */
+ inline double y() const { return _y ; }
+
+ /** Cartesian cartesian z coordinate */
+ inline double z() const { return _z ; }
+
+ /** Assign to cartesian x coordinate */
+ inline double& x() { return _x ; }
+
+ /** Assign to cartesian y coordinate */
+ inline double& y() { return _y ; }
+
+ /** Assign to cartesian z coordinate */
+ inline double& z() { return _z ; }
+
+
+ /** Accessing x,y,z with bracket operator */
+ inline double operator[](int i) const {
+ switch(i) {
+ case 0: return _x ; break ;
+ case 1: return _y ; break ;
+ case 2: return _z ; break ;
+ }
+ return 0.0 ;
+ }
+ /** Accessing x,y,z with bracket operator for assignment */
+ inline double& operator[](int i) {
+ switch(i) {
+ case 0: return _x ; break ;
+ case 1: return _y ; break ;
+ case 2: return _z ; break ;
+ }
+ static double dummy(0.0) ;
+ return dummy ;
+ }
+
+ /** Azimuthal angle - cylindrical and spherical */
+ inline double phi() const {
+
+ return _x == 0.0 && _y == 0.0 ? 0.0 : atan2(_y,_x);
+ }
+
+ /** Transversal component - cylindrical 'r' */
+ inline double rho() const {
+
+ return trans() ;
+ }
+
+ /** Transversal component */
+ inline double trans() const {
+
+ return sqrt( _x*_x + _y*_y ) ;
+ }
+
+ /** Transversal component squared */
+ inline double trans2() const {
+
+ return _x*_x + _y*_y ;
+ }
+
+ /** Spherical r/magnitude */
+ inline double r() const {
+
+ return sqrt( _x*_x + _y*_y + _z*_z ) ;
+ }
+
+
+ /** Spherical r/magnitude, squared */
+ inline double r2() const {
+
+ return _x*_x + _y*_y + _z*_z ;
+ }
+
+ /** Polar angle - spherical */
+ inline double theta() const {
+
+ return _x == 0.0 && _y == 0.0 && _z == 0.0 ? 0.0 : atan2( rho(),_z) ;
+ }
+
+ /** Scalar product */
+ inline double dot( const Vector3D& v) const {
+ return _x * v.x() + _y * v.y() + _z * v.z() ;
+ }
+
+
+ /** Vector product */
+ inline Vector3D cross( const Vector3D& v) const {
+
+ return Vector3D( _y * v.z() - _z * v.y() ,
+ _z * v.x() - _x * v.z() ,
+ _x * v.y() - _y * v.x() ) ;
+ }
+
+ /** Parallel unit vector */
+ inline Vector3D unit() const {
+
+ double n = r() ;
+ return Vector3D( _x / n , _y / n , _z / n ) ;
+ }
+
+
+ /// direct access to data as const double*
+ inline operator const double*() const {
+ return &_x ;
+ }
+
+
+ /** Implicit templated conversion to anything that has a c'tor T(x,y,z)
+ * and accessor functions x(),y(),z(). For safety the result is checked which
+ * causes a small performance penalty.
+ * @see to()
+ *
+ */
+ template <class T>
+ inline operator T() const {
+
+ T t( _x, _y , _z ) ;
+
+ assert( t.x()== _x && t.y()== _y && t.z()== _z ) ;
+
+ return t ;
+
+ // return T( _x, _y, _z ) ;
+ }
+
+
+ /** Explicit, unchecked conversion to anything that has a c'tor T(x,y,z).
+ * Example:
+ * CLHEP::Vector3D clhv = v.to< CLHEP::Vector3D>() ;
+ * @see operator T()
+ */
+ template <class T>
+ inline T to() const { return T( _x, _y, _z ) ; }
+
+
+ protected:
+
+ double _x,_y,_z ;
+
+
+ // helper classes and function to allow
+ // different c'tors selected at compile time
+ public:
+
+ struct Cartesian { } ;
+ struct Cylindrical { } ;
+ struct Spherical { } ;
+
+ static Cartesian cartesian() { return Cartesian() ; }
+ static Cylindrical cylindrical(){ return Cylindrical() ;}
+ static Spherical spherical() { return Spherical() ; }
+
+ } ;
+
+ /** Addition of two vectors */
+ inline Vector3D operator+( const Vector3D& a, const Vector3D& b ) {
+
+ return Vector3D( a.x() + b.x() , a.y() + b.y(), a.z() + b.z() ) ;
+ }
+ /** Subtraction of two vectors */
+ inline Vector3D operator-( const Vector3D& a, const Vector3D& b ) {
+
+ return Vector3D( a.x() - b.x() , a.y() - b.y(), a.z() - b.z() ) ;
+ }
+ /** 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() )
+ return true;
+ else
+ return false;
+ }
+
+ /** Multiplication with scalar */
+ inline Vector3D operator*( double s , const Vector3D& v ) {
+
+ return Vector3D( s * v.x() , s * v.y() , s * v.z() ) ;
+ }
+
+ /// operator for scalar product
+ inline double operator*( const Vector3D& v0, const Vector3D& v1 ){
+ return v0.dot( v1 ) ;
+ }
+
+
+ // template specializations for constructors of different coordinate systems
+
+ /** Cartesian c'tor - example: <br>
+ * Vector3D v( x, y, c , Vector3D::cartesian ) ;
+ */
+ template <>
+ inline Vector3D::Vector3D( double x,double y, double z, Vector3D::Cartesian (&)() ) :
+ _x(x),
+ _y(y),
+ _z(z) {
+ }
+
+ /** Cylindrical c'tor - example: <br>
+ * Vector3D v( rho, phi, z , Vector3D::cylindrical ) ;
+ */
+ template <>
+ inline Vector3D::Vector3D( double rho,double phi, double z, Vector3D::Cylindrical (&)() ) : _z(z) {
+
+ _x = rho * cos( phi ) ;
+ _y = rho * sin( phi ) ;
+ }
+
+
+ /** Spherical c'tor - example: <br>
+ * Vector3D v( r, phi, theta , Vector3D::spherical ) ;
+ */
+ template <>
+ inline Vector3D::Vector3D( double r,double phi, double theta, Vector3D::Spherical (&)() ) {
+ double rst = r * sin( theta ) ;
+ _x = rst * cos( phi ) ;
+ _y = rst * sin( phi ) ;
+ _z = r * cos( theta ) ;
+ }
+
+
+
+ /** Output operator */
+ 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() << " ] " ;
+
+ return os ;
+ }
+
+
+
+} // namespace
-/** Generic vector class */
-class Vector3D {
-public:
- /// Default constructor
- Vector3D(double x = 0., double y = 0., double z = 0.) :
- m_x(x), m_y(y), m_z(z) {
- }
-
- /** Copy constructor from arbitrary foreign vector class
- * Requires x(), y(), z() methods to access Cartesian coordinates.
- */
- template<typename T> Vector3D(const T& vector) :
- m_x(vector.x()), m_y(vector.y()), m_z(vector.z()) {
- }
-
- /// Destructor
- virtual ~Vector3D() {}
-
- /// Access to x
- double x() const {
- return m_x;
- }
-
- /// Access to y
- double y() const {
- return m_y;
- }
-
- /// Access to z
- double z() const {
- return m_z;
- }
-
-protected:
- double m_x;
- double m_y;
- double m_z;
-};
-
-} /* namespace DDSurfaces */
-
-#endif /* VECTOR3D_H_ */
+#endif
--
GitLab