Skip to content
Snippets Groups Projects
Surface.cpp 39 KiB
Newer Older
Frank Gaede's avatar
Frank Gaede committed
#include "DDRec/Surface.h"
#include "DD4hep/objects/DetectorInterna.h"
Markus Frank's avatar
Markus Frank committed
#include "DD4hep/Memory.h"
#include "DDRec/MaterialManager.h"

#include <memory>
#include <exception>

#include "TGeoMatrix.h"
#include "TGeoShape.h"
#include "TRotation.h"
//TGeoTrd1 is apparently not included by defautl
#include "TGeoTrd1.h"
Frank Gaede's avatar
Frank Gaede committed

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 ; }
Frank Gaede's avatar
Frank Gaede committed
    
    long64 VolSurfaceBase::id() const  { return 0 ; } 
Frank Gaede's avatar
Frank Gaede committed

    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 ;}
Frank Gaede's avatar
Frank Gaede committed

    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 {
Frank Gaede's avatar
Frank Gaede committed
      
      const DDSurfaces::Vector3D& o = this->origin() ;
Markus Frank's avatar
Markus Frank committed
      const DDSurfaces::Vector3D& u_val = this->u( o ) ;      
      DDSurfaces::Vector3D  um = -1. * u_val ;
Frank Gaede's avatar
Frank Gaede committed
      
      // 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 ;
      }
	
Frank Gaede's avatar
Frank Gaede committed
      return dist_p + dist_m ;
    double VolSurfaceBase::length_along_v() const {
Frank Gaede's avatar
Frank Gaede committed

      const DDSurfaces::Vector3D& o = this->origin() ;
Markus Frank's avatar
Markus Frank committed
      const DDSurfaces::Vector3D& v_val = this->v( o ) ;      
      DDSurfaces::Vector3D  vm = -1. * v_val ;
Frank Gaede's avatar
Frank Gaede committed
      
      // 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 ;
      }
	
Frank Gaede's avatar
Frank Gaede committed
      return dist_p + dist_m ;
    double VolSurfaceBase::distance(const Vector3D& point ) const  {  point.x() ; return 1.e99 ; }
Frank Gaede's avatar
Frank Gaede committed
    /// Checks if the given point lies within the surface
    bool VolSurfaceBase::insideBounds(const Vector3D& point, double epsilon) const {
Frank Gaede's avatar
Frank Gaede committed
      double dist = std::abs ( distance( point ) ) ;
      
      bool inShape = volume()->GetShape()->Contains( point.const_array() ) ;
Frank Gaede's avatar
Frank Gaede committed
      
      std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << dist 
Markus Frank's avatar
Markus Frank committed
                << " origin = " << origin() << " normal = " << normal() 
                << " p * n = " << point * normal() 
                << " isInShape : " << inShape << std::endl ;
Frank Gaede's avatar
Frank Gaede committed
	
      return dist < epsilon && inShape ;
Frank Gaede's avatar
Frank Gaede committed
#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() )  ) ; 
Frank Gaede's avatar
Frank Gaede committed
#endif
    std::vector< std::pair<Vector3D, Vector3D> > VolSurfaceBase::getLines(unsigned nMax) {
      // dummy implementation returning empty set
      nMax ;
      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) ;
    }
Frank Gaede's avatar
Frank Gaede committed

    //===================================================================
Frank Gaede's avatar
Frank Gaede committed

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

    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. ) ;
Frank Gaede's avatar
Frank Gaede committed
      Vector3D n =  o_rphi.unit() ; 
Markus Frank's avatar
Markus Frank committed
      Vector3D u_val = v_val.cross( n ) ;
Frank Gaede's avatar
Frank Gaede committed

Markus Frank's avatar
Markus Frank committed
      setU( u_val ) ;
      setV( v_val ) ;
Frank Gaede's avatar
Frank Gaede committed
      setNormal( n ) ;

      _type.setProperty( SurfaceType::Plane    , false ) ;
      _type.setProperty( SurfaceType::Cylinder , true ) ;
      _type.checkParallelToZ( *this ) ;
      _type.checkOrthogonalToZ( *this ) ;
Frank Gaede's avatar
Frank Gaede committed
    }      
Frank Gaede's avatar
Frank Gaede committed

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

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

Frank Gaede's avatar
Frank Gaede committed
      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 ) ;
Frank Gaede's avatar
Frank Gaede committed
    }
    Vector2D VolCylinderImpl::globalToLocal( const Vector3D& point) const {
Frank Gaede's avatar
Frank Gaede committed
      // 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 {
      
Frank Gaede's avatar
Frank Gaede committed
      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 )    ;
    }


Frank Gaede's avatar
Frank Gaede committed
    /** Distance to surface */
    double VolCylinderImpl::distance(const Vector3D& point ) const {
      
Frank Gaede's avatar
Frank Gaede committed
      return point.rho() - origin().rho()  ;
    }
    
    //================================================================================================================
Frank Gaede's avatar
Frank Gaede committed
    VolSurfaceList::~VolSurfaceList(){
      // // delete all surfaces attached to this volume
      // for( VolSurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
      //   i->clear() ;
      // }
Frank Gaede's avatar
Frank Gaede committed
    }
    //=======================================================

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

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

    VolSurfaceList* volSurfaceList( DetElement& det ) {
Frank Gaede's avatar
Frank Gaede committed

      
      VolSurfaceList* list = 0 ;

      try {

Markus Frank's avatar
Markus Frank committed
        list = det.extension< VolSurfaceList >() ;
      } catch(const std::runtime_error& e){ 
Markus Frank's avatar
Markus Frank committed
        list = det.addExtension<VolSurfaceList >(  new VolSurfaceList ) ; 
    //======================================================================================================================

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

      volList.push_back( pv ) ;
      
Frank Gaede's avatar
Frank Gaede committed
      //   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() ) { 
	
Markus Frank's avatar
Markus Frank committed
        return true ;
Markus Frank's avatar
Markus Frank committed
        //--------------------------------
Markus Frank's avatar
Markus Frank committed
        const TGeoNode* node = pv.ptr();
Markus Frank's avatar
Markus Frank committed
        if ( !node ) {
Markus Frank's avatar
Markus Frank committed
          //	  std::cout <<  " *** findVolume: Invalid  placement:  - node pointer Null for volume:  " << pv.name() << std::endl ;
Markus Frank's avatar
Markus Frank committed
          throw std::runtime_error("*** findVolume: Invalid  placement:  - node pointer Null ! " + std::string( pv.name()  ) );
        }
        //	Volume vol = pv.volume();
Markus Frank's avatar
Markus Frank committed
        //	std::cout << "              ndau = " << node->GetNdaughters() << std::endl ;
Markus Frank's avatar
Markus Frank committed
        for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
Markus Frank's avatar
Markus Frank committed
          TGeoNode* daughter = node->GetDaughter(idau);
          PlacedVolume placement( daughter );
Markus Frank's avatar
Markus Frank committed
          if ( !placement.data() ) {
            throw std::runtime_error("*** findVolume: Invalid not instrumented placement:"+std::string(daughter->GetName())
                                     +" [Internal error -- bad detector constructor]");
          }
Markus Frank's avatar
Markus Frank committed
          PlacedVolume pv_dau = Ref_t(daughter); // why use a Ref_t here  ???
Markus Frank's avatar
Markus Frank committed
          if( findVolume(  pv_dau , theVol , volList ) ) {
Markus Frank's avatar
Markus Frank committed
            //	    std::cout << "  ----- found in daughter volume !!!  " << std::hex << pv_dau.volume().ptr() << std::endl ;
Markus Frank's avatar
Markus Frank committed
            return true ;
          } 
        }
Markus Frank's avatar
Markus Frank committed
        //  ------- not found:
Markus Frank's avatar
Markus Frank committed
        volList.pop_back() ;
Markus Frank's avatar
Markus Frank committed
        return false ;
        //--------------------------------
Frank Gaede's avatar
Frank Gaede committed
    //======================================================================================================================
    Surface::Surface( Geometry::DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ), 
Markus Frank's avatar
Markus Frank committed
                                                                       _wtM(0) , _id( 0) , _type( _volSurf.type() )  {
Frank Gaede's avatar
Frank Gaede committed
    long64 Surface::id() const { return _id ; }

    const SurfaceType& Surface::type() const { return _type ; }
Frank Gaede's avatar
Frank Gaede committed
    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() ; }
Frank Gaede's avatar
Frank Gaede committed

    /** Thickness of outer material */
    const IMaterial& Surface::innerMaterial() const {
      
      const IMaterial& mat = _volSurf.innerMaterial() ;
      if( ! ( mat.Z() > 0 ) ) {
Markus Frank's avatar
Markus Frank committed
        MaterialManager matMgr ;
        
	Vector3D p = _o - innerThickness() * _n  ;
Markus Frank's avatar
Markus Frank committed
        const MaterialVec& materials = matMgr.materialsBetween( _o , p  ) ;
	_volSurf.setInnerMaterial(  materials.size() > 1  ? 
				    matMgr.createAveragedMaterial( materials ) :
				    materials[0].first  )  ;
    const IMaterial& Surface::outerMaterial() const {
      
      const IMaterial& mat = _volSurf.outerMaterial() ;
      
      if( ! ( mat.Z() > 0 ) ) {
Markus Frank's avatar
Markus Frank committed
        MaterialManager matMgr ;
        
	Vector3D p = _o - outerThickness() * _n  ;
Markus Frank's avatar
Markus Frank committed
        const MaterialVec& materials = matMgr.materialsBetween( _o , p  ) ;
	_volSurf.setOuterMaterial(  materials.size() > 1  ? 
				    matMgr.createAveragedMaterial( materials ) :
				    materials[0].first  )  ;
Frank Gaede's avatar
Frank Gaede committed
      }
    Vector2D Surface::globalToLocal( const Vector3D& point) const {
Frank Gaede's avatar
Frank Gaede committed
      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 ;
    }


Frank Gaede's avatar
Frank Gaede committed
    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 ) ;
      //FG      return ( _volSurf.type().isPlane() ?   VolPlane(_volSurf).distance( localPoint )  : VolCylinder(_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) ;

      //FG      return ( _volSurf.type().isPlane() ?   VolPlane(_volSurf).insideBounds( localPoint, epsilon )  : VolCylinder(_volSurf).insideBounds( localPoint , epsilon) ) ;
Frank Gaede's avatar
Frank Gaede committed
    }

    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 ) ){
Markus Frank's avatar
Markus Frank committed
        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() ) ;
      // std::cout << " **** Surface::initialize() # placements for surface = " << pVList.size() 
      // 		<< " worldTransform : " 
      // 		<< std::endl ; 
      

      //=========== compute and cache world transform for surface ==========
      const TGeoHMatrix& wm = _det.worldTransformation() ;
#if 0 // debug
      for( std::list<PlacedVolume>::iterator it= pVList.begin(), n = pVList.end() ; it != n ; ++it ){
Markus Frank's avatar
Markus Frank committed
        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:
Markus Frank's avatar
Markus Frank committed
      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 ){

Markus Frank's avatar
Markus Frank committed
      	PlacedVolume pvol = *it ;
      	TGeoMatrix* m = pvol->GetMatrix();
      	// std::cout << "  +++ matrix for placed volume : " << std::endl ;
      	// m->Print() ;
Markus Frank's avatar
Markus Frank committed
        //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 ???
Markus Frank's avatar
Markus Frank committed
      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 ) ;

Frank Gaede's avatar
Frank Gaede committed
      // 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 ) ;
    
Frank Gaede's avatar
Frank Gaede committed
      //======== set the unique surface ID from the DetElement ( and placements below ? )
Frank Gaede's avatar
Frank Gaede committed

      // just use the DetElement ID for now ...
      _id = _det.volumeID() ;

      // 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->LocalToMasterVect( local_lines[i].first ,  av.array() ) ;
	  _wtM->LocalToMasterVect( 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() ;
Frank Gaede's avatar
Frank Gaede committed
      //      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() ) {
	
Markus Frank's avatar
Markus Frank committed
        if( shape->IsA() == TGeoBBox::Class() ) {
Markus Frank's avatar
Markus Frank committed
          TGeoBBox* box = ( TGeoBBox* ) shape  ;
Markus Frank's avatar
Markus Frank committed
          DDSurfaces::Vector3D boxDim(  box->GetDX() , box->GetDY() , box->GetDZ() ) ;  
Markus Frank's avatar
Markus Frank committed
          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
Markus Frank's avatar
Markus Frank committed
          if( isYZ || isXZ || isXY ) {  // plane is parallel to one of the box' sides -> need 4 vertices from box dimensions
Markus Frank's avatar
Markus Frank committed
            // if isYZ :
            unsigned uidx = 1 ;
            unsigned vidx = 2 ;
Markus Frank's avatar
Markus Frank committed
            DDSurfaces::Vector3D ubl(  0., 1., 0. ) ; 
            DDSurfaces::Vector3D vbl(  0., 0., 1. ) ;
Markus Frank's avatar
Markus Frank committed
            if( isXZ ) {
Markus Frank's avatar
Markus Frank committed
              ubl.fill( 1., 0., 0. ) ;
              vbl.fill( 0., 0., 1. ) ;
              uidx = 0 ;
              vidx = 2 ;
Markus Frank's avatar
Markus Frank committed
            } else if( isXY ) {
Markus Frank's avatar
Markus Frank committed
              ubl.fill( 1., 0., 0. ) ;
              vbl.fill( 0., 1., 0. ) ;
              uidx = 0 ;
              vidx = 1 ;
            }
Markus Frank's avatar
Markus Frank committed
            DDSurfaces::Vector3D ub ;
            DDSurfaces::Vector3D vb ;
            _wtM->LocalToMasterVect( ubl , ub.array() ) ;
            _wtM->LocalToMasterVect( vbl , vb.array() ) ;
Markus Frank's avatar
Markus Frank committed
            lines.reserve(4) ;
Markus Frank's avatar
Markus Frank committed
            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 ) ) ;
Markus Frank's avatar
Markus Frank committed
            return lines ;
          }	    
Markus Frank's avatar
Markus Frank committed
        } else if( shape->IsA() == TGeoConeSeg::Class() ) {
Markus Frank's avatar
Markus Frank committed
          TGeoCone* cone = ( TGeoCone* ) shape  ;
Markus Frank's avatar
Markus Frank committed
          // can only deal with special case of z-disk and origin in center of cone
          if( type().isZDisk() && lo.rho() < epsilon ) {
Markus Frank's avatar
Markus Frank committed
            double zhalf = cone->GetDZ() ;
            double rmax1 = cone->GetRmax1() ;
            double rmax2 = cone->GetRmax2() ;
            double rmin1 = cone->GetRmin1() ;
            double rmin2 = cone->GetRmin2() ;
Markus Frank's avatar
Markus Frank committed
            // 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()  ) ;  
Markus Frank's avatar
Markus Frank committed
            unsigned n = nMax / 4 ;
            double dPhi = 2.* ROOT::Math::Pi() / double( n ) ; 
Markus Frank's avatar
Markus Frank committed
            for( unsigned i = 0 ; i < n ; ++i ) {
Markus Frank's avatar
Markus Frank committed
              Vector3D rv00(  r0*sin(  i   *dPhi ) , r0*cos(  i   *dPhi )  , 0. ) ;
              Vector3D rv01(  r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi )  , 0. ) ;
Markus Frank's avatar
Markus Frank committed
              Vector3D rv10(  r1*sin(  i   *dPhi ) , r1*cos(  i   *dPhi )  , 0. ) ;
              Vector3D rv11(  r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi )  , 0. ) ;
Markus Frank's avatar
Markus Frank committed
              Vector3D pl0 =  lo + rv00 ;
              Vector3D pl1 =  lo + rv01 ;
Markus Frank's avatar
Markus Frank committed
              Vector3D pl2 =  lo + rv10 ;
              Vector3D pl3 =  lo + rv11 ;
Markus Frank's avatar
Markus Frank committed
              Vector3D pg0,pg1,pg2,pg3 ;
Markus Frank's avatar
Markus Frank committed
              _wtM->LocalToMaster( pl0, pg0.array() ) ;
              _wtM->LocalToMaster( pl1, pg1.array() ) ;
              _wtM->LocalToMaster( pl2, pg2.array() ) ;
              _wtM->LocalToMaster( pl3, pg3.array() ) ;
Markus Frank's avatar
Markus Frank committed
              lines.push_back( std::make_pair( pg0, pg1 ) ) ;
              lines.push_back( std::make_pair( pg2, pg3 ) ) ;
            }
Markus Frank's avatar
Markus Frank committed
            //add some vertical and horizontal lines so that the disc is seen in the rho-z projection
Markus Frank's avatar
Markus Frank committed
            n = 4 ; dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
Markus Frank's avatar
Markus Frank committed
            for( unsigned i = 0 ; i < n ; ++i ) {
Markus Frank's avatar
Markus Frank committed
              Vector3D rv0(  r0*sin( i * dPhi ) , r0*cos(  i * dPhi )  , 0. ) ;
              Vector3D rv1(  r1*sin( i * dPhi ) , r1*cos(  i * dPhi )  , 0. ) ;
Markus Frank's avatar
Markus Frank committed
              Vector3D pl0 =  lo + rv0 ;
              Vector3D pl1 =  lo + rv1 ;
Markus Frank's avatar
Markus Frank committed
              Vector3D pg0,pg1 ;
Markus Frank's avatar
Markus Frank committed
              _wtM->LocalToMaster( pl0, pg0.array() ) ;
              _wtM->LocalToMaster( pl1, pg1.array() ) ;
Markus Frank's avatar
Markus Frank committed
              lines.push_back( std::make_pair( pg0, pg1 ) ) ;
            }
Markus Frank's avatar
Markus Frank committed
          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() ) ;

          //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 AND y
        else if(shape->IsA() == TGeoTrd2::Class()){
          TGeoTrd2* trapezoid = ( TGeoTrd2* ) shape;
          //all lengths are half length
          double dx1 = trapezoid->GetDx1();
          double dx2 = trapezoid->GetDx2();
          //double dy1 = trapezoid->GetDy1();  
          //double dy2 = trapezoid->GetDy2();  
          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() ) ;

          //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;
        }
Markus Frank's avatar
Markus Frank committed
        // ===== default for arbitrary planes in arbitrary shapes ================= 
Markus Frank's avatar
Markus Frank committed
        // We create nMax vertices by rotating the local u vector around the normal
        // and checking the distance to the volume boundary in that direction.
        // This is brute force and not very smart, as many points are created on straight 
        // lines and the edges are still rounded. 
        // The alterative would be to compute the true intersections a plane and the most
        // common shapes - at least for boxes that should be not too hard. To be done...
Markus Frank's avatar
Markus Frank committed
        lines.reserve( nMax ) ;
Markus Frank's avatar
Markus Frank committed
        double dAlpha =  2.* ROOT::Math::Pi() / double( nMax ) ; 
Markus Frank's avatar
Markus Frank committed
        TVector3 norm( ln.x() , ln.y() , ln.z() ) ;
Markus Frank's avatar
Markus Frank committed
        DDSurfaces::Vector3D first, previous ;
Markus Frank's avatar
Markus Frank committed
        for(unsigned i=0 ; i< nMax ; ++i ){ 
Markus Frank's avatar
Markus Frank committed
          double alpha = double(i) * dAlpha ;
Markus Frank's avatar
Markus Frank committed
          TVector3 vec( lu.x() , lu.y() , lu.z() ) ;
Markus Frank's avatar
Markus Frank committed
          TRotation rot ;
          rot.Rotate( alpha , norm );
Markus Frank's avatar
Markus Frank committed
          TVector3 vecR = rot * vec ;
Markus Frank's avatar
Markus Frank committed
          DDSurfaces::Vector3D luRot ;
          luRot.fill( vecR ) ;
Markus Frank's avatar
Markus Frank committed
          double dist = shape->DistFromInside( const_cast<double*> (lo.const_array()) , const_cast<double*> (luRot.const_array())  , 3, 0.1 ) ;
Markus Frank's avatar
Markus Frank committed
          // local point at volume boundary
          DDSurfaces::Vector3D lp = lo + dist * luRot ;
Markus Frank's avatar
Markus Frank committed
          DDSurfaces::Vector3D gp ;
Markus Frank's avatar
Markus Frank committed
          _wtM->LocalToMaster( lp , gp.array() ) ;

          // std::cout << " **** normal:" << ln << " lu:" << lu  << " alpha:" << alpha << " luRot:" << luRot << " lp :" << lp  << " gp:" << gp << " dist : " << dist 
          // 	    << " is point " << gp << " inside : " << shape->Contains( gp.const_array()  )  
          // 	    << " dist from outside for lo,lu " <<  shape->DistFromOutside( lo.const_array()  , lu.const_array()   , 3 )    
          // 	    << " dist from inside for lo,ln " <<  shape->DistFromInside( lo.const_array()  , ln.const_array()   , 3 )    
          // 	    << std::endl;
          //	  shape->Dump() ;
Markus Frank's avatar
Markus Frank committed
          if( i >  0 ) 
            lines.push_back( std::make_pair( previous, gp )  ) ;
          else
            first = gp ;
Markus Frank's avatar
Markus Frank committed
          previous = gp ;
        }
        lines.push_back( std::make_pair( previous, first )  ) ;
      } else if( type().isCylinder() ) {  
Markus Frank's avatar
Markus Frank committed
        //	if( shape->IsA() == TGeoTube::Class() ) {
        if( shape->IsA() == TGeoConeSeg::Class() ) {
Markus Frank's avatar
Markus Frank committed
          lines.reserve( nMax ) ;
Markus Frank's avatar
Markus Frank committed
          TGeoTube* tube = ( TGeoTube* ) shape  ;
Markus Frank's avatar
Markus Frank committed
          double zHalf = tube->GetDZ() ;
Markus Frank's avatar
Markus Frank committed
          Vector3D zv( 0. , 0. , zHalf ) ;
Markus Frank's avatar
Markus Frank committed
          double r = lo.rho() ;
Markus Frank's avatar
Markus Frank committed
          unsigned n = nMax / 4 ;
          double dPhi = 2.* ROOT::Math::Pi() / double( n ) ; 
Markus Frank's avatar
Markus Frank committed
          for( unsigned i = 0 ; i < n ; ++i ) {
Markus Frank's avatar
Markus Frank committed
            Vector3D rv0(  r*sin(  i   *dPhi ) , r*cos(  i   *dPhi )  , 0. ) ;
            Vector3D rv1(  r*sin( (i+1)*dPhi ) , r*cos( (i+1)*dPhi )  , 0. ) ;
Markus Frank's avatar
Markus Frank committed
            // 4 points on local cylinder