Surface.cpp

#include "DDRec/Surface.h"
#include "DD4hep/objects/DetectorInterna.h"
#include "DD4hep/Memory.h"

#include "DDRec/MaterialManager.h"

#include <cmath>
#include <memory>
#include <exception>

#include "TGeoMatrix.h"
#include "TGeoShape.h"
#include "TRotation.h"
//TGeoTrd1 is apparently not included by defautl
#include "TGeoTrd1.h"

namespace DD4hep {
  namespace DDRec {
 
    using namespace Geometry ;


      //======================================================================================================
  
    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 ; } 

    const SurfaceType& VolSurfaceBase::type() const { return _type ; }
    Vector3D VolSurfaceBase::u( const Vector3D& point ) const {  point.x() ; return _u ; }
    Vector3D VolSurfaceBase::v(const Vector3D& point  ) const { point.x() ;  return _v ; }
    Vector3D VolSurfaceBase::normal(const Vector3D& point ) const {  point.x() ; return _n ; }
    const Vector3D& VolSurfaceBase::origin() const { return _o ;}

    Vector2D VolSurfaceBase::globalToLocal( const Vector3D& point) const {

      Vector3D p = point - origin() ;

      // create new orthogonal unit vectors
      // FIXME: these vectors should be cached really ... 

      double uv = u() * v() ;
      Vector3D uprime = ( u() - uv * v() ).unit() ; 
      Vector3D vprime = ( v() - uv * u() ).unit() ; 
      double uup = u() * uprime ;
      double vvp = v() * vprime ;
      
      return  Vector2D(   p*uprime / uup ,  p*vprime / vvp ) ;
    }
    
    Vector3D VolSurfaceBase::localToGlobal( const Vector2D& point) const {

      Vector3D g = origin() + point[0] * u() + point[1] * v() ;

      return g ;
    }

    const IMaterial&  VolSurfaceBase::innerMaterial() const{  return  _innerMat ;  }
    const IMaterial&  VolSurfaceBase::outerMaterial() const { return  _outerMat  ; }
    double VolSurfaceBase::innerThickness() const { return _th_i ; }
    double VolSurfaceBase::outerThickness() const { return _th_o ; }
    

    double VolSurfaceBase::length_along_u() const {
      
      const DDSurfaces::Vector3D& o = this->origin() ;
      const DDSurfaces::Vector3D& u_val = this->u( o ) ;      
      DDSurfaces::Vector3D  um = -1. * u_val ;
      
      double dist_p = 0. ;
      double dist_m = 0. ;
      

      // std::cout << " VolSurfaceBase::length_along_u() : o =  " << o << " u = " <<    this->u( o ) 
      // 		<< " -u = " << um << std::endl ;


      if( volume()->GetShape()->Contains( o.const_array() ) ){

	dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
							      const_cast<double*> ( u_val.const_array() ) ) ;
	dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
							      const_cast<double*> ( um.array()      ) ) ;
	

	// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
	// 	  << " dist_p " <<    dist_p
	// 	  << " dist_m " <<    dist_m
	// 	  << std::endl ;
	

      } else{

	dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
							       const_cast<double*> ( u_val.const_array() ) ) ;
	dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
							       const_cast<double*> ( um.array()      ) ) ;

	dist_p *= 1.0001 ;
	dist_m *= 1.0001 ;

	// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
	// 	  << " dist_p " <<    dist_p
	// 	  << " dist_m " <<    dist_m
	// 	  << std::endl ;

	DDSurfaces::Vector3D o_1 = this->origin() + dist_p * u_val ;
	DDSurfaces::Vector3D o_2 = this->origin() + dist_m * um ;

	dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) , 
							const_cast<double*> ( u_val.const_array() ) ) ;

	dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) , 
							const_cast<double*> ( um.array()      ) ) ;

	// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
	// 	  << " dist_p " <<    dist_p
	// 	  << " dist_m " <<    dist_m
	// 	  << std::endl ;
      }
	
      return dist_p + dist_m ;


    }
    
    double VolSurfaceBase::length_along_v() const {

      const DDSurfaces::Vector3D& o = this->origin() ;
      const DDSurfaces::Vector3D& v_val = this->v( o ) ;      
      DDSurfaces::Vector3D  vm = -1. * v_val ;
      
      double dist_p = 0. ;
      double dist_m = 0. ;
      

      // std::cout << " VolSurfaceBase::length_along_u() : o =  " << o << " u = " <<    this->u( o ) 
      // 		<< " -u = " << vm << std::endl ;


      if( volume()->GetShape()->Contains( o.const_array() ) ){

	dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
							      const_cast<double*> ( v_val.const_array() ) ) ;
	dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
							      const_cast<double*> ( vm.array()      ) ) ;
	

	// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
	// 	  << " dist_p " <<    dist_p
	// 	  << " dist_m " <<    dist_m
	// 	  << std::endl ;
	

      } else{

	dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
							       const_cast<double*> ( v_val.const_array() ) ) ;
	dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
							       const_cast<double*> ( vm.array()      ) ) ;

	dist_p *= 1.0001 ;
	dist_m *= 1.0001 ;

	// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
	// 	  << " dist_p " <<    dist_p
	// 	  << " dist_m " <<    dist_m
	// 	  << std::endl ;

	DDSurfaces::Vector3D o_1 = this->origin() + dist_p * v_val ;
	DDSurfaces::Vector3D o_2 = this->origin() + dist_m * vm ;

	dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) , 
							const_cast<double*> ( v_val.const_array() ) ) ;

	dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) , 
							const_cast<double*> ( vm.array()      ) ) ;

	// std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
	// 	  << " dist_p " <<    dist_p
	// 	  << " dist_m " <<    dist_m
	// 	  << std::endl ;
      }
	
      return dist_p + dist_m ;

    }
    

    double VolSurfaceBase::distance(const Vector3D& point ) const  {  point.x() ; return 1.e99 ; }

    /// Checks if the given point lies within the surface
    bool VolSurfaceBase::insideBounds(const Vector3D& point, double epsilon) const {
      
#if 0
      double dist = std::abs ( distance( point ) ) ;
      
      bool inShape = volume()->GetShape()->Contains( point.const_array() ) ;
      
      std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << dist 
                << " origin = " << origin() << " normal = " << normal() 
                << " p * n = " << point * normal() 
                << " isInShape : " << inShape << std::endl ;
	
      return dist < epsilon && inShape ;
#else
	
      //fixme: older versions of ROOT (~<5.34.10 ) take a non const pointer as argument - therefore use a const cast here for the time being ...
      return ( std::abs ( distance( point ) ) < epsilon )  &&  volume()->GetShape()->Contains( const_cast<double*> (point.const_array() )  ) ; 
#endif
 
    }


    std::vector< std::pair<Vector3D, Vector3D> > VolSurfaceBase::getLines(unsigned ) {
      // dummy implementation returning empty set
      std::vector< std::pair<Vector3D, Vector3D> >  lines ;
      return lines ;
    }

    //===================================================================
    // simple wrapper methods forwarding the call to the implementation object

    long64 VolSurface::id() const  { return _surf->id() ; } 
    const SurfaceType& VolSurface::type() const { return _surf->type() ; }
    Vector3D VolSurface::u( const Vector3D& point ) const {  return _surf->u(point) ; }
    Vector3D VolSurface::v(const Vector3D& point  ) const {  return _surf->v(point) ; }
    Vector3D VolSurface::normal(const Vector3D& point ) const {  return _surf->normal(point) ; }
    const Vector3D& VolSurface::origin() const { return _surf->origin() ;}
    Vector2D VolSurface::globalToLocal( const Vector3D& point) const { return _surf->globalToLocal( point ) ; }
    Vector3D VolSurface::localToGlobal( const Vector2D& point) const { return _surf->localToGlobal( point) ; }
    const IMaterial&  VolSurface::innerMaterial() const{ return _surf->innerMaterial() ; }
    const IMaterial&  VolSurface::outerMaterial() const  { return _surf->outerMaterial()  ; }
    double VolSurface::innerThickness() const { return _surf->innerThickness() ; }
    double VolSurface::outerThickness() const { return _surf->outerThickness() ; }
    double VolSurface::length_along_u() const { return _surf->length_along_u() ; }
    double VolSurface::length_along_v() const { return _surf->length_along_v() ; }
    double VolSurface::distance(const Vector3D& point ) const  {  return _surf->distance( point ) ; }
    bool VolSurface::insideBounds(const Vector3D& point, double epsilon) const {
      return _surf->insideBounds( point, epsilon ) ; 
    }
    std::vector< std::pair<Vector3D, Vector3D> > VolSurface::getLines(unsigned nMax) {
      return _surf->getLines(nMax) ;
    }

    //===================================================================

    /** Distance to planar surface */
    double VolPlaneImpl::distance(const Vector3D& point ) const {
      return ( point - origin() ) *  normal()  ;
    }
   //======================================================================================================

    VolCylinderImpl::VolCylinderImpl( Geometry::Volume vol, SurfaceType typ, 
				      double thickness_inner ,double thickness_outer,  Vector3D o ) :

      VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , Vector3D() , Vector3D() , o , vol, 0) {
      Vector3D v_val( 0., 0., 1. ) ;
      Vector3D o_rphi( o.x() , o.y() , 0. ) ;
      Vector3D n =  o_rphi.unit() ; 
      Vector3D u_val = v_val.cross( n ) ;

      setU( u_val ) ;
      setV( v_val ) ;
      setNormal( n ) ;

      _type.setProperty( SurfaceType::Plane    , false ) ;
      _type.setProperty( SurfaceType::Cylinder , true ) ;
      _type.setProperty( SurfaceType::Cone     , false ) ;
      _type.checkParallelToZ( *this ) ;
      _type.checkOrthogonalToZ( *this ) ;
    }      

    Vector3D VolCylinderImpl::u(const Vector3D& point ) const {  

      Vector3D n( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;

      return v().cross( n ) ; 
    }
    
    Vector3D VolCylinderImpl::normal(const Vector3D& point ) const {  

      // normal is just given by phi of the point 
      return Vector3D( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
    }

    Vector2D VolCylinderImpl::globalToLocal( const Vector3D& point) const {
      
      // cylinder is parallel to v here so u is Z and v is r *Phi
      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() ) ;
    }
    
    
    Vector3D VolCylinderImpl::localToGlobal( const Vector2D& point) const {
      
      double z = point.v() + 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 VolCylinderImpl::distance(const Vector3D& point ) const {
      
      return point.rho() - origin().rho()  ;
    }
    
    //================================================================================================================
    VolConeImpl::VolConeImpl( Geometry::Volume vol, SurfaceType typ, 
                              double thickness_inner ,double thickness_outer, Vector3D v_val,  Vector3D o_val ) :
      
      VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , v_val ,  Vector3D() , Vector3D() , vol, 0) {

      Vector3D o_rphi( o_val.x() , o_val.y() , 0. ) ;

      // sanity check: v and o have to have a common phi
      double dphi = v_val.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_val << " - " << o_val ;
        throw std::runtime_error( sst.str() ) ;
      }
      
      Vector3D n( 1. , v_val.phi() , ( v_val.theta() + M_PI/2. ) , Vector3D::spherical ) ;
      Vector3D u_val = v_val.cross( n ) ;

      setU( u_val ) ;
      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
      // for( VolSurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
      //   i->clear() ;
      // }
    }
    //=======================================================

    SurfaceList::~SurfaceList(){
      
      if( _isOwner ) {
        // delete all surfaces attached to this volume
        for( SurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
          delete (*i) ;
        }
      }
    }

    //================================================================================================================

    VolSurfaceList* volSurfaceList( DetElement& det ) {

      
      VolSurfaceList* list = 0 ;

      try {

        list = det.extension< VolSurfaceList >() ;

      } catch(const std::runtime_error& e){ 
	
        list = det.addExtension<VolSurfaceList >(  new VolSurfaceList ) ; 
      }

      return list ;
    }


    //======================================================================================================================

    bool findVolume( PlacedVolume pv,  Volume theVol, std::list< PlacedVolume >& volList ) {
      

      volList.push_back( pv ) ;
      
      //   unsigned count = volList.size() ;
      //   for(unsigned i=0 ; i < count ; ++i) {
      //   	std::cout << " **" ;
      //   }
      //   std::cout << " searching for volume: " << theVol.name() << " " << std::hex << theVol.ptr() << "  <-> pv.volume : "  << pv.name() << " " <<  pv.volume().ptr() 
      //    		<< " pv.volume().ptr() == theVol.ptr() " <<  (pv.volume().ptr() == theVol.ptr() )
      //    		<< std::endl ;
      

      if( pv.volume().ptr() == theVol.ptr() ) { 
	
        return true ;
	
      } else {
	
        //--------------------------------

        const TGeoNode* node = pv.ptr();
	
        if ( !node ) {
	  
          //	  std::cout <<  " *** findVolume: Invalid  placement:  - node pointer Null for volume:  " << pv.name() << std::endl ;

          throw std::runtime_error("*** findVolume: Invalid  placement:  - node pointer Null ! " + std::string( pv.name()  ) );
        }
        //	Volume vol = pv.volume();
	
        //	std::cout << "              ndau = " << node->GetNdaughters() << std::endl ;

        for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
	  
          TGeoNode* daughter = node->GetDaughter(idau);
          PlacedVolume placement( daughter );
	  
          if ( !placement.data() ) {
            throw std::runtime_error("*** findVolume: Invalid not instrumented placement:"+std::string(daughter->GetName())
                                     +" [Internal error -- bad detector constructor]");
          }
	  
          PlacedVolume pv_dau = Ref_t(daughter); // why use a Ref_t here  ???
	  
          if( findVolume(  pv_dau , theVol , volList ) ) {
	    
            //	    std::cout << "  ----- found in daughter volume !!!  " << std::hex << pv_dau.volume().ptr() << std::endl ;

            return true ;
          } 
        }

        //  ------- not found:

        volList.pop_back() ;

        return false ;
        //--------------------------------

      }
    } 

    //======================================================================================================================

    Surface::Surface( Geometry::DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ), 
                                                                       _wtM(0) , _id( 0) , _type( _volSurf.type() )  {

      initialize() ;
    }      
 
    long64 Surface::id() const { return _id ; }

    const SurfaceType& Surface::type() const { return _type ; }

    Vector3D Surface::u( const Vector3D& point ) const { point.x() ; return _u ; }
    Vector3D Surface::v(const Vector3D& point ) const {  point.x() ; return _v ; }
    Vector3D Surface::normal(const Vector3D& point ) const {  point.x() ; return _n ; }
    const Vector3D& Surface::origin() const { return _o ;}
    double Surface::innerThickness() const { return _volSurf.innerThickness() ; }
    double Surface::outerThickness() const { return _volSurf.outerThickness() ; }
    double Surface::length_along_u() const { return _volSurf.length_along_u() ; }
    double Surface::length_along_v() const { return _volSurf.length_along_v() ; }

    /** Thickness of outer material */
    
    const IMaterial& Surface::innerMaterial() const {
      
      const IMaterial& mat = _volSurf.innerMaterial() ;
      
      if( ! ( mat.Z() > 0 ) ) {
	
        MaterialManager matMgr ;
        
	Vector3D p = _o - innerThickness() * _n  ;

        const MaterialVec& materials = matMgr.materialsBetween( _o , p  ) ;

	_volSurf.setInnerMaterial(  materials.size() > 1  ? 
				    matMgr.createAveragedMaterial( materials ) :
				    materials[0].first  )  ;
      }
      return  mat ;
    }

    const IMaterial& Surface::outerMaterial() const {
      
      const IMaterial& mat = _volSurf.outerMaterial() ;
      
      if( ! ( mat.Z() > 0 ) ) {
	
        MaterialManager matMgr ;
        
	Vector3D p = _o - outerThickness() * _n  ;

        const MaterialVec& materials = matMgr.materialsBetween( _o , p  ) ;

	_volSurf.setOuterMaterial(  materials.size() > 1  ? 
				    matMgr.createAveragedMaterial( materials ) :
				    materials[0].first  )  ;
      }
      return  mat ;
    }
      

    Vector2D Surface::globalToLocal( const Vector3D& point) const {

      Vector3D p = point - origin() ;

      // create new orthogonal unit vectors
      // FIXME: these vectors should be cached really ... 

      double uv = u() * v() ;
      Vector3D uprime = ( u() - uv * v() ).unit() ; 
      Vector3D vprime = ( v() - uv * u() ).unit() ; 
      double uup = u() * uprime ;
      double vvp = v() * vprime ;
      
      return  Vector2D(   p*uprime / uup ,  p*vprime / vvp ) ;
    }
    
    
    Vector3D Surface::localToGlobal( const Vector2D& point) const {

      Vector3D g = origin() + point[0] * u() + point[1] * v() ;
      return g ;
    }


    Vector3D Surface::volumeOrigin() const { 

      double o_array[3] ;

      _wtM->LocalToMaster    ( Vector3D() , o_array ) ;
      
      Vector3D o(o_array) ;
     
      return o ;
    }


    double Surface::distance(const Vector3D& point ) const {

      double pa[3] ;
      _wtM->MasterToLocal( point , pa ) ;
      Vector3D localPoint( pa ) ;
      
      return _volSurf.distance( localPoint ) ;
    }
      
    bool Surface::insideBounds(const Vector3D& point, double epsilon) const {

      double pa[3] ;
      _wtM->MasterToLocal( point , pa ) ;
      Vector3D localPoint( pa ) ;
      
      return _volSurf.insideBounds( localPoint , epsilon) ;
    }

    void Surface::initialize() {
      
      // first we need to find the right volume for the local surface in the DetElement's volumes
      std::list< PlacedVolume > pVList ;
      PlacedVolume pv = _det.placement() ;
      Volume theVol = _volSurf.volume() ;
      
      if( ! findVolume(  pv, theVol , pVList ) ){
        theVol = _volSurf.volume() ;
        std::stringstream sst ; sst << " ***** ERROR: Volume " << theVol.name() << " not found for DetElement " << _det.name()  << " with surface "  ;
        throw std::runtime_error( sst.str() ) ;
      } 

      //=========== compute and cache world transform for surface ==========
      
      const TGeoHMatrix& wm = _det.worldTransformation() ;
      
#if 0 // debug
      wm.Print() ;
      for( std::list<PlacedVolume>::iterator it= pVList.begin(), n = pVList.end() ; it != n ; ++it ){
        PlacedVolume pv = *it ;
        TGeoMatrix* m = pv->GetMatrix();
        std::cout << "  +++ matrix for placed volume : " << std::endl ;
        m->Print() ;
      }
#endif

      // need to get the inverse transformation ( see Detector.cpp )
      // std::auto_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm.Inverse() ) ) ;
      // has been fixed now, no need to get the inverse anymore:
      dd4hep_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm ) ) ;

      //---- if the volSurface is not in the DetElement's volume, we need to mutliply the path to the volume to the
      // DetElements world transform
      for( std::list<PlacedVolume>::iterator it = ++( pVList.begin() ) , n = pVList.end() ; it != n ; ++it ){

      	PlacedVolume pvol = *it ;
      	TGeoMatrix* m = pvol->GetMatrix();
      	// std::cout << "  +++ matrix for placed volume : " << std::endl ;
      	// m->Print() ;
        //wtI->MultiplyLeft( m );

      	wtI->Multiply( m );
      }

      //      std::cout << "  +++ new world transform matrix  : " << std::endl ;

#if 0 //fixme: which convention to use here - the correct should be wtI, however it is the inverse of what is stored in DetElement ???
      dd4hep_ptr<TGeoHMatrix> wt( new TGeoHMatrix( wtI->Inverse() ) ) ;
      wt->Print() ;
      // cache the world transform for the surface
      _wtM = wt.release()  ;
#else
      //      wtI->Print() ;
      // cache the world transform for the surface
      _wtM = wtI.release()  ;
#endif


      //  ============ now fill the global surface vectors ==========================

      double ua[3], va[3], na[3], oa[3] ;

      _wtM->LocalToMasterVect( _volSurf.u()      , ua ) ;
      _wtM->LocalToMasterVect( _volSurf.v()      , va ) ;
      _wtM->LocalToMasterVect( _volSurf.normal() , na ) ;
      _wtM->LocalToMaster    ( _volSurf.origin() , oa ) ;

      _u.fill( ua ) ;
      _v.fill( va ) ;
      _n.fill( na ) ;
      _o.fill( oa ) ;

      // std::cout << " --- local and global surface vectors : ------- " << std::endl 
      // 			<< "    u : " << _volSurf.u()       << "  -  " << _u << std::endl 
      // 			<< "    v : " << _volSurf.v()       << "  -  " << _v << std::endl 
      // 			<< "    n : " << _volSurf.normal()  << "  -  " << _n << std::endl 
      // 			<< "    o : " << _volSurf.origin()  << "  -  " << _o << std::endl ;
      

      //  =========== check parallel and orthogonal to Z ===================
      
      _type.checkParallelToZ( *this ) ;

      _type.checkOrthogonalToZ( *this ) ;
    
      
      //======== set the unique surface ID from the DetElement ( and placements below ? )

      // just use the DetElement ID for now ...
      // or the id set by the user to the VolSurface ...
      _id = ( _volSurf.id()==0 ?  _det.volumeID() : _volSurf.id() ) ;

      // typedef PlacedVolume::VolIDs IDV ;
      // DetElement d = _det ;
      // while( d.isValid() &&  d.parent().isValid() ){
      // 	PlacedVolume pv = d.placement() ;
      // 	if( pv.isValid() ){
      // 	  const IDV& idV = pv.volIDs() ; 
      // 	  std::cout	<< " VolIDs : " << d.name() << std::endl ;
      // 	  for( unsigned i=0, n=idV.size() ; i<n ; ++i){
      // 	    std::cout  << "  " << idV[i].first << " - " << idV[i].second << std::endl ;
      // 	  }
      // 	}
      // 	d = d.parent() ;
      // }
     
    }
    //===================================================================================================================
      
    std::vector< std::pair<Vector3D, Vector3D> > Surface::getLines(unsigned nMax) {


      const static double epsilon = 1e-6 ; 

      std::vector< std::pair<Vector3D, Vector3D> > lines ;

      //--------------------------------------------
      // check if there are lines defined in the VolSurface :
      const std::vector< std::pair<Vector3D, Vector3D> >& local_lines = _volSurf.getLines() ;
      
      if( local_lines.size() > 0 ) {
	unsigned n=local_lines.size() ;
	lines.reserve( n ) ;
	
	for( unsigned i=0;i<n;++i){
	  
	  DDSurfaces::Vector3D av,bv;
	  _wtM->LocalToMaster( local_lines[i].first ,  av.array() ) ;
	  _wtM->LocalToMaster( local_lines[i].second , bv.array() ) ;
	  
	  lines.push_back( std::make_pair( av, bv ) );
	}
	
	return lines ;
      }
      //--------------------------------------------


      // get local and global surface vectors
      const DDSurfaces::Vector3D& lu = _volSurf.u() ;
      //      const DDSurfaces::Vector3D& lv = _volSurf.v() ;
      const DDSurfaces::Vector3D& ln = _volSurf.normal() ;
      const DDSurfaces::Vector3D& lo = _volSurf.origin() ;
      
      Volume vol = volume() ; 
      const TGeoShape* shape = vol->GetShape() ;
      
      
      if( type().isPlane() ) {
	
        if( shape->IsA() == TGeoBBox::Class() ) {
	  
          TGeoBBox* box = ( TGeoBBox* ) shape  ;
	  
          DDSurfaces::Vector3D boxDim(  box->GetDX() , box->GetDY() , box->GetDZ() ) ;  
	  
	  
          bool isYZ = std::fabs(  ln.x() - 1.0 ) < epsilon  ; // normal parallel to x
          bool isXZ = std::fabs(  ln.y() - 1.0 ) < epsilon  ; // normal parallel to y
          bool isXY = std::fabs(  ln.z() - 1.0 ) < epsilon  ; // normal parallel to z
	  
	  
          if( isYZ || isXZ || isXY ) {  // plane is parallel to one of the box' sides -> need 4 vertices from box dimensions
	    
            // if isYZ :
            unsigned uidx = 1 ;
            unsigned vidx = 2 ;
	    
            DDSurfaces::Vector3D ubl(  0., 1., 0. ) ; 
            DDSurfaces::Vector3D vbl(  0., 0., 1. ) ;
	    
            if( isXZ ) {
	      
              ubl.fill( 1., 0., 0. ) ;
              vbl.fill( 0., 0., 1. ) ;
              uidx = 0 ;
              vidx = 2 ;
	      
            } else if( isXY ) {
	      
              ubl.fill( 1., 0., 0. ) ;
              vbl.fill( 0., 1., 0. ) ;
              uidx = 0 ;
              vidx = 1 ;
            }
	    
            DDSurfaces::Vector3D ub ;
            DDSurfaces::Vector3D vb ;
            _wtM->LocalToMasterVect( ubl , ub.array() ) ;
            _wtM->LocalToMasterVect( vbl , vb.array() ) ;
	    
            lines.reserve(4) ;
	    
            lines.push_back( std::make_pair( _o + boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb ,  _o - boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb ) ) ;
            lines.push_back( std::make_pair( _o - boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb ,  _o - boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb ) ) ;
            lines.push_back( std::make_pair( _o - boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb ,  _o + boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb ) ) ;
            lines.push_back( std::make_pair( _o + boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb ,  _o + boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb ) ) ;
	    
            return lines ;
          }	    

        } else if( shape->IsA() == TGeoConeSeg::Class() ) {

          TGeoCone* cone = ( TGeoCone* ) shape  ;

          // can only deal with special case of z-disk and origin in center of cone
          if( type().isZDisk() && lo.rho() < epsilon ) {
	    
            double zhalf = cone->GetDZ() ;
            double rmax1 = cone->GetRmax1() ;
            double rmax2 = cone->GetRmax2() ;
            double rmin1 = cone->GetRmin1() ;
            double rmin2 = cone->GetRmin2() ;
	    
            // two circles around origin 
            // get radii at position of plane 
            double r0 =  rmin1 +  ( rmin2 - rmin1 ) / ( 2. * zhalf )   *  ( zhalf + lo.z()  ) ;  
            double r1 =  rmax1 +  ( rmax2 - rmax1 ) / ( 2. * zhalf )   *  ( zhalf + lo.z()  ) ;  

	    
            unsigned n = nMax / 4 ;
            double dPhi = 2.* ROOT::Math::Pi() / double( n ) ; 
	    
            for( unsigned i = 0 ; i < n ; ++i ) {
	      
              Vector3D rv00(  r0*sin(  i   *dPhi ) , r0*cos(  i   *dPhi )  , 0. ) ;
              Vector3D rv01(  r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi )  , 0. ) ;

              Vector3D rv10(  r1*sin(  i   *dPhi ) , r1*cos(  i   *dPhi )  , 0. ) ;
              Vector3D rv11(  r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi )  , 0. ) ;
	      
	     
              Vector3D pl0 =  lo + rv00 ;
              Vector3D pl1 =  lo + rv01 ;

              Vector3D pl2 =  lo + rv10 ;
              Vector3D pl3 =  lo + rv11 ;
	      

              Vector3D pg0,pg1,pg2,pg3 ;
	      
              _wtM->LocalToMaster( pl0, pg0.array() ) ;
              _wtM->LocalToMaster( pl1, pg1.array() ) ;
              _wtM->LocalToMaster( pl2, pg2.array() ) ;
              _wtM->LocalToMaster( pl3, pg3.array() ) ;
	      
              lines.push_back( std::make_pair( pg0, pg1 ) ) ;
              lines.push_back( std::make_pair( pg2, pg3 ) ) ;
            }

            //add some vertical and horizontal lines so that the disc is seen in the rho-z projection

            n = 4 ; dPhi = 2.* ROOT::Math::Pi() / double( n ) ;

            for( unsigned i = 0 ; i < n ; ++i ) {
	      
              Vector3D rv0(  r0*sin( i * dPhi ) , r0*cos(  i * dPhi )  , 0. ) ;
              Vector3D rv1(  r1*sin( i * dPhi ) , r1*cos(  i * dPhi )  , 0. ) ;
	      
              Vector3D pl0 =  lo + rv0 ;
              Vector3D pl1 =  lo + rv1 ;

              Vector3D pg0,pg1 ;
	      
              _wtM->LocalToMaster( pl0, pg0.array() ) ;
              _wtM->LocalToMaster( pl1, pg1.array() ) ;
	      
              lines.push_back( std::make_pair( pg0, pg1 ) ) ;
            }

          }
	  
          return lines ;
        }

        else if(shape->IsA() == TGeoTrap::Class()) {
          TGeoTrap* trapezoid = ( TGeoTrap* ) shape;
          
          double dx1 = trapezoid->GetBl1();
          double dx2 = trapezoid->GetTl1();
          double dz = trapezoid->GetH1();

          //according to the TGeoTrap definition, the lengths are given such that the normal vector of the surface
          //points in the e_z direction.
          DDSurfaces::Vector3D ubl(  1., 0., 0. ) ; 
          DDSurfaces::Vector3D vbl(  0., 1., 0. ) ; 

          //the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
          DDSurfaces::Vector3D ub ;
          DDSurfaces::Vector3D vb ;
          _wtM->LocalToMasterVect( ubl , ub.array() ) ;
          _wtM->LocalToMasterVect( vbl , vb.array() ) ;

          //the trapezoid is drawn as a set of four lines connecting its four corners
          lines.reserve(4) ;
          //_o is vector to the origin
          lines.push_back( std::make_pair( _o + dx1 * ub  - dz * vb ,  _o + dx2 * ub  + dz * vb ) ) ;
          lines.push_back( std::make_pair( _o + dx2 * ub  + dz * vb ,  _o - dx2 * ub  + dz * vb ) ) ;
          lines.push_back( std::make_pair( _o - dx2 * ub  + dz * vb ,  _o - dx1 * ub  - dz * vb) ) ;
          lines.push_back( std::make_pair( _o - dx1 * ub  - dz * vb , _o + dx1 * ub  - dz * vb) ) ;

          return lines;
        }
        //added code by Thorben Quast for simplified set of lines for trapezoids with unequal lengths in x
        else if(shape->IsA() == TGeoTrd1::Class()){
          TGeoTrd1* trapezoid = ( TGeoTrd1* ) shape;
          //all lengths are half length
          double dx1 = trapezoid->GetDx1();
          double dx2 = trapezoid->GetDx2();
          //double dy = trapezoid->GetDy();
          double dz = trapezoid->GetDz();

          //the normal vector is parallel to e_y for all geometry cases in CLIC
          //if that is at some point not the case anymore, then local plane vectors ubl, vbl
          //must be initialized like it is done for the boxes (line 674 following)
          DDSurfaces::Vector3D ubl(  1., 0., 0. ) ; 
          DDSurfaces::Vector3D vbl(  0., 0., 1. ) ; 
          
          //the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
          DDSurfaces::Vector3D ub ;
          DDSurfaces::Vector3D vb ;
          _wtM->LocalToMasterVect( ubl , ub.array() ) ;
          _wtM->LocalToMasterVect( vbl , vb.array() ) ;