Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-10-05 08:17:21 

0001 //==========================================================================
0002 //  AIDA Detector description implementation 
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : F.Gaede
0011 //
0012 //==========================================================================
0013 #include "DDRec/Surface.h"
0014 #include "DD4hep/detail/DetectorInterna.h"
0015 #include "DD4hep/Memory.h"
0016 
0017 #include "DDRec/MaterialManager.h"
0018 
0019 #include <cmath>
0020 #include <memory>
0021 
0022 #include "TGeoMatrix.h"
0023 #include "TGeoShape.h"
0024 #include "TRotation.h"
0025 //TGeoTrd1 is apparently not included by default
0026 #include "TGeoTrd1.h"
0027 
0028 namespace dd4hep {
0029   namespace rec {
0030  
0031     using namespace detail ;
0032 
0033 
0034     //======================================================================================================
0035   
0036     void VolSurfaceBase::setU(const Vector3D& u_val) {  _u = u_val  ; }
0037     void VolSurfaceBase::setV(const Vector3D& v_val) {  _v = v_val ; }
0038     void VolSurfaceBase::setNormal(const Vector3D& n) { _n = n ; }
0039     void VolSurfaceBase::setOrigin(const Vector3D& o) { _o = o ; }
0040     
0041     long64 VolSurfaceBase::id() const  { return _id ; } 
0042 
0043     const SurfaceType& VolSurfaceBase::type() const { return _type ; }
0044     Vector3D VolSurfaceBase::u(const Vector3D& /*point*/) const { return _u ; }
0045     Vector3D VolSurfaceBase::v(const Vector3D& /*point*/) const { return _v ; }
0046     Vector3D VolSurfaceBase::normal(const Vector3D& /*point*/) const { return _n ; }
0047     const Vector3D& VolSurfaceBase::origin() const { return _o ;}
0048 
0049     Vector2D VolSurfaceBase::globalToLocal( const Vector3D& point) const {
0050 
0051       Vector3D p = point - origin() ;
0052 
0053       // create new orthogonal unit vectors
0054       // FIXME: these vectors should be cached really ... 
0055 
0056       double uv = u() * v() ;
0057       Vector3D uprime = ( u() - uv * v() ).unit() ; 
0058       Vector3D vprime = ( v() - uv * u() ).unit() ; 
0059       double uup = u() * uprime ;
0060       double vvp = v() * vprime ;
0061       
0062       return  Vector2D(   p*uprime / uup ,  p*vprime / vvp ) ;
0063     }
0064     
0065     Vector3D VolSurfaceBase::localToGlobal( const Vector2D& point) const {
0066 
0067       Vector3D g = origin() + point[0] * u() + point[1] * v() ;
0068 
0069       return g ;
0070     }
0071 
0072     const IMaterial&  VolSurfaceBase::innerMaterial() const { return  _innerMat ; }
0073     const IMaterial&  VolSurfaceBase::outerMaterial() const { return  _outerMat ; }
0074     double VolSurfaceBase::innerThickness() const { return _th_i ; }
0075     double VolSurfaceBase::outerThickness() const { return _th_o ; }
0076     
0077 
0078     double VolSurfaceBase::length_along_u() const {
0079       
0080       const Vector3D& o = this->origin() ;
0081       const Vector3D& u_val = this->u( o ) ;      
0082       Vector3D  um = -1. * u_val ;
0083       
0084       double dist_p = 0. ;
0085       double dist_m = 0. ;
0086       
0087 
0088       // std::cout << " VolSurfaceBase::length_along_u() : o =  " << o << " u = " <<    this->u( o ) 
0089       //        << " -u = " << um << std::endl ;
0090 
0091 
0092       if( volume()->GetShape()->Contains( o.const_array() ) ){
0093 
0094         dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
0095                                                        const_cast<double*> ( u_val.const_array() ) ) ;
0096         dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
0097                                                        const_cast<double*> ( um.array()      ) ) ;
0098     
0099 
0100         // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
0101         //    << " dist_p " <<    dist_p
0102         //    << " dist_m " <<    dist_m
0103         //    << std::endl ;
0104     
0105 
0106       } else{
0107 
0108         dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
0109                                                         const_cast<double*> ( u_val.const_array() ) ) ;
0110         dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
0111                                                         const_cast<double*> ( um.array()      ) ) ;
0112 
0113         dist_p *= 1.0001 ;
0114         dist_m *= 1.0001 ;
0115 
0116         // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
0117         //    << " dist_p " <<    dist_p
0118         //    << " dist_m " <<    dist_m
0119         //    << std::endl ;
0120 
0121         Vector3D o_1 = this->origin() + dist_p * u_val ;
0122         Vector3D o_2 = this->origin() + dist_m * um ;
0123 
0124         dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) , 
0125                                                         const_cast<double*> ( u_val.const_array() ) ) ;
0126 
0127         dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) , 
0128                                                         const_cast<double*> ( um.array()      ) ) ;
0129 
0130         // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
0131         //    << " dist_p " <<    dist_p
0132         //    << " dist_m " <<    dist_m
0133         //    << std::endl ;
0134       }
0135     
0136       return dist_p + dist_m ;
0137 
0138 
0139     }
0140     
0141     double VolSurfaceBase::length_along_v() const {
0142 
0143       const Vector3D& o = this->origin() ;
0144       const Vector3D& v_val = this->v( o ) ;      
0145       Vector3D  vm = -1. * v_val ;
0146       
0147       double dist_p = 0. ;
0148       double dist_m = 0. ;
0149       
0150 
0151       // std::cout << " VolSurfaceBase::length_along_u() : o =  " << o << " u = " <<    this->u( o ) 
0152       //        << " -u = " << vm << std::endl ;
0153 
0154 
0155       if( volume()->GetShape()->Contains( o.const_array() ) ){
0156 
0157         dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
0158                                                        const_cast<double*> ( v_val.const_array() ) ) ;
0159         dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) , 
0160                                                        const_cast<double*> ( vm.array()      ) ) ;
0161     
0162 
0163         // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
0164         //    << " dist_p " <<    dist_p
0165         //    << " dist_m " <<    dist_m
0166         //    << std::endl ;
0167     
0168 
0169       } else{
0170 
0171         dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
0172                                                         const_cast<double*> ( v_val.const_array() ) ) ;
0173         dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) , 
0174                                                         const_cast<double*> ( vm.array()      ) ) ;
0175 
0176         dist_p *= 1.0001 ;
0177         dist_m *= 1.0001 ;
0178 
0179         // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
0180         //    << " dist_p " <<    dist_p
0181         //    << " dist_m " <<    dist_m
0182         //    << std::endl ;
0183 
0184         Vector3D o_1 = this->origin() + dist_p * v_val ;
0185         Vector3D o_2 = this->origin() + dist_m * vm ;
0186 
0187         dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) , 
0188                                                         const_cast<double*> ( v_val.const_array() ) ) ;
0189 
0190         dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) , 
0191                                                         const_cast<double*> ( vm.array()      ) ) ;
0192 
0193         // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o)  =  " << volume()->GetShape()->Contains( o.const_array() )
0194         //    << " dist_p " <<    dist_p
0195         //    << " dist_m " <<    dist_m
0196         //    << std::endl ;
0197       }
0198     
0199       return dist_p + dist_m ;
0200 
0201     }
0202     
0203 
0204     double VolSurfaceBase::distance(const Vector3D& /*point*/ ) const { return 1.e99 ; }
0205 
0206     /// Checks if the given point lies within the surface
0207     bool VolSurfaceBase::insideBounds(const Vector3D& point, double epsilon) const {
0208 
0209 #if 0
0210 
0211       bool inShape = ( type().isUnbounded() ? true : volume()->GetShape()->Contains( point.const_array() ) ) ;
0212 
0213       double dist = std::abs ( distance( point ) ) ;
0214 
0215       std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << dist 
0216                 << " origin = " << origin() << " normal = " << normal() 
0217                 << " p * n = " << point * normal() 
0218                 << " isInShape : " << inShape << std::endl ;
0219 
0220       return dist < epsilon && inShape ;
0221 #else
0222 
0223       if(  type().isUnbounded() ){
0224 
0225         return std::abs ( distance( point ) ) < epsilon ;
0226 
0227       } else {
0228 
0229         return (  std::abs ( distance( point ) ) < epsilon &&  volume()->GetShape()->Contains( point.const_array() ) ) ;
0230       }
0231 
0232 #endif
0233  
0234     }
0235 
0236 
0237     std::vector< std::pair<Vector3D, Vector3D> > VolSurfaceBase::getLines(unsigned ) {
0238       // dummy implementation returning empty set
0239       std::vector< std::pair<Vector3D, Vector3D> >  lines ;
0240       return lines ;
0241     }
0242 
0243     //===================================================================
0244     // simple wrapper methods forwarding the call to the implementation object
0245 
0246     long64 VolSurface::id() const  { return _surf->id() ; } 
0247     const SurfaceType& VolSurface::type() const { return _surf->type() ; }
0248     Vector3D VolSurface::u( const Vector3D& point ) const {  return _surf->u(point) ; }
0249     Vector3D VolSurface::v(const Vector3D& point  ) const {  return _surf->v(point) ; }
0250     Vector3D VolSurface::normal(const Vector3D& point ) const {  return _surf->normal(point) ; }
0251     const Vector3D& VolSurface::origin() const { return _surf->origin() ;}
0252     Vector2D VolSurface::globalToLocal( const Vector3D& point) const { return _surf->globalToLocal( point ) ; }
0253     Vector3D VolSurface::localToGlobal( const Vector2D& point) const { return _surf->localToGlobal( point) ; }
0254     const IMaterial&  VolSurface::innerMaterial() const{ return _surf->innerMaterial() ; }
0255     const IMaterial&  VolSurface::outerMaterial() const  { return _surf->outerMaterial()  ; }
0256     double VolSurface::innerThickness() const { return _surf->innerThickness() ; }
0257     double VolSurface::outerThickness() const { return _surf->outerThickness() ; }
0258     double VolSurface::length_along_u() const { return _surf->length_along_u() ; }
0259     double VolSurface::length_along_v() const { return _surf->length_along_v() ; }
0260     double VolSurface::distance(const Vector3D& point ) const  {  return _surf->distance( point ) ; }
0261     bool VolSurface::insideBounds(const Vector3D& point, double epsilon) const {
0262       return _surf->insideBounds( point, epsilon ) ; 
0263     }
0264     std::vector< std::pair<Vector3D, Vector3D> > VolSurface::getLines(unsigned nMax) {
0265       return _surf->getLines(nMax) ;
0266     }
0267 
0268     //===================================================================
0269 
0270     /** Distance to planar surface */
0271     double VolPlaneImpl::distance(const Vector3D& point ) const {
0272       return ( point - origin() ) *  normal()  ;
0273     }
0274     //======================================================================================================
0275 
0276     VolCylinderImpl::VolCylinderImpl( Volume vol, SurfaceType typ, 
0277                                       double thickness_inner ,double thickness_outer,  Vector3D o ) :
0278 
0279       VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , Vector3D() , Vector3D() , o , vol, 0) {
0280       Vector3D v_val( 0., 0., 1. ) ;
0281       Vector3D o_rphi( o.x() , o.y() , 0. ) ;
0282       Vector3D n =  o_rphi.unit() ; 
0283       Vector3D u_val = v_val.cross( n ) ;
0284 
0285       setU( u_val ) ;
0286       setV( v_val ) ;
0287       setNormal( n ) ;
0288 
0289       _type.setProperty( SurfaceType::Plane    , false ) ;
0290       _type.setProperty( SurfaceType::Cylinder , true ) ;
0291       _type.setProperty( SurfaceType::Cone     , false ) ;
0292       _type.checkParallelToZ( *this ) ;
0293       _type.checkOrthogonalToZ( *this ) ;
0294     }      
0295 
0296     Vector3D VolCylinderImpl::u(const Vector3D& point ) const {  
0297 
0298       Vector3D n( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
0299 
0300       return v().cross( n ) ; 
0301     }
0302     
0303     Vector3D VolCylinderImpl::normal(const Vector3D& point ) const {  
0304 
0305       // normal is just given by phi of the point 
0306       return Vector3D( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
0307     }
0308 
0309     Vector2D VolCylinderImpl::globalToLocal( const Vector3D& point) const {
0310       
0311       // cylinder is parallel to v here so u is Z and v is r *Phi
0312       double phi = point.phi() - origin().phi() ;
0313       
0314       while( phi < -M_PI ) phi += 2.*M_PI ;
0315       while( phi >  M_PI ) phi -= 2.*M_PI ;
0316       
0317       return  Vector2D( origin().rho() * phi, point.z() - origin().z() ) ;
0318     }
0319     
0320     
0321     Vector3D VolCylinderImpl::localToGlobal( const Vector2D& point) const {
0322       
0323       double z = point.v() + origin().z() ;
0324       double phi = point.u() / origin().rho() + origin().phi() ;
0325       
0326       while( phi < -M_PI ) phi += 2.*M_PI ;
0327       while( phi >  M_PI ) phi -= 2.*M_PI ;
0328       
0329       return Vector3D( origin().rho() , phi, z  , Vector3D::cylindrical )    ;
0330     }
0331 
0332 
0333     /** Distance to surface */
0334     double VolCylinderImpl::distance(const Vector3D& point ) const {
0335       
0336       return point.rho() - origin().rho()  ;
0337     }
0338     
0339     //================================================================================================================
0340     VolConeImpl::VolConeImpl( Volume vol, SurfaceType typ, 
0341                               double thickness_inner ,double thickness_outer, Vector3D v_val,  Vector3D o_val ) :
0342       
0343       VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , v_val ,  Vector3D() , Vector3D() , vol, 0) {
0344 
0345       Vector3D o_rphi( o_val.x() , o_val.y() , 0. ) ;
0346 
0347       // sanity check: v and o have to have a common phi
0348       double dphi = v_val.phi() - o_rphi.phi() ;
0349       while( dphi < -M_PI ) dphi += 2.*M_PI ;
0350       while( dphi >  M_PI ) dphi -= 2.*M_PI ;
0351 
0352       if( std::fabs( dphi ) > 1e-6 ){
0353         std::stringstream sst ; sst << "VolConeImpl::VolConeImpl() - incompatibel vector v and o given " 
0354                                     << v_val << " - " << o_val ;
0355         throw std::runtime_error( sst.str() ) ;
0356       }
0357       
0358       double theta = v_val.theta() ;
0359 
0360       Vector3D n( 1. , v_val.phi() , ( theta + M_PI/2. ) , Vector3D::spherical ) ;
0361       Vector3D u_val = v_val.cross( n ) ;
0362 
0363       setU( u_val ) ;
0364       setOrigin( o_rphi ) ;
0365       setNormal( n ) ;
0366 
0367       // set helper variable for faster computations (describe cone with tip at origin)
0368       _tanTheta = std::tan( theta ) ; 
0369       double tipoffset = o_val.rho() / _tanTheta ; // distance from tip to origin.z()
0370       _ztip     = o_val.z()  - tipoffset ;
0371 
0372       double dist_p = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) , 
0373                                                        const_cast<double*> ( v_val.const_array() ) ) ;
0374       Vector3D vm = -1. * v_val ;
0375       double dist_m = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) , 
0376                                                        const_cast<double*> ( vm.array()      ) ) ;
0377 
0378       double costh = std::cos( theta) ;
0379       _zt0 = tipoffset - dist_m *  costh ;           
0380       _zt1 = tipoffset + dist_p *  costh ;     
0381 
0382 
0383       _type.setProperty( SurfaceType::Plane   , false ) ;
0384       _type.setProperty( SurfaceType::Cylinder, false ) ;
0385       _type.setProperty( SurfaceType::Cone    , true ) ;
0386       _type.setProperty( SurfaceType::ParallelToZ, true ) ;
0387       _type.setProperty( SurfaceType::OrthogonalToZ, false ) ;
0388     }      
0389 
0390     
0391     Vector3D VolConeImpl::v(const Vector3D& point ) const {  
0392       // just take phi from point
0393       Vector3D av( 1. , point.phi() , _v.theta() , Vector3D::spherical ) ;
0394       return av ; 
0395     }
0396     
0397     Vector3D VolConeImpl::u(const Vector3D& point ) const {  
0398       // compute from v X n 
0399       const Vector3D& av = this->v( point ) ;
0400       const Vector3D& n = normal( point ) ;
0401       return av.cross( n ) ; 
0402     }
0403     
0404     Vector3D VolConeImpl::normal(const Vector3D& point ) const {  
0405       // just take phi from point
0406       Vector3D n( 1. , point.phi() , _n.theta() , Vector3D::spherical ) ;
0407       return n ;
0408     }
0409 
0410     Vector2D VolConeImpl::globalToLocal( const Vector3D& point) const {
0411       
0412       // cone is parallel to z here, so u is r *Phi and v is "along" z
0413       double phi = point.phi() - origin().phi() ;
0414       
0415       while( phi < -M_PI ) phi += 2.*M_PI ;
0416       while( phi >  M_PI ) phi -= 2.*M_PI ;
0417       
0418 
0419       double r = ( point.z() - _ztip ) * _tanTheta ;
0420 
0421       return  Vector2D(  r*phi, ( point.z() - origin().z() ) / cos( _v.theta() ) ) ;
0422     }
0423     
0424     
0425     Vector3D VolConeImpl::localToGlobal( const Vector2D& point) const {
0426       
0427       double z = point.v() * cos( _v.theta() ) + origin().z() ;
0428       
0429       double r =  ( z - _ztip ) * _tanTheta ;
0430 
0431       double phi = point.u() / r + origin().phi() ;
0432       
0433       while( phi < -M_PI ) phi += 2.*M_PI ;
0434       while( phi >  M_PI ) phi -= 2.*M_PI ;
0435       
0436       return Vector3D( r , phi, z  , Vector3D::cylindrical )    ;
0437     }
0438 
0439 
0440     /** Distance to surface */
0441     double VolConeImpl::distance(const Vector3D& point ) const {
0442 
0443       // // if the point is in the other hemispere we return the distance to origin 
0444       // // -> this assumes that the cones do not cross the xy-plane ...
0445       // // otherwise we get the distance to the mirrored part of the cone
0446       // // needs more thought ..
0447       // if( origin().z() * point.z() < 0. ) 
0448       //    return point.r() ;
0449 
0450       //fixme: there are probably faster ways to compute this
0451       // e.g by using the intercept theorem - tbd. ...
0452       // const Vector2D& lp = globalToLocal( point ) ;
0453       // const Vector3D& gp = localToGlobal( lp ) ;
0454 
0455       // Vector3D dz = point - gp ;
0456 
0457       //return dz * normal( point )   ;
0458 
0459       double zp = point.z() - _ztip ;
0460       double r = point.rho() - zp * _tanTheta ;
0461       return r * std::cos( _v.theta() ) ;
0462 
0463     }
0464     
0465     /// create outer bounding lines for the given symmetry of the polyhedron
0466     std::vector< std::pair<Vector3D, Vector3D> >  VolConeImpl::getLines(unsigned nMax){
0467       
0468       std::vector< std::pair<Vector3D, Vector3D> >  lines ;
0469       
0470       lines.reserve( nMax ) ;
0471       
0472       double theta = v().theta() ;
0473       double half_length = 0.5 * length_along_v() * cos( theta ) ;
0474 
0475       Vector3D zv( 0. , 0. , half_length ) ;
0476 
0477       double dr =  half_length * tan( theta ) ;
0478 
0479       double r0 = origin().rho() - dr ;  
0480       double r1 = origin().rho() + dr ;
0481       
0482 
0483       unsigned n = nMax / 4 ;
0484       double dPhi = 2.* ROOT::Math::Pi() / double( n ) ; 
0485       
0486       for( unsigned i = 0 ; i < n ; ++i ) {
0487     
0488         Vector3D r0v0(  r0*sin(  i   *dPhi ) , r0*cos(  i   *dPhi )  , 0. ) ;
0489         Vector3D r0v1(  r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi )  , 0. ) ;
0490 
0491         Vector3D r1v0(  r1*sin(  i   *dPhi ) , r1*cos(  i   *dPhi )  , 0. ) ;
0492         Vector3D r1v1(  r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi )  , 0. ) ;
0493     
0494         Vector3D pl0 =  zv + r1v0 ;
0495         Vector3D pl1 =  zv + r1v1 ;
0496         Vector3D pl2 = -zv + r0v1  ;
0497         Vector3D pl3 = -zv + r0v0 ;
0498     
0499         lines.emplace_back( pl0, pl1 );
0500         lines.emplace_back( pl1, pl2 );
0501         lines.emplace_back( pl2, pl3 );
0502         lines.emplace_back( pl3, pl0 );
0503       } 
0504       return lines; 
0505     }
0506     
0507     //================================================================================================================
0508     
0509     SurfaceList::~SurfaceList(){
0510       if( _isOwner ) {
0511         // delete all surfaces attached to this volume
0512         std::for_each(begin(), end(), detail::deleteObject<ISurface>);
0513       }
0514     }
0515 
0516     //================================================================================================================
0517 
0518     VolSurfaceList* volSurfaceList( const DetElement& det ) {
0519       VolSurfaceList* list = det.extension< VolSurfaceList >(false);
0520       if ( !list )  {
0521         list = det.addExtension<VolSurfaceList >(new VolSurfaceList);
0522       }
0523       return list ;
0524     }
0525 
0526 
0527     //======================================================================================================================
0528 
0529     bool findVolume( PlacedVolume pv,  Volume theVol, std::list< PlacedVolume >& volList ) {
0530       
0531 
0532       volList.emplace_back( pv ) ;
0533       
0534       //   unsigned count = volList.size() ;
0535       //   for(unsigned i=0 ; i < count ; ++i) {
0536       //    std::cout << " **" ;
0537       //   }
0538       //   std::cout << " searching for volume: " << theVol.name() << " " << std::hex << theVol.ptr() << "  <-> pv.volume : "  << pv.name() << " " <<  pv.volume().ptr() 
0539       //            << " pv.volume().ptr() == theVol.ptr() " <<  (pv.volume().ptr() == theVol.ptr() )
0540       //            << std::endl ;
0541       
0542 
0543       if( pv.volume().ptr() == theVol.ptr() ) { 
0544     
0545         return true ;
0546     
0547       } else {
0548     
0549         //--------------------------------
0550 
0551         const TGeoNode* node = pv.ptr();
0552     
0553         if ( !node ) {
0554       
0555           //      std::cout <<  " *** findVolume: Invalid  placement:  - node pointer Null for volume:  " << pv.name() << std::endl ;
0556 
0557           throw std::runtime_error("*** findVolume: Invalid  placement:  - node pointer Null ! " + std::string( pv.name()  ) );
0558         }
0559         //  Volume vol = pv.volume();
0560     
0561         //  std::cout << "              ndau = " << node->GetNdaughters() << std::endl ;
0562 
0563         for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
0564       
0565           TGeoNode* daughter = node->GetDaughter(idau);
0566           PlacedVolume placement( daughter );
0567       
0568           if ( !placement.data() ) {
0569             throw std::runtime_error("*** findVolume: Invalid not instrumented placement:"+std::string(daughter->GetName())
0570                                      +" [Internal error -- bad detector constructor]");
0571           }
0572       
0573           PlacedVolume pv_dau( daughter );
0574 
0575           if( findVolume(  pv_dau , theVol , volList ) ) {
0576         
0577             //      std::cout << "  ----- found in daughter volume !!!  " << std::hex << pv_dau.volume().ptr() << std::endl ;
0578 
0579             return true ;
0580           } 
0581         }
0582 
0583         //  ------- not found:
0584 
0585         volList.pop_back() ;
0586 
0587         return false ;
0588         //--------------------------------
0589 
0590       }
0591     } 
0592 
0593     //======================================================================================================================
0594 
0595     Surface::Surface( DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ), 
0596                                                              _wtM() , _id( 0) , _type( _volSurf.type() )  {
0597 
0598       initialize() ;
0599     }      
0600  
0601     long64 Surface::id() const { return _id ; }
0602 
0603     const SurfaceType& Surface::type() const { return _type ; }
0604 
0605     Vector3D Surface::u(const Vector3D& /*point*/) const { return _u ; }
0606     Vector3D Surface::v(const Vector3D& /*point*/) const { return _v ; }
0607     Vector3D Surface::normal(const Vector3D& /*point*/) const { return _n ; }
0608     const Vector3D& Surface::origin() const { return _o ;}
0609     double Surface::innerThickness() const { return _volSurf.innerThickness() ; }
0610     double Surface::outerThickness() const { return _volSurf.outerThickness() ; }
0611     double Surface::length_along_u() const { return _volSurf.length_along_u() ; }
0612     double Surface::length_along_v() const { return _volSurf.length_along_v() ; }
0613 
0614     /** Thickness of outer material */
0615     
0616     const IMaterial& Surface::innerMaterial() const {
0617       
0618       const IMaterial& mat = _volSurf.innerMaterial() ;
0619       
0620       if( mat.Z() <= 0 ) {
0621     
0622         MaterialManager matMgr( _det.placement().volume() )  ;
0623         
0624         Vector3D p = _o - innerThickness() * _n  ;
0625 
0626         const MaterialVec& materials = matMgr.materialsBetween( _o , p  ) ;
0627 
0628         _volSurf.setInnerMaterial(  materials.size() > 1  ? 
0629                                     matMgr.createAveragedMaterial( materials ) :
0630                                     materials[0].first  )  ;
0631       }
0632       return  mat ;
0633     }
0634 
0635     const IMaterial& Surface::outerMaterial() const {
0636       
0637       const IMaterial& mat = _volSurf.outerMaterial() ;
0638       
0639       if( mat.Z() <= 0 ) {
0640     
0641         MaterialManager matMgr( _det.placement().volume() ) ;
0642         
0643         Vector3D p = _o + outerThickness() * _n  ;
0644 
0645         const MaterialVec& materials = matMgr.materialsBetween( _o , p  ) ;
0646 
0647         _volSurf.setOuterMaterial(  materials.size() > 1  ? 
0648                                     matMgr.createAveragedMaterial( materials ) :
0649                                     materials[0].first  )  ;
0650       }
0651       return  mat ;
0652     }
0653       
0654 
0655     Vector2D Surface::globalToLocal( const Vector3D& point) const {
0656 
0657       Vector3D p = point - origin() ;
0658 
0659       // create new orthogonal unit vectors
0660       // FIXME: these vectors should be cached really ... 
0661 
0662       double uv = u() * v() ;
0663       Vector3D uprime = ( u() - uv * v() ).unit() ; 
0664       Vector3D vprime = ( v() - uv * u() ).unit() ; 
0665       double uup = u() * uprime ;
0666       double vvp = v() * vprime ;
0667       
0668       return  Vector2D(   p*uprime / uup ,  p*vprime / vvp ) ;
0669     }
0670     
0671     
0672     Vector3D Surface::localToGlobal( const Vector2D& point) const {
0673 
0674       Vector3D g = origin() + point[0] * u() + point[1] * v() ;
0675       return g ;
0676     }
0677 
0678 
0679     Vector3D Surface::volumeOrigin() const { 
0680 
0681       double o_array[3] ;
0682 
0683       _wtM->LocalToMaster    ( Vector3D() , o_array ) ;
0684       
0685       Vector3D o(o_array) ;
0686      
0687       return o ;
0688     }
0689 
0690 
0691     double Surface::distance(const Vector3D& point ) const {
0692 
0693       double pa[3] ;
0694       _wtM->MasterToLocal( point , pa ) ;
0695       Vector3D localPoint( pa ) ;
0696       
0697       return _volSurf.distance( localPoint ) ;
0698     }
0699       
0700     bool Surface::insideBounds(const Vector3D& point, double epsilon) const {
0701 
0702       double pa[3] ;
0703       _wtM->MasterToLocal( point , pa ) ;
0704       Vector3D localPoint( pa ) ;
0705       
0706       return _volSurf.insideBounds( localPoint , epsilon) ;
0707     }
0708 
0709     void Surface::initialize() {
0710       
0711       // first we need to find the right volume for the local surface in the DetElement's volumes
0712       std::list< PlacedVolume > pVList ;
0713       PlacedVolume pv = _det.placement() ;
0714       Volume   theVol = _volSurf.volume() ;
0715       
0716       if( ! findVolume(  pv, theVol , pVList ) ){
0717         theVol = _volSurf.volume() ;
0718         std::stringstream sst ; sst << " ***** ERROR: Volume " << theVol.name() << " not found for DetElement " << _det.name()  << " with surface "  ;
0719         throw std::runtime_error( sst.str() ) ;
0720       } 
0721 
0722       //=========== compute and cache world transform for surface ==========
0723       Alignment nominal = _det.nominal();
0724       const TGeoHMatrix& wm = nominal.worldTransformation() ;
0725       
0726 #if 0 // debug
0727       wm.Print() ;
0728       for( std::list<PlacedVolume>::iterator it= pVList.begin(), n = pVList.end() ; it != n ; ++it ){
0729         PlacedVolume pv = *it ;
0730         TGeoMatrix* m = pv->GetMatrix();
0731         std::cout << "  +++ matrix for placed volume : " << std::endl ;
0732         m->Print() ;
0733       }
0734 #endif
0735 
0736       // need to get the inverse transformation ( see Detector.cpp )
0737       // std::auto_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm.Inverse() ) ) ;
0738       // has been fixed now, no need to get the inverse anymore:
0739       std::unique_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm ) ) ;
0740 
0741       //---- if the volSurface is not in the DetElement's volume, we need to mutliply the path to the volume to the
0742       // DetElements world transform
0743       for( auto it = std::next(pVList.begin()) ; it != pVList.end() ; ++it ) {
0744 
0745         PlacedVolume pvol = *it ;
0746         TGeoMatrix* m = pvol->GetMatrix();
0747         // std::cout << "  +++ matrix for placed volume : " << std::endl ;
0748         // m->Print() ;
0749         //wtI->MultiplyLeft( m );
0750 
0751         wtI->Multiply( m );
0752       }
0753 
0754       //      std::cout << "  +++ new world transform matrix  : " << std::endl ;
0755 
0756 #if 0 //fixme: which convention to use here - the correct should be wtI, however it is the inverse of what is stored in DetElement ???
0757       dd4hep_ptr<TGeoHMatrix> wt( new TGeoHMatrix( wtI->Inverse() ) ) ;
0758       wt->Print() ;
0759       // cache the world transform for the surface
0760       _wtM = std::move(wtI);
0761 #else
0762       //      wtI->Print() ;
0763       // cache the world transform for the surface
0764       _wtM = std::move(wtI);
0765 #endif
0766 
0767 
0768       //  ============ now fill the global surface vectors ==========================
0769 
0770       double ua[3], va[3], na[3], oa[3] ;
0771 
0772       _wtM->LocalToMasterVect( _volSurf.u()      , ua ) ;
0773       _wtM->LocalToMasterVect( _volSurf.v()      , va ) ;
0774       _wtM->LocalToMasterVect( _volSurf.normal() , na ) ;
0775       _wtM->LocalToMaster    ( _volSurf.origin() , oa ) ;
0776 
0777       _u.fill( ua ) ;
0778       _v.fill( va ) ;
0779       _n.fill( na ) ;
0780       _o.fill( oa ) ;
0781 
0782       // std::cout << " --- local and global surface vectors : ------- " << std::endl 
0783       //            << "    u : " << _volSurf.u()       << "  -  " << _u << std::endl 
0784       //            << "    v : " << _volSurf.v()       << "  -  " << _v << std::endl 
0785       //            << "    n : " << _volSurf.normal()  << "  -  " << _n << std::endl 
0786       //            << "    o : " << _volSurf.origin()  << "  -  " << _o << std::endl ;
0787       
0788 
0789       //  =========== check parallel and orthogonal to Z ===================
0790       
0791       if( ! _type.isCone() ) { 
0792         //fixme: workaround for conical surfaces that should always be parallel to z
0793         //       however the check with the normal does not work here ...
0794 
0795         _type.checkParallelToZ( *this ) ;
0796     
0797         _type.checkOrthogonalToZ( *this ) ;
0798       }
0799       
0800       //======== set the unique surface ID from the DetElement ( and placements below ? )
0801 
0802       // just use the DetElement ID for now ...
0803       // or the id set by the user to the VolSurface ...
0804       _id = ( _volSurf.id()==0 ?  _det.volumeID() : _volSurf.id() ) ;
0805 
0806       // typedef PlacedVolume::VolIDs IDV ;
0807       // DetElement d = _det ;
0808       // while( d.isValid() &&  d.parent().isValid() ){
0809       //    PlacedVolume pv = d.placement() ;
0810       //    if( pv.isValid() ){
0811       //      const IDV& idV = pv.volIDs() ; 
0812       //      std::cout << " VolIDs : " << d.name() << std::endl ;
0813       //      for( unsigned i=0, n=idV.size() ; i<n ; ++i){
0814       //        std::cout  << "  " << idV[i].first << " - " << idV[i].second << std::endl ;
0815       //      }
0816       //    }
0817       //    d = d.parent() ;
0818       // }
0819      
0820     }
0821     //===================================================================================================================
0822       
0823     std::vector< std::pair<Vector3D, Vector3D> > Surface::getLines(unsigned nMax) {
0824 
0825 
0826       const static double epsilon = 1e-6 ; 
0827 
0828       std::vector< std::pair<Vector3D, Vector3D> > lines ;
0829 
0830       //--------------------------------------------
0831       // check if there are lines defined in the VolSurface :
0832       const std::vector< std::pair<Vector3D, Vector3D> >& local_lines = _volSurf.getLines() ;
0833       
0834       if( local_lines.size() > 0 ) {
0835         unsigned n=local_lines.size() ;
0836         lines.reserve( n ) ;
0837     
0838         for( unsigned i=0;i<n;++i){
0839       
0840           Vector3D av,bv;
0841           _wtM->LocalToMaster( local_lines[i].first ,  av.array() ) ;
0842           _wtM->LocalToMaster( local_lines[i].second , bv.array() ) ;
0843       
0844           lines.emplace_back( av, bv );
0845         }
0846     
0847         return lines ;
0848       }
0849       //--------------------------------------------
0850 
0851 
0852       // get local and global surface vectors
0853       const Vector3D& lu = _volSurf.u() ;
0854       //      const Vector3D& lv = _volSurf.v() ;
0855       const Vector3D& ln = _volSurf.normal() ;
0856       Vector3D lo = _volSurf.origin() ;
0857       
0858       Volume vol = volume() ; 
0859       const TGeoShape* shape = vol->GetShape() ;
0860       
0861       
0862       if( type().isPlane() ) {
0863     
0864         if( shape->IsA() == TGeoBBox::Class() ) {
0865       
0866           TGeoBBox* box = ( TGeoBBox* ) shape  ;
0867       
0868           Vector3D boxDim(  box->GetDX() , box->GetDY() , box->GetDZ() ) ;  
0869       
0870       
0871           bool isYZ = std::fabs(  ln.x() - 1.0 ) < epsilon  ; // normal parallel to x
0872           bool isXZ = std::fabs(  ln.y() - 1.0 ) < epsilon  ; // normal parallel to y
0873           bool isXY = std::fabs(  ln.z() - 1.0 ) < epsilon  ; // normal parallel to z
0874       
0875       
0876           if( isYZ || isXZ || isXY ) {  // plane is parallel to one of the box' sides -> need 4 vertices from box dimensions
0877         
0878             // if isYZ :
0879             unsigned uidx = 1 ;
0880             unsigned vidx = 2 ;
0881         
0882             Vector3D ubl(  0., 1., 0. ) ; 
0883             Vector3D vbl(  0., 0., 1. ) ;
0884         
0885             if( isXZ ) {
0886           
0887               ubl.fill( 1., 0., 0. ) ;
0888               vbl.fill( 0., 0., 1. ) ;
0889               uidx = 0 ;
0890               vidx = 2 ;
0891           
0892             } else if( isXY ) {
0893           
0894               ubl.fill( 1., 0., 0. ) ;
0895               vbl.fill( 0., 1., 0. ) ;
0896               uidx = 0 ;
0897               vidx = 1 ;
0898             }
0899         
0900             Vector3D ub ;
0901             Vector3D vb ;
0902             _wtM->LocalToMasterVect( ubl , ub.array() ) ;
0903             _wtM->LocalToMasterVect( vbl , vb.array() ) ;
0904         
0905             lines.reserve(4) ;
0906         
0907             lines.emplace_back(_o + boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb ,  _o - boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb );
0908             lines.emplace_back(_o - boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb ,  _o - boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb );
0909             lines.emplace_back(_o - boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb ,  _o + boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb );
0910             lines.emplace_back(_o + boxDim[ uidx ] * ub  - boxDim[ vidx ] * vb ,  _o + boxDim[ uidx ] * ub  + boxDim[ vidx ] * vb );
0911         
0912             return lines ;
0913           }     
0914 
0915         }  else if( shape->InheritsFrom("TGeoTube") ||  shape->InheritsFrom("TGeoCone") ) {
0916 
0917 
0918           // can only deal with special case of z-disk and origin in center of cone
0919           if( type().isZDisk() ) { // && lo.rho() < epsilon ) {
0920         
0921             if( lo.rho() > epsilon ) {
0922               // move origin to z-axis 
0923               lo.x() = 0. ; 
0924               lo.y() = 0. ;
0925             }
0926 
0927             double zhalf = 0 ;
0928             double rmax1 = 0 ;
0929             double rmax2 = 0 ;
0930             double rmin1 = 0 ;
0931             double rmin2 = 0 ;
0932             if( shape->InheritsFrom("TGeoTube") ) {
0933               TGeoTube* tube = ( TGeoTube* ) shape  ;
0934               zhalf = tube->GetDZ() ;
0935               rmax1 = tube->GetRmax() ;
0936               rmax2 = tube->GetRmax() ;
0937               rmin1 = tube->GetRmin() ;
0938               rmin2 = tube->GetRmin() ;
0939             } else {   // shape->InheritsFrom("TGeoCone") )
0940               TGeoCone* cone = ( TGeoCone* ) shape  ;
0941               zhalf = cone->GetDZ() ;
0942               rmax1 = cone->GetRmax1() ;
0943               rmax2 = cone->GetRmax2() ;
0944               rmin1 = cone->GetRmin1() ;
0945               rmin2 = cone->GetRmin2() ;
0946             }
0947 
0948             // two circles around origin 
0949             // get radii at position of plane 
0950             double r0 =  rmin1 +  ( rmin2 - rmin1 ) / ( 2. * zhalf )   *  ( zhalf + lo.z()  ) ;  
0951             double r1 =  rmax1 +  ( rmax2 - rmax1 ) / ( 2. * zhalf )   *  ( zhalf + lo.z()  ) ;  
0952 
0953         
0954             unsigned n = nMax / 4 ;
0955             double dPhi = 2.* ROOT::Math::Pi() / double( n ) ; 
0956         
0957             for( unsigned i = 0 ; i < n ; ++i ) {
0958           
0959               Vector3D rv00(  r0*sin(  i   *dPhi ) , r0*cos(  i   *dPhi )  , 0. ) ;
0960               Vector3D rv01(  r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi )  , 0. ) ;
0961 
0962               Vector3D rv10(  r1*sin(  i   *dPhi ) , r1*cos(  i   *dPhi )  , 0. ) ;
0963               Vector3D rv11(  r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi )  , 0. ) ;
0964           
0965          
0966               Vector3D pl0 =  lo + rv00 ;
0967               Vector3D pl1 =  lo + rv01 ;
0968 
0969               Vector3D pl2 =  lo + rv10 ;
0970               Vector3D pl3 =  lo + rv11 ;
0971           
0972 
0973               Vector3D pg0,pg1,pg2,pg3 ;
0974           
0975               _wtM->LocalToMaster( pl0, pg0.array() ) ;
0976               _wtM->LocalToMaster( pl1, pg1.array() ) ;
0977               _wtM->LocalToMaster( pl2, pg2.array() ) ;
0978               _wtM->LocalToMaster( pl3, pg3.array() ) ;
0979           
0980               lines.emplace_back( pg0, pg1 );
0981               lines.emplace_back( pg2, pg3 );
0982             }
0983 
0984             //add some vertical and horizontal lines so that the disc is seen in the rho-z projection
0985 
0986             n = 4 ; dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
0987 
0988             for( unsigned i = 0 ; i < n ; ++i ) {
0989           
0990               Vector3D rv0(  r0*sin( i * dPhi ) , r0*cos(  i * dPhi )  , 0. ) ;
0991               Vector3D rv1(  r1*sin( i * dPhi ) , r1*cos(  i * dPhi )  , 0. ) ;
0992           
0993               Vector3D pl0 =  lo + rv0 ;
0994               Vector3D pl1 =  lo + rv1 ;
0995 
0996               Vector3D pg0,pg1 ;
0997           
0998               _wtM->LocalToMaster( pl0, pg0.array() ) ;
0999               _wtM->LocalToMaster( pl1, pg1.array() ) ;
1000           
1001               lines.emplace_back(pg0, pg1);
1002             }
1003 
1004           }
1005       
1006           return lines ;
1007         }
1008 
1009         else if(shape->IsA() == TGeoTrap::Class()) {
1010           TGeoTrap* trapezoid = ( TGeoTrap* ) shape;
1011           
1012           double dx1 = trapezoid->GetBl1();
1013           double dx2 = trapezoid->GetTl1();
1014           double dz = trapezoid->GetH1();
1015 
1016           //according to the TGeoTrap definition, the lengths are given such that the normal vector of the surface
1017           //points in the e_z direction.
1018           Vector3D ubl(  1., 0., 0. ) ; 
1019           Vector3D vbl(  0., 1., 0. ) ; 
1020 
1021           //the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
1022           Vector3D ub ;
1023           Vector3D vb ;
1024           _wtM->LocalToMasterVect( ubl , ub.array() ) ;
1025           _wtM->LocalToMasterVect( vbl , vb.array() ) ;
1026 
1027           //the trapezoid is drawn as a set of four lines connecting its four corners
1028           lines.reserve(4) ;
1029           //_o is vector to the origin
1030           lines.emplace_back( _o + dx1 * ub  - dz * vb ,  _o + dx2 * ub  + dz * vb);
1031           lines.emplace_back( _o + dx2 * ub  + dz * vb ,  _o - dx2 * ub  + dz * vb);
1032           lines.emplace_back( _o - dx2 * ub  + dz * vb ,  _o - dx1 * ub  - dz * vb);
1033           lines.emplace_back( _o - dx1 * ub  - dz * vb ,  _o + dx1 * ub  - dz * vb);
1034 
1035           return lines;
1036         }
1037         //added code by Thorben Quast for simplified set of lines for trapezoids with unequal lengths in x
1038         else if(shape->IsA() == TGeoTrd1::Class()){
1039           TGeoTrd1* trapezoid = ( TGeoTrd1* ) shape;
1040           //all lengths are half length
1041           double dx1 = trapezoid->GetDx1();
1042           double dx2 = trapezoid->GetDx2();
1043           double dy = trapezoid->GetDy();
1044           double dz = trapezoid->GetDz();
1045 
1046           bool isYZ = std::fabs(  ln.x() - 1.0 ) < epsilon  ; // normal parallel to x
1047           bool isXZ = std::fabs(  ln.y() - 1.0 ) < epsilon  ; // normal parallel to y
1048           bool isXY = std::fabs(  ln.z() - 1.0 ) < epsilon  ; // normal parallel to z
1049       
1050           if(not (isYZ || isXZ || isXY)) {
1051             std::stringstream sst ; 
1052             sst << " ***** ERROR: Trapezoid surface cannot be defined, normal not parallel to x, y, or z axis";
1053             throw std::runtime_error( sst.str() ) ;
1054           }
1055           
1056           Vector3D ubl, vbl;
1057 
1058           if(isYZ) {
1059             ubl.fill( 0., 1., 0. ) ;
1060             vbl.fill( 0., 0., 1. ) ;
1061           } else if( isXZ ) {
1062             ubl.fill( 1., 0., 0. ) ;
1063             vbl.fill( 0., 0., 1. ) ;        
1064           } else if( isXY ) {
1065             ubl.fill( 1., 0., 0. ) ;
1066             vbl.fill( 0., 1., 0. ) ;
1067           } 
1068           
1069           Vector3D ub ;
1070           Vector3D vb ;
1071           _wtM->LocalToMasterVect( ubl , ub.array() ) ;
1072           _wtM->LocalToMasterVect( vbl , vb.array() ) ;
1073           //the trapezoid is drawn as a set of four lines connecting its four corners
1074           lines.reserve(4) ;
1075           //_o is vector to the origin
1076           
1077           if( isYZ ) {
1078             lines.emplace_back( _o + dy * ub  - dz * vb ,  _o + dy * ub  + dz * vb);
1079             lines.emplace_back( _o + dy * ub  + dz * vb ,  _o - dy * ub  + dz * vb);
1080             lines.emplace_back( _o - dy * ub  + dz * vb ,  _o - dy * ub  - dz * vb);
1081             lines.emplace_back( _o - dy * ub  - dz * vb ,  _o + dy * ub  - dz * vb);          
1082           } else if( isXZ ) {
1083             lines.emplace_back( _o + dx1 * ub  - dz * vb ,  _o + dx2 * ub  + dz * vb);
1084             lines.emplace_back( _o + dx2 * ub  + dz * vb ,  _o - dx2 * ub  + dz * vb);
1085             lines.emplace_back( _o - dx2 * ub  + dz * vb ,  _o - dx1 * ub  - dz * vb);
1086             lines.emplace_back( _o - dx1 * ub  - dz * vb ,  _o + dx1 * ub  - dz * vb);        
1087           } else if( isXY ) {
1088             lines.emplace_back( _o + dx1 * ub  - dy * vb ,  _o + dx2 * ub  + dy * vb);
1089             lines.emplace_back( _o + dx2 * ub  + dy * vb ,  _o - dx2 * ub  + dy * vb);
1090             lines.emplace_back( _o - dx2 * ub  + dy * vb ,  _o - dx1 * ub  - dy * vb);
1091             lines.emplace_back( _o - dx1 * ub  - dy * vb ,  _o + dx1 * ub  - dy * vb);        
1092           }
1093           return lines ;
1094         }
1095         //added code by Thorben Quast for simplified set of lines for trapezoids with unequal lengths in x AND y
1096         else if(shape->IsA() == TGeoTrd2::Class()){
1097           TGeoTrd2* trapezoid = ( TGeoTrd2* ) shape;
1098           //all lengths are half length
1099           double dx1 = trapezoid->GetDx1();
1100           double dx2 = trapezoid->GetDx2();
1101           //double dy1 = trapezoid->GetDy1();  
1102           //double dy2 = trapezoid->GetDy2();  
1103           double dz = trapezoid->GetDz();
1104 
1105           //the normal vector is parallel to e_y for all geometry cases in CLIC
1106           //if that is at some point not the case anymore, then local plane vectors ubl, vbl
1107           //must be initialized like it is done for the boxes (line 674 following)
1108           Vector3D ubl(  1., 0., 0. ) ; 
1109           Vector3D vbl(  0., 0., 1. ) ; 
1110           
1111           //the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
1112           Vector3D ub ;
1113           Vector3D vb ;
1114           _wtM->LocalToMasterVect( ubl , ub.array() ) ;
1115           _wtM->LocalToMasterVect( vbl , vb.array() ) ;
1116 
1117           //the trapezoid is drawn as a set of four lines connecting its four corners
1118           lines.reserve(4) ;
1119           //_o is vector to the origin
1120           lines.emplace_back( _o + dx1 * ub  - dz * vb ,  _o + dx2 * ub  + dz * vb);
1121           lines.emplace_back( _o + dx2 * ub  + dz * vb ,  _o - dx2 * ub  + dz * vb);
1122           lines.emplace_back( _o - dx2 * ub  + dz * vb ,  _o - dx1 * ub  - dz * vb);
1123           lines.emplace_back( _o - dx1 * ub  - dz * vb ,  _o + dx1 * ub  - dz * vb);
1124 
1125           return lines;
1126         }
1127         // ===== default for arbitrary planes in arbitrary shapes ================= 
1128     
1129         // We create nMax vertices by rotating the local u vector around the normal
1130         // and checking the distance to the volume boundary in that direction.
1131         // This is brute force and not very smart, as many points are created on straight 
1132         // lines and the edges are still rounded. 
1133         // The alterative would be to compute the true intersections a plane and the most
1134         // common shapes - at least for boxes that should be not too hard. To be done...
1135     
1136         lines.reserve( nMax ) ;
1137     
1138         double dAlpha =  2.* ROOT::Math::Pi() / double( nMax ) ; 
1139 
1140         TVector3 norm( ln.x() , ln.y() , ln.z() ) ;
1141 
1142     
1143         Vector3D first, previous ;
1144 
1145         for(unsigned i=0 ; i< nMax ; ++i ){ 
1146       
1147           double alpha = double(i) * dAlpha ;
1148       
1149           TVector3 vec( lu.x() , lu.y() , lu.z() ) ;
1150 
1151           TRotation rot ;
1152           rot.Rotate( alpha , norm );
1153     
1154           TVector3 vecR = rot * vec ;
1155       
1156           Vector3D luRot ;
1157           luRot.fill( vecR ) ;
1158       
1159           double dist = shape->DistFromInside( const_cast<double*> (lo.const_array()) , const_cast<double*> (luRot.const_array())  , 3, 0.1 ) ;
1160       
1161           // local point at volume boundary
1162           Vector3D lp = lo + dist * luRot ;
1163 
1164           Vector3D gp ;
1165       
1166           _wtM->LocalToMaster( lp , gp.array() ) ;
1167 
1168           // std::cout << " **** normal:" << ln << " lu:" << lu  << " alpha:" << alpha << " luRot:" << luRot << " lp :" << lp  << " gp:" << gp << " dist : " << dist 
1169           //        << " is point " << gp << " inside : " << shape->Contains( gp.const_array()  )  
1170           //        << " dist from outside for lo,lu " <<  shape->DistFromOutside( lo.const_array()  , lu.const_array()   , 3 )    
1171           //        << " dist from inside for lo,ln " <<  shape->DistFromInside( lo.const_array()  , ln.const_array()   , 3 )    
1172           //        << std::endl;
1173           //      shape->Dump() ;
1174       
1175 
1176           if( i >  0 ) 
1177             lines.emplace_back(previous, gp);
1178           else
1179             first = gp ;
1180 
1181           previous = gp ;
1182         }
1183         lines.emplace_back(previous, first);
1184 
1185 
1186       } else if( type().isCylinder() ) {  
1187 
1188         if( shape->InheritsFrom("TGeoTube") || shape->InheritsFrom("TGeoCone") ) {
1189 
1190           lines.reserve( nMax ) ;
1191 
1192           TGeoBBox* tube = ( TGeoBBox* ) shape  ;
1193 
1194           double zHalf = tube->GetDZ() ;
1195 
1196           Vector3D zv( 0. , 0. , zHalf ) ;
1197       
1198           double r = lo.rho() ;
1199 
1200 
1201           unsigned n = nMax / 4 ;
1202           double dPhi = 2.* ROOT::Math::Pi() / double( n ) ; 
1203 
1204           for( unsigned i = 0 ; i < n ; ++i ) {
1205 
1206             Vector3D rv0(  r*sin(  i   *dPhi ) , r*cos(  i   *dPhi )  , 0. ) ;
1207             Vector3D rv1(  r*sin( (i+1)*dPhi ) , r*cos( (i+1)*dPhi )  , 0. ) ;
1208 
1209             // 4 points on local cylinder
1210 
1211             Vector3D pl0 =  zv + rv0 ;
1212             Vector3D pl1 =  zv + rv1 ;
1213             Vector3D pl2 = -zv + rv1  ;
1214             Vector3D pl3 = -zv + rv0 ;
1215 
1216             Vector3D pg0,pg1,pg2,pg3 ;
1217 
1218             _wtM->LocalToMaster( pl0, pg0.array() ) ;
1219             _wtM->LocalToMaster( pl1, pg1.array() ) ;
1220             _wtM->LocalToMaster( pl2, pg2.array() ) ;
1221             _wtM->LocalToMaster( pl3, pg3.array() ) ;
1222 
1223             lines.emplace_back( pg0, pg1 );
1224             lines.emplace_back( pg1, pg2 );
1225             lines.emplace_back( pg2, pg3 );
1226             lines.emplace_back( pg3, pg0 );
1227           }
1228         }
1229       }
1230       return lines ;
1231 
1232     }
1233 
1234     //================================================================================================================
1235 
1236  
1237     Vector3D CylinderSurface::u( const Vector3D& point  ) const { 
1238  
1239       Vector3D lp , u_val ;
1240       _wtM->MasterToLocal( point , lp.array() ) ;
1241       const Vector3D& lu = _volSurf.u( lp  ) ;
1242       _wtM->LocalToMasterVect( lu , u_val.array() ) ;
1243       return u_val ; 
1244     }
1245     
1246     Vector3D CylinderSurface::v(const Vector3D& point ) const {  
1247       Vector3D lp , v_val ;
1248       _wtM->MasterToLocal( point , lp.array() ) ;
1249       const Vector3D& lv =  _volSurf.v( lp  ) ;
1250       _wtM->LocalToMasterVect( lv , v_val.array() ) ;
1251       return v_val ; 
1252     }
1253     
1254     Vector3D CylinderSurface::normal(const Vector3D& point ) const {  
1255       Vector3D lp , n ;
1256       _wtM->MasterToLocal( point , lp.array() ) ;
1257       const Vector3D& ln =  _volSurf.normal( lp  ) ;
1258       _wtM->LocalToMasterVect( ln , n.array() ) ;
1259       return n ; 
1260     }
1261  
1262     Vector2D CylinderSurface::globalToLocal( const Vector3D& point) const {
1263       
1264       Vector3D lp;
1265       _wtM->MasterToLocal( point , lp.array() ) ;
1266  
1267       return _volSurf.globalToLocal( lp )  ;
1268     }
1269     
1270     
1271     Vector3D CylinderSurface::localToGlobal( const Vector2D& point) const {
1272  
1273       Vector3D lp = _volSurf.localToGlobal( point ) ;
1274       Vector3D p ;
1275       _wtM->LocalToMaster( lp , p.array() ) ;
1276  
1277       return p ;
1278     }
1279 
1280     double CylinderSurface::radius() const {    return _volSurf.origin().rho() ;  }
1281 
1282     Vector3D CylinderSurface::center() const {  return volumeOrigin() ;  }
1283 
1284 
1285     //================================================================================================================
1286 
1287 
1288     double ConeSurface::radius0() const {
1289       
1290       double theta = _volSurf.v().theta() ;
1291       double l = length_along_v() * cos( theta ) ;
1292       
1293       return origin().rho() - 0.5 * l * tan( theta ) ;  
1294     }
1295 
1296     double ConeSurface::radius1() const {
1297       
1298       double theta = _volSurf.v().theta() ;
1299       double l = length_along_v() * cos( theta ) ;
1300       
1301       return origin().rho() + 0.5 * l * tan( theta ) ;  
1302     }
1303 
1304     double ConeSurface::z0() const {
1305       
1306       double theta = _volSurf.v().theta() ;
1307       double l = length_along_v() * cos( theta ) ;
1308       
1309       return origin().z() - 0.5 * l ;  
1310     }
1311 
1312     double ConeSurface::z1() const {
1313       
1314       double theta = _volSurf.v().theta() ;
1315       double l = length_along_v() * cos( theta ) ;
1316       
1317       return origin().z() + 0.5 * l ;
1318     }
1319 
1320     Vector3D ConeSurface::center() const {  return volumeOrigin() ;  }
1321 
1322 
1323     //================================================================================================================
1324 
1325   } // namespace
1326 } // namespace
1327 
1328
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!