DDRec/src/Surface.cpp

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