DD4hep - The AIDA detector description toolkit for high energy physics experiments
DD4hep  Rev:Unversioneddirectory
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Surface.cpp
Go to the documentation of this file.
1 #include "DDRec/Surface.h"
3 #include "DD4hep/Memory.h"
4 
6 
7 #include <cmath>
8 #include <memory>
9 #include <exception>
10 
11 #include "TGeoMatrix.h"
12 #include "TGeoShape.h"
13 #include "TRotation.h"
14 //TGeoTrd1 is apparently not included by defautl
15 #include "TGeoTrd1.h"
16 
17 namespace DD4hep {
18  namespace DDRec {
19 
20  using namespace Geometry ;
21 
22 
23  //======================================================================================================
24 
25  void VolSurfaceBase::setU(const Vector3D& u_val) { _u = u_val ; }
26  void VolSurfaceBase::setV(const Vector3D& v_val) { _v = v_val ; }
27  void VolSurfaceBase::setNormal(const Vector3D& n) { _n = n ; }
28  void VolSurfaceBase::setOrigin(const Vector3D& o) { _o = o ; }
29 
30  long64 VolSurfaceBase::id() const { return _id ; }
31 
32  const SurfaceType& VolSurfaceBase::type() const { return _type ; }
33  Vector3D VolSurfaceBase::u( const Vector3D& point ) const { point.x() ; return _u ; }
34  Vector3D VolSurfaceBase::v(const Vector3D& point ) const { point.x() ; return _v ; }
35  Vector3D VolSurfaceBase::normal(const Vector3D& point ) const { point.x() ; return _n ; }
36  const Vector3D& VolSurfaceBase::origin() const { return _o ;}
37 
39 
40  Vector3D p = point - origin() ;
41 
42  // create new orthogonal unit vectors
43  // FIXME: these vectors should be cached really ...
44 
45  double uv = u() * v() ;
46  Vector3D uprime = ( u() - uv * v() ).unit() ;
47  Vector3D vprime = ( v() - uv * u() ).unit() ;
48  double uup = u() * uprime ;
49  double vvp = v() * vprime ;
50 
51  return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
52  }
53 
55 
56  Vector3D g = origin() + point[0] * u() + point[1] * v() ;
57 
58  return g ;
59  }
60 
61  const IMaterial& VolSurfaceBase::innerMaterial() const{ return _innerMat ; }
62  const IMaterial& VolSurfaceBase::outerMaterial() const { return _outerMat ; }
63  double VolSurfaceBase::innerThickness() const { return _th_i ; }
64  double VolSurfaceBase::outerThickness() const { return _th_o ; }
65 
66 
68 
69  const DDSurfaces::Vector3D& o = this->origin() ;
70  const DDSurfaces::Vector3D& u_val = this->u( o ) ;
71  DDSurfaces::Vector3D um = -1. * u_val ;
72 
73  double dist_p = 0. ;
74  double dist_m = 0. ;
75 
76 
77  // std::cout << " VolSurfaceBase::length_along_u() : o = " << o << " u = " << this->u( o )
78  // << " -u = " << um << std::endl ;
79 
80 
81  if( volume()->GetShape()->Contains( o.const_array() ) ){
82 
83  dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
84  const_cast<double*> ( u_val.const_array() ) ) ;
85  dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
86  const_cast<double*> ( um.array() ) ) ;
87 
88 
89  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
90  // << " dist_p " << dist_p
91  // << " dist_m " << dist_m
92  // << std::endl ;
93 
94 
95  } else{
96 
97  dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
98  const_cast<double*> ( u_val.const_array() ) ) ;
99  dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
100  const_cast<double*> ( um.array() ) ) ;
101 
102  dist_p *= 1.0001 ;
103  dist_m *= 1.0001 ;
104 
105  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
106  // << " dist_p " << dist_p
107  // << " dist_m " << dist_m
108  // << std::endl ;
109 
110  DDSurfaces::Vector3D o_1 = this->origin() + dist_p * u_val ;
111  DDSurfaces::Vector3D o_2 = this->origin() + dist_m * um ;
112 
113  dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) ,
114  const_cast<double*> ( u_val.const_array() ) ) ;
115 
116  dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) ,
117  const_cast<double*> ( um.array() ) ) ;
118 
119  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
120  // << " dist_p " << dist_p
121  // << " dist_m " << dist_m
122  // << std::endl ;
123  }
124 
125  return dist_p + dist_m ;
126 
127 
128  }
129 
131 
132  const DDSurfaces::Vector3D& o = this->origin() ;
133  const DDSurfaces::Vector3D& v_val = this->v( o ) ;
134  DDSurfaces::Vector3D vm = -1. * v_val ;
135 
136  double dist_p = 0. ;
137  double dist_m = 0. ;
138 
139 
140  // std::cout << " VolSurfaceBase::length_along_u() : o = " << o << " u = " << this->u( o )
141  // << " -u = " << vm << std::endl ;
142 
143 
144  if( volume()->GetShape()->Contains( o.const_array() ) ){
145 
146  dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
147  const_cast<double*> ( v_val.const_array() ) ) ;
148  dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
149  const_cast<double*> ( vm.array() ) ) ;
150 
151 
152  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
153  // << " dist_p " << dist_p
154  // << " dist_m " << dist_m
155  // << std::endl ;
156 
157 
158  } else{
159 
160  dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
161  const_cast<double*> ( v_val.const_array() ) ) ;
162  dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
163  const_cast<double*> ( vm.array() ) ) ;
164 
165  dist_p *= 1.0001 ;
166  dist_m *= 1.0001 ;
167 
168  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
169  // << " dist_p " << dist_p
170  // << " dist_m " << dist_m
171  // << std::endl ;
172 
173  DDSurfaces::Vector3D o_1 = this->origin() + dist_p * v_val ;
174  DDSurfaces::Vector3D o_2 = this->origin() + dist_m * vm ;
175 
176  dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) ,
177  const_cast<double*> ( v_val.const_array() ) ) ;
178 
179  dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) ,
180  const_cast<double*> ( vm.array() ) ) ;
181 
182  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
183  // << " dist_p " << dist_p
184  // << " dist_m " << dist_m
185  // << std::endl ;
186  }
187 
188  return dist_p + dist_m ;
189 
190  }
191 
192 
193  double VolSurfaceBase::distance(const Vector3D& point ) const { point.x() ; return 1.e99 ; }
194 
196  bool VolSurfaceBase::insideBounds(const Vector3D& point, double epsilon) const {
197 
198 #if 0
199  double dist = std::abs ( distance( point ) ) ;
200 
201  bool inShape = volume()->GetShape()->Contains( point.const_array() ) ;
202 
203  std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << dist
204  << " origin = " << origin() << " normal = " << normal()
205  << " p * n = " << point * normal()
206  << " isInShape : " << inShape << std::endl ;
207 
208  return dist < epsilon && inShape ;
209 #else
210 
211  //fixme: older versions of ROOT (~<5.34.10 ) take a non const pointer as argument - therefore use a const cast here for the time being ...
212  return ( std::abs ( distance( point ) ) < epsilon ) && volume()->GetShape()->Contains( const_cast<double*> (point.const_array() ) ) ;
213 #endif
214 
215  }
216 
217 
218  std::vector< std::pair<Vector3D, Vector3D> > VolSurfaceBase::getLines(unsigned ) {
219  // dummy implementation returning empty set
220  std::vector< std::pair<Vector3D, Vector3D> > lines ;
221  return lines ;
222  }
223 
224  //===================================================================
225  // simple wrapper methods forwarding the call to the implementation object
226 
227  long64 VolSurface::id() const { return _surf->id() ; }
228  const SurfaceType& VolSurface::type() const { return _surf->type() ; }
229  Vector3D VolSurface::u( const Vector3D& point ) const { return _surf->u(point) ; }
230  Vector3D VolSurface::v(const Vector3D& point ) const { return _surf->v(point) ; }
231  Vector3D VolSurface::normal(const Vector3D& point ) const { return _surf->normal(point) ; }
232  const Vector3D& VolSurface::origin() const { return _surf->origin() ;}
233  Vector2D VolSurface::globalToLocal( const Vector3D& point) const { return _surf->globalToLocal( point ) ; }
234  Vector3D VolSurface::localToGlobal( const Vector2D& point) const { return _surf->localToGlobal( point) ; }
235  const IMaterial& VolSurface::innerMaterial() const{ return _surf->innerMaterial() ; }
236  const IMaterial& VolSurface::outerMaterial() const { return _surf->outerMaterial() ; }
237  double VolSurface::innerThickness() const { return _surf->innerThickness() ; }
238  double VolSurface::outerThickness() const { return _surf->outerThickness() ; }
239  double VolSurface::length_along_u() const { return _surf->length_along_u() ; }
240  double VolSurface::length_along_v() const { return _surf->length_along_v() ; }
241  double VolSurface::distance(const Vector3D& point ) const { return _surf->distance( point ) ; }
242  bool VolSurface::insideBounds(const Vector3D& point, double epsilon) const {
243  return _surf->insideBounds( point, epsilon ) ;
244  }
245  std::vector< std::pair<Vector3D, Vector3D> > VolSurface::getLines(unsigned nMax) {
246  return _surf->getLines(nMax) ;
247  }
248 
249  //===================================================================
250 
252  double VolPlaneImpl::distance(const Vector3D& point ) const {
253  return ( point - origin() ) * normal() ;
254  }
255  //======================================================================================================
256 
258  double thickness_inner ,double thickness_outer, Vector3D o ) :
259 
260  VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , Vector3D() , Vector3D() , o , vol, 0) {
261  Vector3D v_val( 0., 0., 1. ) ;
262  Vector3D o_rphi( o.x() , o.y() , 0. ) ;
263  Vector3D n = o_rphi.unit() ;
264  Vector3D u_val = v_val.cross( n ) ;
265 
266  setU( u_val ) ;
267  setV( v_val ) ;
268  setNormal( n ) ;
269 
272  _type.setProperty( SurfaceType::Cone , false ) ;
273  _type.checkParallelToZ( *this ) ;
274  _type.checkOrthogonalToZ( *this ) ;
275  }
276 
277  Vector3D VolCylinderImpl::u(const Vector3D& point ) const {
278 
279  Vector3D n( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
280 
281  return v().cross( n ) ;
282  }
283 
284  Vector3D VolCylinderImpl::normal(const Vector3D& point ) const {
285 
286  // normal is just given by phi of the point
287  return Vector3D( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
288  }
289 
291 
292  // cylinder is parallel to v here so u is Z and v is r *Phi
293  double phi = point.phi() - origin().phi() ;
294 
295  while( phi < -M_PI ) phi += 2.*M_PI ;
296  while( phi > M_PI ) phi -= 2.*M_PI ;
297 
298  return Vector2D( origin().rho() * phi, point.z() - origin().z() ) ;
299  }
300 
301 
303 
304  double z = point.v() + origin().z() ;
305  double phi = point.u() / origin().rho() + origin().phi() ;
306 
307  while( phi < -M_PI ) phi += 2.*M_PI ;
308  while( phi > M_PI ) phi -= 2.*M_PI ;
309 
310  return Vector3D( origin().rho() , phi, z , Vector3D::cylindrical ) ;
311  }
312 
313 
315  double VolCylinderImpl::distance(const Vector3D& point ) const {
316 
317  return point.rho() - origin().rho() ;
318  }
319 
320  //================================================================================================================
322  double thickness_inner ,double thickness_outer, Vector3D v_val, Vector3D o_val ) :
323 
324  VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , v_val , Vector3D() , Vector3D() , vol, 0) {
325 
326  Vector3D o_rphi( o_val.x() , o_val.y() , 0. ) ;
327 
328  // sanity check: v and o have to have a common phi
329  double dphi = v_val.phi() - o_rphi.phi() ;
330  while( dphi < -M_PI ) dphi += 2.*M_PI ;
331  while( dphi > M_PI ) dphi -= 2.*M_PI ;
332 
333  if( std::fabs( dphi ) > 1e-6 ){
334  std::stringstream sst ; sst << "VolConeImpl::VolConeImpl() - incompatibel vector v and o given "
335  << v_val << " - " << o_val ;
336  throw std::runtime_error( sst.str() ) ;
337  }
338 
339  double theta = v_val.theta() ;
340 
341  Vector3D n( 1. , v_val.phi() , ( theta + M_PI/2. ) , Vector3D::spherical ) ;
342  Vector3D u_val = v_val.cross( n ) ;
343 
344  setU( u_val ) ;
345  setOrigin( o_rphi ) ;
346  setNormal( n ) ;
347 
348  // set helper variable for faster computations (describe cone with tip at origin)
349  _tanTheta = std::tan( theta ) ;
350  double tipoffset = o_val.rho() / _tanTheta ; // distance from tip to origin.z()
351  _ztip = o_val.z() - tipoffset ;
352 
353  double dist_p = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) ,
354  const_cast<double*> ( v_val.const_array() ) ) ;
355  Vector3D vm = -1. * v_val ;
356  double dist_m = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) ,
357  const_cast<double*> ( vm.array() ) ) ;
358 
359  double costh = std::cos( theta) ;
360  _zt0 = tipoffset - dist_m * costh ;
361  _zt1 = tipoffset + dist_p * costh ;
362 
363 
369  }
370 
371 
372  Vector3D VolConeImpl::v(const Vector3D& point ) const {
373  // just take phi from point
374  Vector3D av( 1. , point.phi() , _v.theta() , Vector3D::spherical ) ;
375  return av ;
376  }
377 
378  Vector3D VolConeImpl::u(const Vector3D& point ) const {
379  // compute from v X n
380  const Vector3D& av = this->v( point ) ;
381  const Vector3D& n = normal( point ) ;
382  return av.cross( n ) ;
383  }
384 
385  Vector3D VolConeImpl::normal(const Vector3D& point ) const {
386  // just take phi from point
387  Vector3D n( 1. , point.phi() , _n.theta() , Vector3D::spherical ) ;
388  return n ;
389  }
390 
392 
393  // cone is parallel to z here, so u is r *Phi and v is "along" z
394  double phi = point.phi() - origin().phi() ;
395 
396  while( phi < -M_PI ) phi += 2.*M_PI ;
397  while( phi > M_PI ) phi -= 2.*M_PI ;
398 
399 
400  double r = ( point.z() - _ztip ) * _tanTheta ;
401 
402  return Vector2D( r*phi, ( point.z() - origin().z() ) / cos( _v.theta() ) ) ;
403  }
404 
405 
407 
408  double z = point.v() * cos( _v.theta() ) + origin().z() ;
409 
410  double r = ( z - _ztip ) * _tanTheta ;
411 
412  double phi = point.u() / r + origin().phi() ;
413 
414  while( phi < -M_PI ) phi += 2.*M_PI ;
415  while( phi > M_PI ) phi -= 2.*M_PI ;
416 
417  return Vector3D( r , phi, z , Vector3D::cylindrical ) ;
418  }
419 
420 
422  double VolConeImpl::distance(const Vector3D& point ) const {
423 
424  // // if the point is in the other hemispere we return the distance to origin
425  // // -> this assumes that the cones do not cross the xy-plane ...
426  // // otherwise we get the distance to the mirrored part of the cone
427  // // needs more thought ..
428  // if( origin().z() * point.z() < 0. )
429  // return point.r() ;
430 
431  //fixme: there are probably faster ways to compute this
432  // e.g by using the intercept theorem - tbd. ...
433  // const Vector2D& lp = globalToLocal( point ) ;
434  // const Vector3D& gp = localToGlobal( lp ) ;
435 
436  // Vector3D dz = point - gp ;
437 
438  //return dz * normal( point ) ;
439 
440  double zp = point.z() - _ztip ;
441  double r = point.rho() - zp * _tanTheta ;
442  return r * std::cos( _v.theta() ) ;
443 
444  }
445 
447  std::vector< std::pair<DDSurfaces::Vector3D, DDSurfaces::Vector3D> > VolConeImpl::getLines(unsigned nMax){
448 
449  std::vector< std::pair<DDSurfaces::Vector3D, DDSurfaces::Vector3D> > lines ;
450 
451  lines.reserve( nMax ) ;
452 
453  double theta = v().theta() ;
454  double half_length = 0.5 * length_along_v() * cos( theta ) ;
455 
456  Vector3D zv( 0. , 0. , half_length ) ;
457 
458  double dr = half_length * tan( theta ) ;
459 
460  double r0 = origin().rho() - dr ;
461  double r1 = origin().rho() + dr ;
462 
463 
464  unsigned n = nMax / 4 ;
465  double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
466 
467  for( unsigned i = 0 ; i < n ; ++i ) {
468 
469  Vector3D r0v0( r0*sin( i *dPhi ) , r0*cos( i *dPhi ) , 0. ) ;
470  Vector3D r0v1( r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi ) , 0. ) ;
471 
472  Vector3D r1v0( r1*sin( i *dPhi ) , r1*cos( i *dPhi ) , 0. ) ;
473  Vector3D r1v1( r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi ) , 0. ) ;
474 
475  Vector3D pl0 = zv + r1v0 ;
476  Vector3D pl1 = zv + r1v1 ;
477  Vector3D pl2 = -zv + r0v1 ;
478  Vector3D pl3 = -zv + r0v0 ;
479 
480  lines.push_back( std::make_pair( pl0, pl1 ) ) ;
481  lines.push_back( std::make_pair( pl1, pl2 ) ) ;
482  lines.push_back( std::make_pair( pl2, pl3 ) ) ;
483  lines.push_back( std::make_pair( pl3, pl0 ) ) ;
484  }
485  return lines;
486  }
487 
488  //================================================================================================================
489 
490 
492  // // delete all surfaces attached to this volume
493  // for( VolSurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
494  // i->clear() ;
495  // }
496  }
497  //=======================================================
498 
500 
501  if( _isOwner ) {
502  // delete all surfaces attached to this volume
503  for( SurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
504  delete (*i) ;
505  }
506  }
507  }
508 
509  //================================================================================================================
510 
512 
513 
514  VolSurfaceList* list = 0 ;
515 
516  try {
517 
518  list = det.extension< VolSurfaceList >() ;
519 
520  } catch(const std::runtime_error& e){
521 
522  list = det.addExtension<VolSurfaceList >( new VolSurfaceList ) ;
523  }
524 
525  return list ;
526  }
527 
528 
529  //======================================================================================================================
530 
531  bool findVolume( PlacedVolume pv, Volume theVol, std::list< PlacedVolume >& volList ) {
532 
533 
534  volList.push_back( pv ) ;
535 
536  // unsigned count = volList.size() ;
537  // for(unsigned i=0 ; i < count ; ++i) {
538  // std::cout << " **" ;
539  // }
540  // std::cout << " searching for volume: " << theVol.name() << " " << std::hex << theVol.ptr() << " <-> pv.volume : " << pv.name() << " " << pv.volume().ptr()
541  // << " pv.volume().ptr() == theVol.ptr() " << (pv.volume().ptr() == theVol.ptr() )
542  // << std::endl ;
543 
544 
545  if( pv.volume().ptr() == theVol.ptr() ) {
546 
547  return true ;
548 
549  } else {
550 
551  //--------------------------------
552 
553  const TGeoNode* node = pv.ptr();
554 
555  if ( !node ) {
556 
557  // std::cout << " *** findVolume: Invalid placement: - node pointer Null for volume: " << pv.name() << std::endl ;
558 
559  throw std::runtime_error("*** findVolume: Invalid placement: - node pointer Null ! " + std::string( pv.name() ) );
560  }
561  // Volume vol = pv.volume();
562 
563  // std::cout << " ndau = " << node->GetNdaughters() << std::endl ;
564 
565  for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
566 
567  TGeoNode* daughter = node->GetDaughter(idau);
568  PlacedVolume placement( daughter );
569 
570  if ( !placement.data() ) {
571  throw std::runtime_error("*** findVolume: Invalid not instrumented placement:"+std::string(daughter->GetName())
572  +" [Internal error -- bad detector constructor]");
573  }
574 
575  PlacedVolume pv_dau = Ref_t(daughter); // why use a Ref_t here ???
576 
577  if( findVolume( pv_dau , theVol , volList ) ) {
578 
579  // std::cout << " ----- found in daughter volume !!! " << std::hex << pv_dau.volume().ptr() << std::endl ;
580 
581  return true ;
582  }
583  }
584 
585  // ------- not found:
586 
587  volList.pop_back() ;
588 
589  return false ;
590  //--------------------------------
591 
592  }
593  }
594 
595  //======================================================================================================================
596 
597  Surface::Surface( Geometry::DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ),
598  _wtM(0) , _id( 0) , _type( _volSurf.type() ) {
599 
600  initialize() ;
601  }
602 
603  long64 Surface::id() const { return _id ; }
604 
605  const SurfaceType& Surface::type() const { return _type ; }
606 
607  Vector3D Surface::u( const Vector3D& point ) const { point.x() ; return _u ; }
608  Vector3D Surface::v(const Vector3D& point ) const { point.x() ; return _v ; }
609  Vector3D Surface::normal(const Vector3D& point ) const { point.x() ; return _n ; }
610  const Vector3D& Surface::origin() const { return _o ;}
611  double Surface::innerThickness() const { return _volSurf.innerThickness() ; }
612  double Surface::outerThickness() const { return _volSurf.outerThickness() ; }
613  double Surface::length_along_u() const { return _volSurf.length_along_u() ; }
614  double Surface::length_along_v() const { return _volSurf.length_along_v() ; }
615 
619 
620  const IMaterial& mat = _volSurf.innerMaterial() ;
621 
622  if( ! ( mat.Z() > 0 ) ) {
623 
624  MaterialManager matMgr ;
625 
626  Vector3D p = _o - innerThickness() * _n ;
627 
628  const MaterialVec& materials = matMgr.materialsBetween( _o , p ) ;
629 
630  _volSurf.setInnerMaterial( materials.size() > 1 ?
631  matMgr.createAveragedMaterial( materials ) :
632  materials[0].first ) ;
633  }
634  return mat ;
635  }
636 
638 
639  const IMaterial& mat = _volSurf.outerMaterial() ;
640 
641  if( ! ( mat.Z() > 0 ) ) {
642 
643  MaterialManager matMgr ;
644 
645  Vector3D p = _o + outerThickness() * _n ;
646 
647  const MaterialVec& materials = matMgr.materialsBetween( _o , p ) ;
648 
649  _volSurf.setOuterMaterial( materials.size() > 1 ?
650  matMgr.createAveragedMaterial( materials ) :
651  materials[0].first ) ;
652  }
653  return mat ;
654  }
655 
656 
657  Vector2D Surface::globalToLocal( const Vector3D& point) const {
658 
659  Vector3D p = point - origin() ;
660 
661  // create new orthogonal unit vectors
662  // FIXME: these vectors should be cached really ...
663 
664  double uv = u() * v() ;
665  Vector3D uprime = ( u() - uv * v() ).unit() ;
666  Vector3D vprime = ( v() - uv * u() ).unit() ;
667  double uup = u() * uprime ;
668  double vvp = v() * vprime ;
669 
670  return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
671  }
672 
673 
674  Vector3D Surface::localToGlobal( const Vector2D& point) const {
675 
676  Vector3D g = origin() + point[0] * u() + point[1] * v() ;
677  return g ;
678  }
679 
680 
682 
683  double o_array[3] ;
684 
685  _wtM->LocalToMaster ( Vector3D() , o_array ) ;
686 
687  Vector3D o(o_array) ;
688 
689  return o ;
690  }
691 
692 
693  double Surface::distance(const Vector3D& point ) const {
694 
695  double pa[3] ;
696  _wtM->MasterToLocal( point , pa ) ;
697  Vector3D localPoint( pa ) ;
698 
699  return _volSurf.distance( localPoint ) ;
700  }
701 
702  bool Surface::insideBounds(const Vector3D& point, double epsilon) const {
703 
704  double pa[3] ;
705  _wtM->MasterToLocal( point , pa ) ;
706  Vector3D localPoint( pa ) ;
707 
708  return _volSurf.insideBounds( localPoint , epsilon) ;
709  }
710 
712 
713  // first we need to find the right volume for the local surface in the DetElement's volumes
714  std::list< PlacedVolume > pVList ;
715  PlacedVolume pv = _det.placement() ;
716  Volume theVol = _volSurf.volume() ;
717 
718  if( ! findVolume( pv, theVol , pVList ) ){
719  theVol = _volSurf.volume() ;
720  std::stringstream sst ; sst << " ***** ERROR: Volume " << theVol.name() << " not found for DetElement " << _det.name() << " with surface " ;
721  throw std::runtime_error( sst.str() ) ;
722  }
723 
724  //=========== compute and cache world transform for surface ==========
725 
726  const TGeoHMatrix& wm = _det.worldTransformation() ;
727 
728 #if 0 // debug
729  wm.Print() ;
730  for( std::list<PlacedVolume>::iterator it= pVList.begin(), n = pVList.end() ; it != n ; ++it ){
731  PlacedVolume pv = *it ;
732  TGeoMatrix* m = pv->GetMatrix();
733  std::cout << " +++ matrix for placed volume : " << std::endl ;
734  m->Print() ;
735  }
736 #endif
737 
738  // need to get the inverse transformation ( see Detector.cpp )
739  // std::auto_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm.Inverse() ) ) ;
740  // has been fixed now, no need to get the inverse anymore:
741  dd4hep_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm ) ) ;
742 
743  //---- if the volSurface is not in the DetElement's volume, we need to mutliply the path to the volume to the
744  // DetElements world transform
745  for( std::list<PlacedVolume>::iterator it = ++( pVList.begin() ) , n = pVList.end() ; it != n ; ++it ){
746 
747  PlacedVolume pvol = *it ;
748  TGeoMatrix* m = pvol->GetMatrix();
749  // std::cout << " +++ matrix for placed volume : " << std::endl ;
750  // m->Print() ;
751  //wtI->MultiplyLeft( m );
752 
753  wtI->Multiply( m );
754  }
755 
756  // std::cout << " +++ new world transform matrix : " << std::endl ;
757 
758 #if 0 //fixme: which convention to use here - the correct should be wtI, however it is the inverse of what is stored in DetElement ???
759  dd4hep_ptr<TGeoHMatrix> wt( new TGeoHMatrix( wtI->Inverse() ) ) ;
760  wt->Print() ;
761  // cache the world transform for the surface
762  _wtM = wt.release() ;
763 #else
764  // wtI->Print() ;
765  // cache the world transform for the surface
766  _wtM = wtI.release() ;
767 #endif
768 
769 
770  // ============ now fill the global surface vectors ==========================
771 
772  double ua[3], va[3], na[3], oa[3] ;
773 
774  _wtM->LocalToMasterVect( _volSurf.u() , ua ) ;
775  _wtM->LocalToMasterVect( _volSurf.v() , va ) ;
776  _wtM->LocalToMasterVect( _volSurf.normal() , na ) ;
777  _wtM->LocalToMaster ( _volSurf.origin() , oa ) ;
778 
779  _u.fill( ua ) ;
780  _v.fill( va ) ;
781  _n.fill( na ) ;
782  _o.fill( oa ) ;
783 
784  // std::cout << " --- local and global surface vectors : ------- " << std::endl
785  // << " u : " << _volSurf.u() << " - " << _u << std::endl
786  // << " v : " << _volSurf.v() << " - " << _v << std::endl
787  // << " n : " << _volSurf.normal() << " - " << _n << std::endl
788  // << " o : " << _volSurf.origin() << " - " << _o << std::endl ;
789 
790 
791  // =========== check parallel and orthogonal to Z ===================
792 
793  if( ! _type.isCone() ) {
794  //fixme: workaround for conical surfaces that should always be parallel to z
795  // however the check with the normal does not work here ...
796 
797  _type.checkParallelToZ( *this ) ;
798 
799  _type.checkOrthogonalToZ( *this ) ;
800  }
801 
802  //======== set the unique surface ID from the DetElement ( and placements below ? )
803 
804  // just use the DetElement ID for now ...
805  // or the id set by the user to the VolSurface ...
806  _id = ( _volSurf.id()==0 ? _det.volumeID() : _volSurf.id() ) ;
807 
808  // typedef PlacedVolume::VolIDs IDV ;
809  // DetElement d = _det ;
810  // while( d.isValid() && d.parent().isValid() ){
811  // PlacedVolume pv = d.placement() ;
812  // if( pv.isValid() ){
813  // const IDV& idV = pv.volIDs() ;
814  // std::cout << " VolIDs : " << d.name() << std::endl ;
815  // for( unsigned i=0, n=idV.size() ; i<n ; ++i){
816  // std::cout << " " << idV[i].first << " - " << idV[i].second << std::endl ;
817  // }
818  // }
819  // d = d.parent() ;
820  // }
821 
822  }
823  //===================================================================================================================
824 
825  std::vector< std::pair<Vector3D, Vector3D> > Surface::getLines(unsigned nMax) {
826 
827 
828  const static double epsilon = 1e-6 ;
829 
830  std::vector< std::pair<Vector3D, Vector3D> > lines ;
831 
832  //--------------------------------------------
833  // check if there are lines defined in the VolSurface :
834  const std::vector< std::pair<Vector3D, Vector3D> >& local_lines = _volSurf.getLines() ;
835 
836  if( local_lines.size() > 0 ) {
837  unsigned n=local_lines.size() ;
838  lines.reserve( n ) ;
839 
840  for( unsigned i=0;i<n;++i){
841 
842  DDSurfaces::Vector3D av,bv;
843  _wtM->LocalToMaster( local_lines[i].first , av.array() ) ;
844  _wtM->LocalToMaster( local_lines[i].second , bv.array() ) ;
845 
846  lines.push_back( std::make_pair( av, bv ) );
847  }
848 
849  return lines ;
850  }
851  //--------------------------------------------
852 
853 
854  // get local and global surface vectors
855  const DDSurfaces::Vector3D& lu = _volSurf.u() ;
856  // const DDSurfaces::Vector3D& lv = _volSurf.v() ;
857  const DDSurfaces::Vector3D& ln = _volSurf.normal() ;
859 
860  Volume vol = volume() ;
861  const TGeoShape* shape = vol->GetShape() ;
862 
863 
864  if( type().isPlane() ) {
865 
866  if( shape->IsA() == TGeoBBox::Class() ) {
867 
868  TGeoBBox* box = ( TGeoBBox* ) shape ;
869 
870  DDSurfaces::Vector3D boxDim( box->GetDX() , box->GetDY() , box->GetDZ() ) ;
871 
872 
873  bool isYZ = std::fabs( ln.x() - 1.0 ) < epsilon ; // normal parallel to x
874  bool isXZ = std::fabs( ln.y() - 1.0 ) < epsilon ; // normal parallel to y
875  bool isXY = std::fabs( ln.z() - 1.0 ) < epsilon ; // normal parallel to z
876 
877 
878  if( isYZ || isXZ || isXY ) { // plane is parallel to one of the box' sides -> need 4 vertices from box dimensions
879 
880  // if isYZ :
881  unsigned uidx = 1 ;
882  unsigned vidx = 2 ;
883 
884  DDSurfaces::Vector3D ubl( 0., 1., 0. ) ;
885  DDSurfaces::Vector3D vbl( 0., 0., 1. ) ;
886 
887  if( isXZ ) {
888 
889  ubl.fill( 1., 0., 0. ) ;
890  vbl.fill( 0., 0., 1. ) ;
891  uidx = 0 ;
892  vidx = 2 ;
893 
894  } else if( isXY ) {
895 
896  ubl.fill( 1., 0., 0. ) ;
897  vbl.fill( 0., 1., 0. ) ;
898  uidx = 0 ;
899  vidx = 1 ;
900  }
901 
904  _wtM->LocalToMasterVect( ubl , ub.array() ) ;
905  _wtM->LocalToMasterVect( vbl , vb.array() ) ;
906 
907  lines.reserve(4) ;
908 
909  lines.push_back( std::make_pair( _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ) ) ;
910  lines.push_back( std::make_pair( _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ) ) ;
911  lines.push_back( std::make_pair( _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb ) ) ;
912  lines.push_back( std::make_pair( _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb ) ) ;
913 
914  return lines ;
915  }
916 
917  } else if( shape->IsA() == TGeoConeSeg::Class() ) {
918 
919  TGeoCone* cone = ( TGeoCone* ) shape ;
920 
921  // can only deal with special case of z-disk and origin in center of cone
922  if( type().isZDisk() ) { // && lo.rho() < epsilon ) {
923 
924  if( lo.rho() > epsilon ) {
925  // move origin to z-axis
926  lo.x() = 0. ;
927  lo.y() = 0. ;
928  }
929 
930 
931  double zhalf = cone->GetDZ() ;
932  double rmax1 = cone->GetRmax1() ;
933  double rmax2 = cone->GetRmax2() ;
934  double rmin1 = cone->GetRmin1() ;
935  double rmin2 = cone->GetRmin2() ;
936 
937  // two circles around origin
938  // get radii at position of plane
939  double r0 = rmin1 + ( rmin2 - rmin1 ) / ( 2. * zhalf ) * ( zhalf + lo.z() ) ;
940  double r1 = rmax1 + ( rmax2 - rmax1 ) / ( 2. * zhalf ) * ( zhalf + lo.z() ) ;
941 
942 
943  unsigned n = nMax / 4 ;
944  double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
945 
946  for( unsigned i = 0 ; i < n ; ++i ) {
947 
948  Vector3D rv00( r0*sin( i *dPhi ) , r0*cos( i *dPhi ) , 0. ) ;
949  Vector3D rv01( r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi ) , 0. ) ;
950 
951  Vector3D rv10( r1*sin( i *dPhi ) , r1*cos( i *dPhi ) , 0. ) ;
952  Vector3D rv11( r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi ) , 0. ) ;
953 
954 
955  Vector3D pl0 = lo + rv00 ;
956  Vector3D pl1 = lo + rv01 ;
957 
958  Vector3D pl2 = lo + rv10 ;
959  Vector3D pl3 = lo + rv11 ;
960 
961 
962  Vector3D pg0,pg1,pg2,pg3 ;
963 
964  _wtM->LocalToMaster( pl0, pg0.array() ) ;
965  _wtM->LocalToMaster( pl1, pg1.array() ) ;
966  _wtM->LocalToMaster( pl2, pg2.array() ) ;
967  _wtM->LocalToMaster( pl3, pg3.array() ) ;
968 
969  lines.push_back( std::make_pair( pg0, pg1 ) ) ;
970  lines.push_back( std::make_pair( pg2, pg3 ) ) ;
971  }
972 
973  //add some vertical and horizontal lines so that the disc is seen in the rho-z projection
974 
975  n = 4 ; dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
976 
977  for( unsigned i = 0 ; i < n ; ++i ) {
978 
979  Vector3D rv0( r0*sin( i * dPhi ) , r0*cos( i * dPhi ) , 0. ) ;
980  Vector3D rv1( r1*sin( i * dPhi ) , r1*cos( i * dPhi ) , 0. ) ;
981 
982  Vector3D pl0 = lo + rv0 ;
983  Vector3D pl1 = lo + rv1 ;
984 
985  Vector3D pg0,pg1 ;
986 
987  _wtM->LocalToMaster( pl0, pg0.array() ) ;
988  _wtM->LocalToMaster( pl1, pg1.array() ) ;
989 
990  lines.push_back( std::make_pair( pg0, pg1 ) ) ;
991  }
992 
993  }
994 
995  return lines ;
996  }
997 
998  else if(shape->IsA() == TGeoTrap::Class()) {
999  TGeoTrap* trapezoid = ( TGeoTrap* ) shape;
1000 
1001  double dx1 = trapezoid->GetBl1();
1002  double dx2 = trapezoid->GetTl1();
1003  double dz = trapezoid->GetH1();
1004 
1005  //according to the TGeoTrap definition, the lengths are given such that the normal vector of the surface
1006  //points in the e_z direction.
1007  DDSurfaces::Vector3D ubl( 1., 0., 0. ) ;
1008  DDSurfaces::Vector3D vbl( 0., 1., 0. ) ;
1009 
1010  //the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
1013  _wtM->LocalToMasterVect( ubl , ub.array() ) ;
1014  _wtM->LocalToMasterVect( vbl , vb.array() ) ;
1015 
1016  //the trapezoid is drawn as a set of four lines connecting its four corners
1017  lines.reserve(4) ;
1018  //_o is vector to the origin
1019  lines.push_back( std::make_pair( _o + dx1 * ub - dz * vb , _o + dx2 * ub + dz * vb ) ) ;
1020  lines.push_back( std::make_pair( _o + dx2 * ub + dz * vb , _o - dx2 * ub + dz * vb ) ) ;
1021  lines.push_back( std::make_pair( _o - dx2 * ub + dz * vb , _o - dx1 * ub - dz * vb) ) ;
1022  lines.push_back( std::make_pair( _o - dx1 * ub - dz * vb , _o + dx1 * ub - dz * vb) ) ;
1023 
1024  return lines;
1025  }
1026  //added code by Thorben Quast for simplified set of lines for trapezoids with unequal lengths in x
1027  else if(shape->IsA() == TGeoTrd1::Class()){
1028  TGeoTrd1* trapezoid = ( TGeoTrd1* ) shape;
1029  //all lengths are half length
1030  double dx1 = trapezoid->GetDx1();
1031  double dx2 = trapezoid->GetDx2();
1032  //double dy = trapezoid->GetDy();
1033  double dz = trapezoid->GetDz();
1034 
1035  //the normal vector is parallel to e_y for all geometry cases in CLIC
1036  //if that is at some point not the case anymore, then local plane vectors ubl, vbl
1037  //must be initialized like it is done for the boxes (line 674 following)
1038  DDSurfaces::Vector3D ubl( 1., 0., 0. ) ;
1039  DDSurfaces::Vector3D vbl( 0., 0., 1. ) ;
1040 
1041  //the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
1044  _wtM->LocalToMasterVect( ubl , ub.array() ) ;
1045  _wtM->LocalToMasterVect( vbl , vb.array() ) ;
1046 
1047  //the trapezoid is drawn as a set of four lines connecting its four corners
1048  lines.reserve(4) ;
1049  //_o is vector to the origin
1050  lines.push_back( std::make_pair( _o + dx1 * ub - dz * vb , _o + dx2 * ub + dz * vb ) ) ;
1051  lines.push_back( std::make_pair( _o + dx2 * ub + dz * vb , _o - dx2 * ub + dz * vb ) ) ;
1052  lines.push_back( std::make_pair( _o - dx2 * ub + dz * vb , _o - dx1 * ub - dz * vb) ) ;
1053  lines.push_back( std::make_pair( _o - dx1 * ub - dz * vb , _o + dx1 * ub - dz * vb) ) ;
1054 
1055  return lines;
1056  }
1057  //added code by Thorben Quast for simplified set of lines for trapezoids with unequal lengths in x AND y
1058  else if(shape->IsA() == TGeoTrd2::Class()){
1059  TGeoTrd2* trapezoid = ( TGeoTrd2* ) shape;
1060  //all lengths are half length
1061  double dx1 = trapezoid->GetDx1();
1062  double dx2 = trapezoid->GetDx2();
1063  //double dy1 = trapezoid->GetDy1();
1064  //double dy2 = trapezoid->GetDy2();
1065  double dz = trapezoid->GetDz();
1066 
1067  //the normal vector is parallel to e_y for all geometry cases in CLIC
1068  //if that is at some point not the case anymore, then local plane vectors ubl, vbl
1069  //must be initialized like it is done for the boxes (line 674 following)
1070  DDSurfaces::Vector3D ubl( 1., 0., 0. ) ;
1071  DDSurfaces::Vector3D vbl( 0., 0., 1. ) ;
1072 
1073  //the local span vectors are transformed into the main coordinate system (in LocalToMasterVect())
1076  _wtM->LocalToMasterVect( ubl , ub.array() ) ;
1077  _wtM->LocalToMasterVect( vbl , vb.array() ) ;
1078 
1079  //the trapezoid is drawn as a set of four lines connecting its four corners
1080  lines.reserve(4) ;
1081  //_o is vector to the origin
1082  lines.push_back( std::make_pair( _o + dx1 * ub - dz * vb , _o + dx2 * ub + dz * vb ) ) ;
1083  lines.push_back( std::make_pair( _o + dx2 * ub + dz * vb , _o - dx2 * ub + dz * vb ) ) ;
1084  lines.push_back( std::make_pair( _o - dx2 * ub + dz * vb , _o - dx1 * ub - dz * vb) ) ;
1085  lines.push_back( std::make_pair( _o - dx1 * ub - dz * vb , _o + dx1 * ub - dz * vb) ) ;
1086 
1087  return lines;
1088  }
1089  // ===== default for arbitrary planes in arbitrary shapes =================
1090 
1091  // We create nMax vertices by rotating the local u vector around the normal
1092  // and checking the distance to the volume boundary in that direction.
1093  // This is brute force and not very smart, as many points are created on straight
1094  // lines and the edges are still rounded.
1095  // The alterative would be to compute the true intersections a plane and the most
1096  // common shapes - at least for boxes that should be not too hard. To be done...
1097 
1098  lines.reserve( nMax ) ;
1099 
1100  double dAlpha = 2.* ROOT::Math::Pi() / double( nMax ) ;
1101 
1102  TVector3 norm( ln.x() , ln.y() , ln.z() ) ;
1103 
1104 
1105  DDSurfaces::Vector3D first, previous ;
1106 
1107  for(unsigned i=0 ; i< nMax ; ++i ){
1108 
1109  double alpha = double(i) * dAlpha ;
1110 
1111  TVector3 vec( lu.x() , lu.y() , lu.z() ) ;
1112 
1113  TRotation rot ;
1114  rot.Rotate( alpha , norm );
1115 
1116  TVector3 vecR = rot * vec ;
1117 
1118  DDSurfaces::Vector3D luRot ;
1119  luRot.fill( vecR ) ;
1120 
1121  double dist = shape->DistFromInside( const_cast<double*> (lo.const_array()) , const_cast<double*> (luRot.const_array()) , 3, 0.1 ) ;
1122 
1123  // local point at volume boundary
1124  DDSurfaces::Vector3D lp = lo + dist * luRot ;
1125 
1127 
1128  _wtM->LocalToMaster( lp , gp.array() ) ;
1129 
1130  // std::cout << " **** normal:" << ln << " lu:" << lu << " alpha:" << alpha << " luRot:" << luRot << " lp :" << lp << " gp:" << gp << " dist : " << dist
1131  // << " is point " << gp << " inside : " << shape->Contains( gp.const_array() )
1132  // << " dist from outside for lo,lu " << shape->DistFromOutside( lo.const_array() , lu.const_array() , 3 )
1133  // << " dist from inside for lo,ln " << shape->DistFromInside( lo.const_array() , ln.const_array() , 3 )
1134  // << std::endl;
1135  // shape->Dump() ;
1136 
1137 
1138  if( i > 0 )
1139  lines.push_back( std::make_pair( previous, gp ) ) ;
1140  else
1141  first = gp ;
1142 
1143  previous = gp ;
1144  }
1145  lines.push_back( std::make_pair( previous, first ) ) ;
1146 
1147 
1148  } else if( type().isCylinder() ) {
1149 
1150  // if( shape->IsA() == TGeoTube::Class() ) {
1151  if( shape->IsA() == TGeoConeSeg::Class() ) {
1152 
1153  lines.reserve( nMax ) ;
1154 
1155  TGeoTube* tube = ( TGeoTube* ) shape ;
1156 
1157  double zHalf = tube->GetDZ() ;
1158 
1159  Vector3D zv( 0. , 0. , zHalf ) ;
1160 
1161  double r = lo.rho() ;
1162 
1163 
1164 
1165  unsigned n = nMax / 4 ;
1166  double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
1167 
1168  for( unsigned i = 0 ; i < n ; ++i ) {
1169 
1170  Vector3D rv0( r*sin( i *dPhi ) , r*cos( i *dPhi ) , 0. ) ;
1171  Vector3D rv1( r*sin( (i+1)*dPhi ) , r*cos( (i+1)*dPhi ) , 0. ) ;
1172 
1173  // 4 points on local cylinder
1174 
1175  Vector3D pl0 = zv + rv0 ;
1176  Vector3D pl1 = zv + rv1 ;
1177  Vector3D pl2 = -zv + rv1 ;
1178  Vector3D pl3 = -zv + rv0 ;
1179 
1180  Vector3D pg0,pg1,pg2,pg3 ;
1181 
1182  _wtM->LocalToMaster( pl0, pg0.array() ) ;
1183  _wtM->LocalToMaster( pl1, pg1.array() ) ;
1184  _wtM->LocalToMaster( pl2, pg2.array() ) ;
1185  _wtM->LocalToMaster( pl3, pg3.array() ) ;
1186 
1187  lines.push_back( std::make_pair( pg0, pg1 ) ) ;
1188  lines.push_back( std::make_pair( pg1, pg2 ) ) ;
1189  lines.push_back( std::make_pair( pg2, pg3 ) ) ;
1190  lines.push_back( std::make_pair( pg3, pg0 ) ) ;
1191  }
1192  }
1193  }
1194  return lines ;
1195 
1196  }
1197 
1198  //================================================================================================================
1199 
1200 
1201  Vector3D CylinderSurface::u( const Vector3D& point ) const {
1202 
1203  Vector3D lp , u_val ;
1204  _wtM->MasterToLocal( point , lp.array() ) ;
1205  const DDSurfaces::Vector3D& lu = _volSurf.u( lp ) ;
1206  _wtM->LocalToMasterVect( lu , u_val.array() ) ;
1207  return u_val ;
1208  }
1209 
1210  Vector3D CylinderSurface::v(const Vector3D& point ) const {
1211  Vector3D lp , v_val ;
1212  _wtM->MasterToLocal( point , lp.array() ) ;
1213  const DDSurfaces::Vector3D& lv = _volSurf.v( lp ) ;
1214  _wtM->LocalToMasterVect( lv , v_val.array() ) ;
1215  return v_val ;
1216  }
1217 
1219  Vector3D lp , n ;
1220  _wtM->MasterToLocal( point , lp.array() ) ;
1221  const DDSurfaces::Vector3D& ln = _volSurf.normal( lp ) ;
1222  _wtM->LocalToMasterVect( ln , n.array() ) ;
1223  return n ;
1224  }
1225 
1227 
1228  Vector3D lp;
1229  _wtM->MasterToLocal( point , lp.array() ) ;
1230 
1231  return _volSurf.globalToLocal( lp ) ;
1232  }
1233 
1234 
1236 
1237  Vector3D lp = _volSurf.localToGlobal( point ) ;
1238  Vector3D p ;
1239  _wtM->LocalToMaster( lp , p.array() ) ;
1240 
1241  return p ;
1242  }
1243 
1244  double CylinderSurface::radius() const { return _volSurf.origin().rho() ; }
1245 
1247 
1248 
1249  //================================================================================================================
1250 
1251 
1252  double ConeSurface::radius0() const {
1253 
1254  double theta = _volSurf.v().theta() ;
1255  double l = length_along_v() * cos( theta ) ;
1256 
1257  return origin().rho() - 0.5 * l * tan( theta ) ;
1258  }
1259 
1260  double ConeSurface::radius1() const {
1261 
1262  double theta = _volSurf.v().theta() ;
1263  double l = length_along_v() * cos( theta ) ;
1264 
1265  return origin().rho() + 0.5 * l * tan( theta ) ;
1266  }
1267 
1268  double ConeSurface::z0() const {
1269 
1270  double theta = _volSurf.v().theta() ;
1271  double l = length_along_v() * cos( theta ) ;
1272 
1273  return origin().z() - 0.5 * l ;
1274  }
1275 
1276  double ConeSurface::z1() const {
1277 
1278  double theta = _volSurf.v().theta() ;
1279  double l = length_along_v() * cos( theta ) ;
1280 
1281  return origin().z() + 0.5 * l ;
1282  }
1283 
1285 
1286 
1287  //================================================================================================================
1288 
1289  } // namespace
1290 } // namespace
1291 
1292 
virtual double innerThickness() const
Definition: Surface.cpp:611
virtual std::vector< std::pair< Vector3D, Vector3D > > getLines(unsigned nMax=100)
Definition: Surface.cpp:245
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:674
Handle class holding a placed volume (also called physical volume)
Definition: Volumes.h:135
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:241
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:33
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:34
Volume volume() const
Logical volume of this placement.
Definition: Volumes.cpp:389
const Vector3D & fill(const T &v)
fill vector from arbitrary class that defines operator[]
Definition: Vector3D.h:70
bool isZDisk() const
true if this is a plane orthogonal to Z
Definition: ISurface.h:220
virtual double length_along_v() const
Definition: Surface.cpp:240
virtual long64 id() const
The id of this surface - corresponds to DetElement id.
Definition: Surface.cpp:603
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:233
double u() const
Definition: Vector2D.h:21
VolCylinderImpl()
default c'tor
Definition: Surface.h:379
const char * name() const
Access the object name (or "" if not supported by the object)
Definition: Handle.inl:36
virtual bool insideBounds(const Vector3D &point, double epsilon=1.e-4) const
Checks if the given point lies within the surface.
Definition: Surface.cpp:702
static Cylindrical cylindrical()
Definition: Vector3D.h:266
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:608
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:657
#define M_PI
Definition: Handle.h:39
double x() const
Definition: Vector3D.h:91
virtual std::vector< std::pair< Vector3D, Vector3D > > getLines(unsigned nMax=100)
Definition: Surface.cpp:218
TGeoMatrix * _wtM
Definition: Surface.h:517
static Spherical spherical()
Definition: Vector3D.h:267
IFACE * addExtension(CONCRETE *c) const
Extend the detector element with an arbitrary structure accessible by the type.
Definition: Detector.h:297
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:406
virtual void setU(const Vector3D &u)
setter for daughter classes
Definition: Surface.cpp:25
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:391
virtual long64 id() const
The id of this surface.
Definition: Surface.cpp:30
Geometry::DetElement _det
Definition: Surface.h:515
VolSurface _volSurf
Definition: Surface.h:516
Object * data() const
Check if placement is properly instrumented.
Definition: Volumes.cpp:371
virtual const IMaterial & outerMaterial() const
Access to the material in direction of the normal.
Definition: Surface.cpp:637
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:193
double phi() const
Definition: Vector3D.h:130
Out version of the std auto_ptr implementation base either on auto_ptr or unique_ptr.
Definition: Memory.h:43
virtual double length_along_u() const
Definition: Surface.cpp:613
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:284
virtual ~SurfaceList()
d'tor deletes all owned surfaces
Definition: Surface.cpp:499
return e
Definition: Volumes.cpp:297
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:315
std::vector< std::pair< Material, double > > MaterialVec
virtual const Vector3D & origin() const
Definition: Surface.cpp:232
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Access to the normal direction at the given point.
Definition: Surface.cpp:35
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:1210
T * ptr() const
Access to the held object.
Definition: Handle.h:149
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Access to the normal direction at the given point.
Definition: Surface.cpp:609
long long int long64
Definition: ISurface.h:14
virtual const IMaterial & outerMaterial() const
Access to the material in direction of the normal.
Definition: Surface.cpp:236
virtual double z0() const
the start z of the cone
Definition: Surface.cpp:1268
Handle class holding a placed volume (also called physical volume)
Definition: Volumes.h:237
virtual const IMaterial & innerMaterial() const
Access to the material in opposite direction of the normal.
Definition: Surface.cpp:618
double theta() const
Definition: Vector3D.h:167
virtual void setV(const Vector3D &v)
setter for daughter classes
Definition: Surface.cpp:26
virtual double radius() const
the radius of the cylinder (rho of the origin vector)
Definition: Surface.cpp:1244
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:252
virtual long64 id() const
The id of this surface - always 0 for VolSurfaces.
Definition: Surface.cpp:227
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:38
virtual double length_along_v() const
Definition: Surface.cpp:614
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:1201
VolumeID volumeID() const
The cached VolumeID of this subdetector element.
Definition: Detector.cpp:303
virtual double length_along_v() const
Definition: Surface.cpp:130
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:234
static const double um
Definition: DD4hepUnits.h:67
virtual Vector3D volumeOrigin() const
Definition: Surface.cpp:681
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:372
VolConeImpl()
default c'tor
Definition: Surface.h:425
virtual double z1() const
the end z of the cone
Definition: Surface.cpp:1276
bool checkParallelToZ(const ISurface &surf, double epsilon=1.e-6) const
Definition: ISurface.h:236
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:422
const MaterialVec & materialsBetween(const DDSurfaces::Vector3D &p0, const DDSurfaces::Vector3D &p1, double epsilon=1e-4)
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:54
MaterialData createAveragedMaterial(const MaterialVec &materials)
virtual bool insideBounds(const Vector3D &point, double epsilon=1e-4) const
Checks if the given point lies within the surface.
Definition: Surface.cpp:242
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:290
virtual double innerThickness() const
Definition: Surface.cpp:237
virtual bool insideBounds(const Vector3D &point, double epsilon=1e-4) const
Checks if the given point lies within the surface.
Definition: Surface.cpp:196
virtual const SurfaceType & type() const
Definition: Surface.cpp:605
double z() const
Definition: Vector3D.h:97
PlacedVolume placement() const
Access to the physical volume of this detector element.
Definition: Detector.cpp:279
void setOuterMaterial(const IMaterial &mat)
set the outer Materal
Definition: Surface.h:295
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:1218
static const double g
Definition: DD4hepUnits.h:162
Surface()
default c'tor
Definition: Surface.h:528
void setProperty(unsigned prop, bool val=true)
set the given peorperty
Definition: ISurface.h:183
virtual double outerThickness() const
Definition: Surface.cpp:238
virtual double Z() const =0
averaged proton number
double y() const
Definition: Vector3D.h:94
virtual double length_along_u() const
Definition: Surface.cpp:239
virtual Vector3D center() const
the center of the cylinder
Definition: Surface.cpp:1246
Handle< NamedObject > Ref_t
Default Ref_t definition describing named objects.
Definition: Handle.h:176
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:607
Geometry::Volume volume() const
the volume to which this surface is attached.
Definition: Surface.h:227
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:230
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:1235
virtual double outerThickness() const
Definition: Surface.cpp:612
virtual std::vector< std::pair< DDSurfaces::Vector3D, DDSurfaces::Vector3D > > getLines(unsigned nMax=100)
create outer bounding lines for the given symmetry of the polyhedron
Definition: Surface.cpp:447
virtual double innerThickness() const
Definition: Surface.cpp:63
virtual const Vector3D & origin() const
Definition: Surface.cpp:610
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:1226
virtual const SurfaceType & type() const
Definition: Surface.cpp:228
Handle class describing a detector element.
Definition: Detector.h:172
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:229
VolSurfaceList * volSurfaceList(DetElement &det)
Definition: Surface.cpp:511
SurfaceType _type
Definition: Surface.h:521
virtual double length_along_u() const
Definition: Surface.cpp:67
Vector3D unit() const
Definition: Vector3D.h:187
virtual double outerThickness() const
Definition: Surface.cpp:64
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:693
View * v
Definition: MultiView.cpp:30
IFACE * extension() const
Access extension element by the type.
Definition: Detector.h:302
virtual const Vector3D & origin() const
Definition: Surface.cpp:36
bool checkOrthogonalToZ(const ISurface &surf, double epsilon=1.e-6) const
Definition: ISurface.h:253
virtual std::vector< std::pair< Vector3D, Vector3D > > getLines(unsigned nMax=100)
Definition: Surface.cpp:825
Vector3D cross(const Vector3D &v) const
Definition: Vector3D.h:179
static const double second
Definition: DD4hepUnits.h:112
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:378
virtual const IMaterial & outerMaterial() const
Access to the material in direction of the normal.
Definition: Surface.cpp:62
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:385
virtual void setNormal(const Vector3D &n)
setter for daughter classes
Definition: Surface.cpp:27
const double * const_array() const
direct access to data as const double*
Definition: Vector3D.h:199
double * array()
direct access to data as double* - allows modification
Definition: Vector3D.h:204
virtual const IMaterial & innerMaterial() const
Access to the material in opposite direction of the normal.
Definition: Surface.cpp:235
virtual Vector3D center() const
the center of the cone
Definition: Surface.cpp:1284
bool findVolume(PlacedVolume pv, Volume theVol, std::list< PlacedVolume > &volList)
Definition: Surface.cpp:531
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:277
double v() const
Definition: Vector2D.h:23
virtual const IMaterial & innerMaterial() const
Access to the material in opposite direction of the normal.
Definition: Surface.cpp:61
virtual const SurfaceType & type() const
Definition: Surface.cpp:32
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Access to the normal direction at the given point.
Definition: Surface.cpp:231
virtual void setOrigin(const Vector3D &o)
setter for daughter classes
Definition: Surface.cpp:28
void setInnerMaterial(const IMaterial &mat)
set the innerMaterial
Definition: Surface.h:292
virtual double radius0() const
the start radius of the cone
Definition: Surface.cpp:1252
double rho() const
Definition: Vector3D.h:136
Geometry::Volume volume() const
The volume that has the surface attached.
Definition: Surface.h:548
virtual double radius1() const
the end radius of the cone
Definition: Surface.cpp:1260
const TGeoHMatrix & worldTransformation() const
Set detector element for reference transformations. Will delete existing reference trafo...
Definition: Detector.cpp:375
TGeoShape TGeoMedium * m
Definition: Volumes.cpp:294
bool isCone() const
true if this a conical surface
Definition: ISurface.h:198
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:302