DD4hep  1.28.0
Detector Description Toolkit for High Energy Physics
Surface.cpp
Go to the documentation of this file.
1 //==========================================================================
2 // AIDA Detector description implementation
3 //--------------------------------------------------------------------------
4 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
5 // All rights reserved.
6 //
7 // For the licensing terms see $DD4hepINSTALL/LICENSE.
8 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
9 //
10 // Author : F.Gaede
11 //
12 //==========================================================================
13 #include "DDRec/Surface.h"
15 #include "DD4hep/Memory.h"
16 
17 #include "DDRec/MaterialManager.h"
18 
19 #include <cmath>
20 #include <memory>
21 #include <exception>
22 
23 #include "TGeoMatrix.h"
24 #include "TGeoShape.h"
25 #include "TRotation.h"
26 //TGeoTrd1 is apparently not included by defautl
27 #include "TGeoTrd1.h"
28 
29 namespace dd4hep {
30  namespace rec {
31 
32  using namespace detail ;
33 
34 
35  //======================================================================================================
36 
37  void VolSurfaceBase::setU(const Vector3D& u_val) { _u = u_val ; }
38  void VolSurfaceBase::setV(const Vector3D& v_val) { _v = v_val ; }
39  void VolSurfaceBase::setNormal(const Vector3D& n) { _n = n ; }
40  void VolSurfaceBase::setOrigin(const Vector3D& o) { _o = o ; }
41 
42  long64 VolSurfaceBase::id() const { return _id ; }
43 
44  const SurfaceType& VolSurfaceBase::type() const { return _type ; }
45  Vector3D VolSurfaceBase::u(const Vector3D& /*point*/) const { return _u ; }
46  Vector3D VolSurfaceBase::v(const Vector3D& /*point*/) const { return _v ; }
47  Vector3D VolSurfaceBase::normal(const Vector3D& /*point*/) const { return _n ; }
48  const Vector3D& VolSurfaceBase::origin() const { return _o ;}
49 
51 
52  Vector3D p = point - origin() ;
53 
54  // create new orthogonal unit vectors
55  // FIXME: these vectors should be cached really ...
56 
57  double uv = u() * v() ;
58  Vector3D uprime = ( u() - uv * v() ).unit() ;
59  Vector3D vprime = ( v() - uv * u() ).unit() ;
60  double uup = u() * uprime ;
61  double vvp = v() * vprime ;
62 
63  return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
64  }
65 
67 
68  Vector3D g = origin() + point[0] * u() + point[1] * v() ;
69 
70  return g ;
71  }
72 
73  const IMaterial& VolSurfaceBase::innerMaterial() const{ return _innerMat ; }
74  const IMaterial& VolSurfaceBase::outerMaterial() const { return _outerMat ; }
75  double VolSurfaceBase::innerThickness() const { return _th_i ; }
76  double VolSurfaceBase::outerThickness() const { return _th_o ; }
77 
78 
80 
81  const Vector3D& o = this->origin() ;
82  const Vector3D& u_val = this->u( o ) ;
83  Vector3D um = -1. * u_val ;
84 
85  double dist_p = 0. ;
86  double dist_m = 0. ;
87 
88 
89  // std::cout << " VolSurfaceBase::length_along_u() : o = " << o << " u = " << this->u( o )
90  // << " -u = " << um << std::endl ;
91 
92 
93  if( volume()->GetShape()->Contains( o.const_array() ) ){
94 
95  dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
96  const_cast<double*> ( u_val.const_array() ) ) ;
97  dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
98  const_cast<double*> ( um.array() ) ) ;
99 
100 
101  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
102  // << " dist_p " << dist_p
103  // << " dist_m " << dist_m
104  // << std::endl ;
105 
106 
107  } else{
108 
109  dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
110  const_cast<double*> ( u_val.const_array() ) ) ;
111  dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
112  const_cast<double*> ( um.array() ) ) ;
113 
114  dist_p *= 1.0001 ;
115  dist_m *= 1.0001 ;
116 
117  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
118  // << " dist_p " << dist_p
119  // << " dist_m " << dist_m
120  // << std::endl ;
121 
122  Vector3D o_1 = this->origin() + dist_p * u_val ;
123  Vector3D o_2 = this->origin() + dist_m * um ;
124 
125  dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) ,
126  const_cast<double*> ( u_val.const_array() ) ) ;
127 
128  dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) ,
129  const_cast<double*> ( um.array() ) ) ;
130 
131  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
132  // << " dist_p " << dist_p
133  // << " dist_m " << dist_m
134  // << std::endl ;
135  }
136 
137  return dist_p + dist_m ;
138 
139 
140  }
141 
143 
144  const Vector3D& o = this->origin() ;
145  const Vector3D& v_val = this->v( o ) ;
146  Vector3D vm = -1. * v_val ;
147 
148  double dist_p = 0. ;
149  double dist_m = 0. ;
150 
151 
152  // std::cout << " VolSurfaceBase::length_along_u() : o = " << o << " u = " << this->u( o )
153  // << " -u = " << vm << std::endl ;
154 
155 
156  if( volume()->GetShape()->Contains( o.const_array() ) ){
157 
158  dist_p = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
159  const_cast<double*> ( v_val.const_array() ) ) ;
160  dist_m = volume()->GetShape()->DistFromInside( const_cast<double*> ( o.const_array() ) ,
161  const_cast<double*> ( vm.array() ) ) ;
162 
163 
164  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
165  // << " dist_p " << dist_p
166  // << " dist_m " << dist_m
167  // << std::endl ;
168 
169 
170  } else{
171 
172  dist_p = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
173  const_cast<double*> ( v_val.const_array() ) ) ;
174  dist_m = volume()->GetShape()->DistFromOutside( const_cast<double*> ( o.const_array() ) ,
175  const_cast<double*> ( vm.array() ) ) ;
176 
177  dist_p *= 1.0001 ;
178  dist_m *= 1.0001 ;
179 
180  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
181  // << " dist_p " << dist_p
182  // << " dist_m " << dist_m
183  // << std::endl ;
184 
185  Vector3D o_1 = this->origin() + dist_p * v_val ;
186  Vector3D o_2 = this->origin() + dist_m * vm ;
187 
188  dist_p += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_1.const_array() ) ,
189  const_cast<double*> ( v_val.const_array() ) ) ;
190 
191  dist_m += volume()->GetShape()->DistFromInside( const_cast<double*> ( o_2.const_array() ) ,
192  const_cast<double*> ( vm.array() ) ) ;
193 
194  // std::cout << " VolSurfaceBase::length_along_u() : shape contains(o) = " << volume()->GetShape()->Contains( o.const_array() )
195  // << " dist_p " << dist_p
196  // << " dist_m " << dist_m
197  // << std::endl ;
198  }
199 
200  return dist_p + dist_m ;
201 
202  }
203 
204 
205  double VolSurfaceBase::distance(const Vector3D& /*point*/ ) const { return 1.e99 ; }
206 
208  bool VolSurfaceBase::insideBounds(const Vector3D& point, double epsilon) const {
209 
210 #if 0
211 
212  bool inShape = ( type().isUnbounded() ? true : volume()->GetShape()->Contains( point.const_array() ) ) ;
213 
214  double dist = std::abs ( distance( point ) ) ;
215 
216  std::cout << " ** Surface::insideBound( " << point << " ) - distance = " << dist
217  << " origin = " << origin() << " normal = " << normal()
218  << " p * n = " << point * normal()
219  << " isInShape : " << inShape << std::endl ;
220 
221  return dist < epsilon && inShape ;
222 #else
223 
224  if( type().isUnbounded() ){
225 
226  return std::abs ( distance( point ) ) < epsilon ;
227 
228  } else {
229 
230  return ( std::abs ( distance( point ) ) < epsilon && volume()->GetShape()->Contains( point.const_array() ) ) ;
231  }
232 
233 #endif
234 
235  }
236 
237 
238  std::vector< std::pair<Vector3D, Vector3D> > VolSurfaceBase::getLines(unsigned ) {
239  // dummy implementation returning empty set
240  std::vector< std::pair<Vector3D, Vector3D> > lines ;
241  return lines ;
242  }
243 
244  //===================================================================
245  // simple wrapper methods forwarding the call to the implementation object
246 
247  long64 VolSurface::id() const { return _surf->id() ; }
248  const SurfaceType& VolSurface::type() const { return _surf->type() ; }
249  Vector3D VolSurface::u( const Vector3D& point ) const { return _surf->u(point) ; }
250  Vector3D VolSurface::v(const Vector3D& point ) const { return _surf->v(point) ; }
251  Vector3D VolSurface::normal(const Vector3D& point ) const { return _surf->normal(point) ; }
252  const Vector3D& VolSurface::origin() const { return _surf->origin() ;}
253  Vector2D VolSurface::globalToLocal( const Vector3D& point) const { return _surf->globalToLocal( point ) ; }
254  Vector3D VolSurface::localToGlobal( const Vector2D& point) const { return _surf->localToGlobal( point) ; }
255  const IMaterial& VolSurface::innerMaterial() const{ return _surf->innerMaterial() ; }
256  const IMaterial& VolSurface::outerMaterial() const { return _surf->outerMaterial() ; }
257  double VolSurface::innerThickness() const { return _surf->innerThickness() ; }
258  double VolSurface::outerThickness() const { return _surf->outerThickness() ; }
259  double VolSurface::length_along_u() const { return _surf->length_along_u() ; }
260  double VolSurface::length_along_v() const { return _surf->length_along_v() ; }
261  double VolSurface::distance(const Vector3D& point ) const { return _surf->distance( point ) ; }
262  bool VolSurface::insideBounds(const Vector3D& point, double epsilon) const {
263  return _surf->insideBounds( point, epsilon ) ;
264  }
265  std::vector< std::pair<Vector3D, Vector3D> > VolSurface::getLines(unsigned nMax) {
266  return _surf->getLines(nMax) ;
267  }
268 
269  //===================================================================
270 
272  double VolPlaneImpl::distance(const Vector3D& point ) const {
273  return ( point - origin() ) * normal() ;
274  }
275  //======================================================================================================
276 
278  double thickness_inner ,double thickness_outer, Vector3D o ) :
279 
280  VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , Vector3D() , Vector3D() , o , vol, 0) {
281  Vector3D v_val( 0., 0., 1. ) ;
282  Vector3D o_rphi( o.x() , o.y() , 0. ) ;
283  Vector3D n = o_rphi.unit() ;
284  Vector3D u_val = v_val.cross( n ) ;
285 
286  setU( u_val ) ;
287  setV( v_val ) ;
288  setNormal( n ) ;
289 
292  _type.setProperty( SurfaceType::Cone , false ) ;
293  _type.checkParallelToZ( *this ) ;
294  _type.checkOrthogonalToZ( *this ) ;
295  }
296 
297  Vector3D VolCylinderImpl::u(const Vector3D& point ) const {
298 
299  Vector3D n( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
300 
301  return v().cross( n ) ;
302  }
303 
304  Vector3D VolCylinderImpl::normal(const Vector3D& point ) const {
305 
306  // normal is just given by phi of the point
307  return Vector3D( 1. , point.phi() , 0. , Vector3D::cylindrical ) ;
308  }
309 
311 
312  // cylinder is parallel to v here so u is Z and v is r *Phi
313  double phi = point.phi() - origin().phi() ;
314 
315  while( phi < -M_PI ) phi += 2.*M_PI ;
316  while( phi > M_PI ) phi -= 2.*M_PI ;
317 
318  return Vector2D( origin().rho() * phi, point.z() - origin().z() ) ;
319  }
320 
321 
323 
324  double z = point.v() + origin().z() ;
325  double phi = point.u() / origin().rho() + origin().phi() ;
326 
327  while( phi < -M_PI ) phi += 2.*M_PI ;
328  while( phi > M_PI ) phi -= 2.*M_PI ;
329 
330  return Vector3D( origin().rho() , phi, z , Vector3D::cylindrical ) ;
331  }
332 
333 
335  double VolCylinderImpl::distance(const Vector3D& point ) const {
336 
337  return point.rho() - origin().rho() ;
338  }
339 
340  //================================================================================================================
342  double thickness_inner ,double thickness_outer, Vector3D v_val, Vector3D o_val ) :
343 
344  VolSurfaceBase(typ, thickness_inner, thickness_outer, Vector3D() , v_val , Vector3D() , Vector3D() , vol, 0) {
345 
346  Vector3D o_rphi( o_val.x() , o_val.y() , 0. ) ;
347 
348  // sanity check: v and o have to have a common phi
349  double dphi = v_val.phi() - o_rphi.phi() ;
350  while( dphi < -M_PI ) dphi += 2.*M_PI ;
351  while( dphi > M_PI ) dphi -= 2.*M_PI ;
352 
353  if( std::fabs( dphi ) > 1e-6 ){
354  std::stringstream sst ; sst << "VolConeImpl::VolConeImpl() - incompatibel vector v and o given "
355  << v_val << " - " << o_val ;
356  throw std::runtime_error( sst.str() ) ;
357  }
358 
359  double theta = v_val.theta() ;
360 
361  Vector3D n( 1. , v_val.phi() , ( theta + M_PI/2. ) , Vector3D::spherical ) ;
362  Vector3D u_val = v_val.cross( n ) ;
363 
364  setU( u_val ) ;
365  setOrigin( o_rphi ) ;
366  setNormal( n ) ;
367 
368  // set helper variable for faster computations (describe cone with tip at origin)
369  _tanTheta = std::tan( theta ) ;
370  double tipoffset = o_val.rho() / _tanTheta ; // distance from tip to origin.z()
371  _ztip = o_val.z() - tipoffset ;
372 
373  double dist_p = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) ,
374  const_cast<double*> ( v_val.const_array() ) ) ;
375  Vector3D vm = -1. * v_val ;
376  double dist_m = vol->GetShape()->DistFromInside( const_cast<double*> ( o_val.const_array() ) ,
377  const_cast<double*> ( vm.array() ) ) ;
378 
379  double costh = std::cos( theta) ;
380  _zt0 = tipoffset - dist_m * costh ;
381  _zt1 = tipoffset + dist_p * costh ;
382 
383 
389  }
390 
391 
392  Vector3D VolConeImpl::v(const Vector3D& point ) const {
393  // just take phi from point
394  Vector3D av( 1. , point.phi() , _v.theta() , Vector3D::spherical ) ;
395  return av ;
396  }
397 
398  Vector3D VolConeImpl::u(const Vector3D& point ) const {
399  // compute from v X n
400  const Vector3D& av = this->v( point ) ;
401  const Vector3D& n = normal( point ) ;
402  return av.cross( n ) ;
403  }
404 
405  Vector3D VolConeImpl::normal(const Vector3D& point ) const {
406  // just take phi from point
407  Vector3D n( 1. , point.phi() , _n.theta() , Vector3D::spherical ) ;
408  return n ;
409  }
410 
412 
413  // cone is parallel to z here, so u is r *Phi and v is "along" z
414  double phi = point.phi() - origin().phi() ;
415 
416  while( phi < -M_PI ) phi += 2.*M_PI ;
417  while( phi > M_PI ) phi -= 2.*M_PI ;
418 
419 
420  double r = ( point.z() - _ztip ) * _tanTheta ;
421 
422  return Vector2D( r*phi, ( point.z() - origin().z() ) / cos( _v.theta() ) ) ;
423  }
424 
425 
427 
428  double z = point.v() * cos( _v.theta() ) + origin().z() ;
429 
430  double r = ( z - _ztip ) * _tanTheta ;
431 
432  double phi = point.u() / r + origin().phi() ;
433 
434  while( phi < -M_PI ) phi += 2.*M_PI ;
435  while( phi > M_PI ) phi -= 2.*M_PI ;
436 
437  return Vector3D( r , phi, z , Vector3D::cylindrical ) ;
438  }
439 
440 
442  double VolConeImpl::distance(const Vector3D& point ) const {
443 
444  // // if the point is in the other hemispere we return the distance to origin
445  // // -> this assumes that the cones do not cross the xy-plane ...
446  // // otherwise we get the distance to the mirrored part of the cone
447  // // needs more thought ..
448  // if( origin().z() * point.z() < 0. )
449  // return point.r() ;
450 
451  //fixme: there are probably faster ways to compute this
452  // e.g by using the intercept theorem - tbd. ...
453  // const Vector2D& lp = globalToLocal( point ) ;
454  // const Vector3D& gp = localToGlobal( lp ) ;
455 
456  // Vector3D dz = point - gp ;
457 
458  //return dz * normal( point ) ;
459 
460  double zp = point.z() - _ztip ;
461  double r = point.rho() - zp * _tanTheta ;
462  return r * std::cos( _v.theta() ) ;
463 
464  }
465 
467  std::vector< std::pair<Vector3D, Vector3D> > VolConeImpl::getLines(unsigned nMax){
468 
469  std::vector< std::pair<Vector3D, Vector3D> > lines ;
470 
471  lines.reserve( nMax ) ;
472 
473  double theta = v().theta() ;
474  double half_length = 0.5 * length_along_v() * cos( theta ) ;
475 
476  Vector3D zv( 0. , 0. , half_length ) ;
477 
478  double dr = half_length * tan( theta ) ;
479 
480  double r0 = origin().rho() - dr ;
481  double r1 = origin().rho() + dr ;
482 
483 
484  unsigned n = nMax / 4 ;
485  double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
486 
487  for( unsigned i = 0 ; i < n ; ++i ) {
488 
489  Vector3D r0v0( r0*sin( i *dPhi ) , r0*cos( i *dPhi ) , 0. ) ;
490  Vector3D r0v1( r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi ) , 0. ) ;
491 
492  Vector3D r1v0( r1*sin( i *dPhi ) , r1*cos( i *dPhi ) , 0. ) ;
493  Vector3D r1v1( r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi ) , 0. ) ;
494 
495  Vector3D pl0 = zv + r1v0 ;
496  Vector3D pl1 = zv + r1v1 ;
497  Vector3D pl2 = -zv + r0v1 ;
498  Vector3D pl3 = -zv + r0v0 ;
499 
500  lines.emplace_back( pl0, pl1 );
501  lines.emplace_back( pl1, pl2 );
502  lines.emplace_back( pl2, pl3 );
503  lines.emplace_back( pl3, pl0 );
504  }
505  return lines;
506  }
507 
508  //================================================================================================================
509 
510 
512  // // delete all surfaces attached to this volume
513  // for( VolSurfaceList::iterator i=begin(), n=end() ; i !=n ; ++i ) {
514  // i->clear() ;
515  // }
516  }
517  //=======================================================
518 
520  if( _isOwner ) {
521  // delete all surfaces attached to this volume
522  std::for_each(begin(), end(), detail::deleteObject<ISurface>);
523  }
524  }
525 
526  //================================================================================================================
527 
529  VolSurfaceList* list = det.extension< VolSurfaceList >(false);
530  if ( !list ) {
531  list = det.addExtension<VolSurfaceList >(new VolSurfaceList);
532  }
533  return list ;
534  }
535 
536 
537  //======================================================================================================================
538 
539  bool findVolume( PlacedVolume pv, Volume theVol, std::list< PlacedVolume >& volList ) {
540 
541 
542  volList.emplace_back( pv ) ;
543 
544  // unsigned count = volList.size() ;
545  // for(unsigned i=0 ; i < count ; ++i) {
546  // std::cout << " **" ;
547  // }
548  // std::cout << " searching for volume: " << theVol.name() << " " << std::hex << theVol.ptr() << " <-> pv.volume : " << pv.name() << " " << pv.volume().ptr()
549  // << " pv.volume().ptr() == theVol.ptr() " << (pv.volume().ptr() == theVol.ptr() )
550  // << std::endl ;
551 
552 
553  if( pv.volume().ptr() == theVol.ptr() ) {
554 
555  return true ;
556 
557  } else {
558 
559  //--------------------------------
560 
561  const TGeoNode* node = pv.ptr();
562 
563  if ( !node ) {
564 
565  // std::cout << " *** findVolume: Invalid placement: - node pointer Null for volume: " << pv.name() << std::endl ;
566 
567  throw std::runtime_error("*** findVolume: Invalid placement: - node pointer Null ! " + std::string( pv.name() ) );
568  }
569  // Volume vol = pv.volume();
570 
571  // std::cout << " ndau = " << node->GetNdaughters() << std::endl ;
572 
573  for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
574 
575  TGeoNode* daughter = node->GetDaughter(idau);
576  PlacedVolume placement( daughter );
577 
578  if ( !placement.data() ) {
579  throw std::runtime_error("*** findVolume: Invalid not instrumented placement:"+std::string(daughter->GetName())
580  +" [Internal error -- bad detector constructor]");
581  }
582 
583  PlacedVolume pv_dau( daughter );
584 
585  if( findVolume( pv_dau , theVol , volList ) ) {
586 
587  // std::cout << " ----- found in daughter volume !!! " << std::hex << pv_dau.volume().ptr() << std::endl ;
588 
589  return true ;
590  }
591  }
592 
593  // ------- not found:
594 
595  volList.pop_back() ;
596 
597  return false ;
598  //--------------------------------
599 
600  }
601  }
602 
603  //======================================================================================================================
604 
605  Surface::Surface( DetElement det, VolSurface volSurf ) : _det( det) , _volSurf( volSurf ),
606  _wtM() , _id( 0) , _type( _volSurf.type() ) {
607 
608  initialize() ;
609  }
610 
611  long64 Surface::id() const { return _id ; }
612 
613  const SurfaceType& Surface::type() const { return _type ; }
614 
615  Vector3D Surface::u(const Vector3D& /*point*/) const { return _u ; }
616  Vector3D Surface::v(const Vector3D& /*point*/) const { return _v ; }
617  Vector3D Surface::normal(const Vector3D& /*point*/) const { return _n ; }
618  const Vector3D& Surface::origin() const { return _o ;}
619  double Surface::innerThickness() const { return _volSurf.innerThickness() ; }
620  double Surface::outerThickness() const { return _volSurf.outerThickness() ; }
621  double Surface::length_along_u() const { return _volSurf.length_along_u() ; }
622  double Surface::length_along_v() const { return _volSurf.length_along_v() ; }
623 
627 
628  const IMaterial& mat = _volSurf.innerMaterial() ;
629 
630  if( ! ( mat.Z() > 0 ) ) {
631 
632  MaterialManager matMgr( _det.placement().volume() ) ;
633 
634  Vector3D p = _o - innerThickness() * _n ;
635 
636  const MaterialVec& materials = matMgr.materialsBetween( _o , p ) ;
637 
638  _volSurf.setInnerMaterial( materials.size() > 1 ?
639  matMgr.createAveragedMaterial( materials ) :
640  materials[0].first ) ;
641  }
642  return mat ;
643  }
644 
646 
647  const IMaterial& mat = _volSurf.outerMaterial() ;
648 
649  if( ! ( mat.Z() > 0 ) ) {
650 
651  MaterialManager matMgr( _det.placement().volume() ) ;
652 
653  Vector3D p = _o + outerThickness() * _n ;
654 
655  const MaterialVec& materials = matMgr.materialsBetween( _o , p ) ;
656 
657  _volSurf.setOuterMaterial( materials.size() > 1 ?
658  matMgr.createAveragedMaterial( materials ) :
659  materials[0].first ) ;
660  }
661  return mat ;
662  }
663 
664 
665  Vector2D Surface::globalToLocal( const Vector3D& point) const {
666 
667  Vector3D p = point - origin() ;
668 
669  // create new orthogonal unit vectors
670  // FIXME: these vectors should be cached really ...
671 
672  double uv = u() * v() ;
673  Vector3D uprime = ( u() - uv * v() ).unit() ;
674  Vector3D vprime = ( v() - uv * u() ).unit() ;
675  double uup = u() * uprime ;
676  double vvp = v() * vprime ;
677 
678  return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
679  }
680 
681 
682  Vector3D Surface::localToGlobal( const Vector2D& point) const {
683 
684  Vector3D g = origin() + point[0] * u() + point[1] * v() ;
685  return g ;
686  }
687 
688 
690 
691  double o_array[3] ;
692 
693  _wtM->LocalToMaster ( Vector3D() , o_array ) ;
694 
695  Vector3D o(o_array) ;
696 
697  return o ;
698  }
699 
700 
701  double Surface::distance(const Vector3D& point ) const {
702 
703  double pa[3] ;
704  _wtM->MasterToLocal( point , pa ) ;
705  Vector3D localPoint( pa ) ;
706 
707  return _volSurf.distance( localPoint ) ;
708  }
709 
710  bool Surface::insideBounds(const Vector3D& point, double epsilon) const {
711 
712  double pa[3] ;
713  _wtM->MasterToLocal( point , pa ) ;
714  Vector3D localPoint( pa ) ;
715 
716  return _volSurf.insideBounds( localPoint , epsilon) ;
717  }
718 
720 
721  // first we need to find the right volume for the local surface in the DetElement's volumes
722  std::list< PlacedVolume > pVList ;
723  PlacedVolume pv = _det.placement() ;
724  Volume theVol = _volSurf.volume() ;
725 
726  if( ! findVolume( pv, theVol , pVList ) ){
727  theVol = _volSurf.volume() ;
728  std::stringstream sst ; sst << " ***** ERROR: Volume " << theVol.name() << " not found for DetElement " << _det.name() << " with surface " ;
729  throw std::runtime_error( sst.str() ) ;
730  }
731 
732  //=========== compute and cache world transform for surface ==========
733  Alignment nominal = _det.nominal();
734  const TGeoHMatrix& wm = nominal.worldTransformation() ;
735 
736 #if 0 // debug
737  wm.Print() ;
738  for( std::list<PlacedVolume>::iterator it= pVList.begin(), n = pVList.end() ; it != n ; ++it ){
739  PlacedVolume pv = *it ;
740  TGeoMatrix* m = pv->GetMatrix();
741  std::cout << " +++ matrix for placed volume : " << std::endl ;
742  m->Print() ;
743  }
744 #endif
745 
746  // need to get the inverse transformation ( see Detector.cpp )
747  // std::auto_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm.Inverse() ) ) ;
748  // has been fixed now, no need to get the inverse anymore:
749  std::unique_ptr<TGeoHMatrix> wtI( new TGeoHMatrix( wm ) ) ;
750 
751  //---- if the volSurface is not in the DetElement's volume, we need to mutliply the path to the volume to the
752  // DetElements world transform
753  for( std::list<PlacedVolume>::iterator it = ++( pVList.begin() ) , n = pVList.end() ; it != n ; ++it ){
754 
755  PlacedVolume pvol = *it ;
756  TGeoMatrix* m = pvol->GetMatrix();
757  // std::cout << " +++ matrix for placed volume : " << std::endl ;
758  // m->Print() ;
759  //wtI->MultiplyLeft( m );
760 
761  wtI->Multiply( m );
762  }
763 
764  // std::cout << " +++ new world transform matrix : " << std::endl ;
765 
766 #if 0 //fixme: which convention to use here - the correct should be wtI, however it is the inverse of what is stored in DetElement ???
767  dd4hep_ptr<TGeoHMatrix> wt( new TGeoHMatrix( wtI->Inverse() ) ) ;
768  wt->Print() ;
769  // cache the world transform for the surface
770  _wtM = std::move(wtI);
771 #else
772  // wtI->Print() ;
773  // cache the world transform for the surface
774  _wtM = std::move(wtI);
775 #endif
776 
777 
778  // ============ now fill the global surface vectors ==========================
779 
780  double ua[3], va[3], na[3], oa[3] ;
781 
782  _wtM->LocalToMasterVect( _volSurf.u() , ua ) ;
783  _wtM->LocalToMasterVect( _volSurf.v() , va ) ;
784  _wtM->LocalToMasterVect( _volSurf.normal() , na ) ;
785  _wtM->LocalToMaster ( _volSurf.origin() , oa ) ;
786 
787  _u.fill( ua ) ;
788  _v.fill( va ) ;
789  _n.fill( na ) ;
790  _o.fill( oa ) ;
791 
792  // std::cout << " --- local and global surface vectors : ------- " << std::endl
793  // << " u : " << _volSurf.u() << " - " << _u << std::endl
794  // << " v : " << _volSurf.v() << " - " << _v << std::endl
795  // << " n : " << _volSurf.normal() << " - " << _n << std::endl
796  // << " o : " << _volSurf.origin() << " - " << _o << std::endl ;
797 
798 
799  // =========== check parallel and orthogonal to Z ===================
800 
801  if( ! _type.isCone() ) {
802  //fixme: workaround for conical surfaces that should always be parallel to z
803  // however the check with the normal does not work here ...
804 
805  _type.checkParallelToZ( *this ) ;
806 
807  _type.checkOrthogonalToZ( *this ) ;
808  }
809 
810  //======== set the unique surface ID from the DetElement ( and placements below ? )
811 
812  // just use the DetElement ID for now ...
813  // or the id set by the user to the VolSurface ...
814  _id = ( _volSurf.id()==0 ? _det.volumeID() : _volSurf.id() ) ;
815 
816  // typedef PlacedVolume::VolIDs IDV ;
817  // DetElement d = _det ;
818  // while( d.isValid() && d.parent().isValid() ){
819  // PlacedVolume pv = d.placement() ;
820  // if( pv.isValid() ){
821  // const IDV& idV = pv.volIDs() ;
822  // std::cout << " VolIDs : " << d.name() << std::endl ;
823  // for( unsigned i=0, n=idV.size() ; i<n ; ++i){
824  // std::cout << " " << idV[i].first << " - " << idV[i].second << std::endl ;
825  // }
826  // }
827  // d = d.parent() ;
828  // }
829 
830  }
831  //===================================================================================================================
832 
833  std::vector< std::pair<Vector3D, Vector3D> > Surface::getLines(unsigned nMax) {
834 
835 
836  const static double epsilon = 1e-6 ;
837 
838  std::vector< std::pair<Vector3D, Vector3D> > lines ;
839 
840  //--------------------------------------------
841  // check if there are lines defined in the VolSurface :
842  const std::vector< std::pair<Vector3D, Vector3D> >& local_lines = _volSurf.getLines() ;
843 
844  if( local_lines.size() > 0 ) {
845  unsigned n=local_lines.size() ;
846  lines.reserve( n ) ;
847 
848  for( unsigned i=0;i<n;++i){
849 
850  Vector3D av,bv;
851  _wtM->LocalToMaster( local_lines[i].first , av.array() ) ;
852  _wtM->LocalToMaster( local_lines[i].second , bv.array() ) ;
853 
854  lines.emplace_back( av, bv );
855  }
856 
857  return lines ;
858  }
859  //--------------------------------------------
860 
861 
862  // get local and global surface vectors
863  const Vector3D& lu = _volSurf.u() ;
864  // const Vector3D& lv = _volSurf.v() ;
865  const Vector3D& ln = _volSurf.normal() ;
866  Vector3D lo = _volSurf.origin() ;
867 
868  Volume vol = volume() ;
869  const TGeoShape* shape = vol->GetShape() ;
870 
871 
872  if( type().isPlane() ) {
873 
874  if( shape->IsA() == TGeoBBox::Class() ) {
875 
876  TGeoBBox* box = ( TGeoBBox* ) shape ;
877 
878  Vector3D boxDim( box->GetDX() , box->GetDY() , box->GetDZ() ) ;
879 
880 
881  bool isYZ = std::fabs( ln.x() - 1.0 ) < epsilon ; // normal parallel to x
882  bool isXZ = std::fabs( ln.y() - 1.0 ) < epsilon ; // normal parallel to y
883  bool isXY = std::fabs( ln.z() - 1.0 ) < epsilon ; // normal parallel to z
884 
885 
886  if( isYZ || isXZ || isXY ) { // plane is parallel to one of the box' sides -> need 4 vertices from box dimensions
887 
888  // if isYZ :
889  unsigned uidx = 1 ;
890  unsigned vidx = 2 ;
891 
892  Vector3D ubl( 0., 1., 0. ) ;
893  Vector3D vbl( 0., 0., 1. ) ;
894 
895  if( isXZ ) {
896 
897  ubl.fill( 1., 0., 0. ) ;
898  vbl.fill( 0., 0., 1. ) ;
899  uidx = 0 ;
900  vidx = 2 ;
901 
902  } else if( isXY ) {
903 
904  ubl.fill( 1., 0., 0. ) ;
905  vbl.fill( 0., 1., 0. ) ;
906  uidx = 0 ;
907  vidx = 1 ;
908  }
909 
910  Vector3D ub ;
911  Vector3D vb ;
912  _wtM->LocalToMasterVect( ubl , ub.array() ) ;
913  _wtM->LocalToMasterVect( vbl , vb.array() ) ;
914 
915  lines.reserve(4) ;
916 
917  lines.emplace_back(_o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb );
918  lines.emplace_back(_o - boxDim[ uidx ] * ub + boxDim[ vidx ] * vb , _o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb );
919  lines.emplace_back(_o - boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb );
920  lines.emplace_back(_o + boxDim[ uidx ] * ub - boxDim[ vidx ] * vb , _o + boxDim[ uidx ] * ub + boxDim[ vidx ] * vb );
921 
922  return lines ;
923  }
924 
925  } else if( shape->InheritsFrom("TGeoTube") || shape->InheritsFrom("TGeoCone") ) {
926 
927 
928  // can only deal with special case of z-disk and origin in center of cone
929  if( type().isZDisk() ) { // && lo.rho() < epsilon ) {
930 
931  if( lo.rho() > epsilon ) {
932  // move origin to z-axis
933  lo.x() = 0. ;
934  lo.y() = 0. ;
935  }
936 
937  double zhalf = 0 ;
938  double rmax1 = 0 ;
939  double rmax2 = 0 ;
940  double rmin1 = 0 ;
941  double rmin2 = 0 ;
942  if( shape->InheritsFrom("TGeoTube") ) {
943  TGeoTube* tube = ( TGeoTube* ) shape ;
944  zhalf = tube->GetDZ() ;
945  rmax1 = tube->GetRmax() ;
946  rmax2 = tube->GetRmax() ;
947  rmin1 = tube->GetRmin() ;
948  rmin2 = tube->GetRmin() ;
949  } else { // shape->InheritsFrom("TGeoCone") )
950  TGeoCone* cone = ( TGeoCone* ) shape ;
951  zhalf = cone->GetDZ() ;
952  rmax1 = cone->GetRmax1() ;
953  rmax2 = cone->GetRmax2() ;
954  rmin1 = cone->GetRmin1() ;
955  rmin2 = cone->GetRmin2() ;
956  }
957 
958  // two circles around origin
959  // get radii at position of plane
960  double r0 = rmin1 + ( rmin2 - rmin1 ) / ( 2. * zhalf ) * ( zhalf + lo.z() ) ;
961  double r1 = rmax1 + ( rmax2 - rmax1 ) / ( 2. * zhalf ) * ( zhalf + lo.z() ) ;
962 
963 
964  unsigned n = nMax / 4 ;
965  double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
966 
967  for( unsigned i = 0 ; i < n ; ++i ) {
968 
969  Vector3D rv00( r0*sin( i *dPhi ) , r0*cos( i *dPhi ) , 0. ) ;
970  Vector3D rv01( r0*sin( (i+1)*dPhi ) , r0*cos( (i+1)*dPhi ) , 0. ) ;
971 
972  Vector3D rv10( r1*sin( i *dPhi ) , r1*cos( i *dPhi ) , 0. ) ;
973  Vector3D rv11( r1*sin( (i+1)*dPhi ) , r1*cos( (i+1)*dPhi ) , 0. ) ;
974 
975 
976  Vector3D pl0 = lo + rv00 ;
977  Vector3D pl1 = lo + rv01 ;
978 
979  Vector3D pl2 = lo + rv10 ;
980  Vector3D pl3 = lo + rv11 ;
981 
982 
983  Vector3D pg0,pg1,pg2,pg3 ;
984 
985  _wtM->LocalToMaster( pl0, pg0.array() ) ;
986  _wtM->LocalToMaster( pl1, pg1.array() ) ;
987  _wtM->LocalToMaster( pl2, pg2.array() ) ;
988  _wtM->LocalToMaster( pl3, pg3.array() ) ;
989 
990  lines.emplace_back( pg0, pg1 );
991  lines.emplace_back( pg2, pg3 );
992  }
993 
994  //add some vertical and horizontal lines so that the disc is seen in the rho-z projection
995 
996  n = 4 ; dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
997 
998  for( unsigned i = 0 ; i < n ; ++i ) {
999 
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 
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 
1246 
1247  Vector3D lp;
1248  _wtM->MasterToLocal( point , lp.array() ) ;
1249 
1250  return _volSurf.globalToLocal( lp ) ;
1251  }
1252 
1253 
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 
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 
1304 
1305 
1306  //================================================================================================================
1307 
1308  } // namespace
1309 } // namespace
1310 
1311 
dd4hep::rec::VolSurface::v
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:250
dd4hep::rec::VolSurfaceBase::_type
SurfaceType _type
Definition: Surface.h:48
dd4hep::rec::Vector3D::fill
const Vector3D & fill(const T &v)
fill vector from arbitrary class that defines operator[]
Definition: Vector3D.h:82
dd4hep::rec::VolSurfaceBase::length_along_v
virtual double length_along_v() const
Definition: Surface.cpp:142
dd4hep::rec::Surface::initialize
void initialize()
Definition: Surface.cpp:719
dd4hep::rec::Vector3D::rho
double rho() const
Definition: Vector3D.h:148
dd4hep::rec::Surface::localToGlobal
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:682
dd4hep::rec::Vector3D::z
double z() const
Definition: Vector3D.h:109
dd4hep::rec::VolCylinderImpl::localToGlobal
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:322
dd4hep::rec::SurfaceType::checkOrthogonalToZ
bool checkOrthogonalToZ(const ISurface &surf, double epsilon=1.e-6) const
Definition: ISurface.h:269
dd4hep::rec::VolSurface::type
virtual const SurfaceType & type() const
Definition: Surface.cpp:248
dd4hep::rec::SurfaceType::OrthogonalToZ
@ OrthogonalToZ
Definition: ISurface.h:150
dd4hep::rec::VolCylinderImpl::VolCylinderImpl
VolCylinderImpl()
default c'tor
Definition: Surface.h:366
dd4hep::rec::VolConeImpl::_tanTheta
double _tanTheta
Definition: Surface.h:407
v
View * v
Definition: MultiView.cpp:28
dd4hep::rec::Surface::_u
Vector3D _u
Definition: Surface.h:509
dd4hep::rec::Vector3D::array
double * array()
direct access to data as double* - allows modification
Definition: Vector3D.h:216
dd4hep::rec::Vector3D
Definition: Vector3D.h:32
dd4hep::rec::Surface::length_along_v
virtual double length_along_v() const
Definition: Surface.cpp:622
dd4hep::rec::ConeSurface::z1
virtual double z1() const
the end z of the cone
Definition: Surface.cpp:1295
dd4hep::rec::Surface::_det
DetElement _det
Definition: Surface.h:502
dd4hep::rec::VolSurface::getLines
virtual std::vector< std::pair< Vector3D, Vector3D > > getLines(unsigned nMax=100)
Definition: Surface.cpp:265
dd4hep::rec::VolSurfaceBase::normal
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Access to the normal direction at the given point.
Definition: Surface.cpp:47
dd4hep::rec::Surface::_wtM
std::unique_ptr< TGeoMatrix > _wtM
Definition: Surface.h:504
dd4hep::rec::MaterialManager::createAveragedMaterial
MaterialData createAveragedMaterial(const MaterialVec &materials)
Definition: MaterialManager.cpp:196
M_PI
#define M_PI
Definition: Handle.h:33
dd4hep::rec::SurfaceType::Plane
@ Plane
Definition: ISurface.h:146
dd4hep::rec::findVolume
bool findVolume(PlacedVolume pv, Volume theVol, std::list< PlacedVolume > &volList)
Definition: Surface.cpp:539
dd4hep::PlacedVolume
Handle class holding a placed volume (also called physical volume)
Definition: Volumes.h:173
dd4hep::DetElement::placement
PlacedVolume placement() const
Access to the physical volume of this detector element.
Definition: DetElement.cpp:320
dd4hep::rec::VolSurface::localToGlobal
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:254
dd4hep::rec::VolSurface::length_along_v
virtual double length_along_v() const
Definition: Surface.cpp:260
dd4hep::rec::VolConeImpl::localToGlobal
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:426
dd4hep::rec::MaterialVec
std::vector< std::pair< Material, double > > MaterialVec
Definition: MaterialManager.h:30
dd4hep::rec::Vector3D::cylindrical
static Cylindrical cylindrical()
Definition: Vector3D.h:278
dd4hep::rec::VolSurface::volume
Volume volume() const
the volume to which this surface is attached.
Definition: Surface.h:214
dd4hep::rec::VolSurface::u
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:249
DetectorInterna.h
dd4hep::rec::Surface::_type
SurfaceType _type
Definition: Surface.h:508
dd4hep::rec::VolSurfaceBase::distance
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:205
dd4hep::rec::VolSurfaceBase::id
virtual long64 id() const
The id of this surface.
Definition: Surface.cpp:42
dd4hep::rec::VolSurface::outerThickness
virtual double outerThickness() const
Definition: Surface.cpp:258
dd4hep::rec::Vector3D::x
double x() const
Definition: Vector3D.h:103
dd4hep::DetElement::volumeID
VolumeID volumeID() const
The cached VolumeID of this subdetector element.
Definition: DetElement.cpp:343
dd4hep::Handle::name
const char * name() const
Access the object name (or "" if not supported by the object)
dd4hep::rec::Surface::id
virtual long64 id() const
The id of this surface - corresponds to DetElement id.
Definition: Surface.cpp:611
dd4hep::rec::SurfaceType::checkParallelToZ
bool checkParallelToZ(const ISurface &surf, double epsilon=1.e-6) const
Definition: ISurface.h:252
dd4hep::rec::Vector3D::phi
double phi() const
Definition: Vector3D.h:142
dd4hep::rec::VolSurfaceBase::u
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:45
dd4hep::rec::VolSurfaceBase::innerThickness
virtual double innerThickness() const
Definition: Surface.cpp:75
dd4hep::rec::ConeSurface::center
virtual Vector3D center() const
the center of the cone
Definition: Surface.cpp:1303
dd4hep::rec::CylinderSurface::localToGlobal
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:1254
dd4hep::rec::VolSurfaceBase::_v
Vector3D _v
Definition: Surface.h:50
epsilon
const double epsilon
Definition: test_cellid_position_converter.cpp:42
dd4hep::rec::Surface::outerMaterial
virtual const IMaterial & outerMaterial() const
Access to the material in direction of the normal.
Definition: Surface.cpp:645
dd4hep::rec::VolCylinderImpl::globalToLocal
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:310
dd4hep::rec::VolSurfaceBase::setNormal
virtual void setNormal(const Vector3D &n)
setter for daughter classes
Definition: Surface.cpp:39
dd4hep::rec::VolSurfaceList::~VolSurfaceList
virtual ~VolSurfaceList()
Definition: Surface.cpp:511
dd4hep::rec::VolCylinderImpl::normal
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:304
dd4hep::rec::Surface::innerThickness
virtual double innerThickness() const
Definition: Surface.cpp:619
dd4hep::rec::Surface::_n
Vector3D _n
Definition: Surface.h:511
dd4hep::rec::SurfaceList::_isOwner
bool _isOwner
Definition: Surface.h:684
dd4hep::rec::SurfaceType::isCone
bool isCone() const
true if this a conical surface
Definition: ISurface.h:211
Surface.h
dd4hep::rec::VolSurface::innerMaterial
virtual const IMaterial & innerMaterial() const
Access to the material in opposite direction of the normal.
Definition: Surface.cpp:255
dd4hep::rec::Surface::globalToLocal
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:665
dd4hep::rec::VolSurfaceBase::v
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:46
dd4hep::rec::Vector2D::v
double v() const
Definition: Vector2D.h:36
dd4hep::rec::volSurfaceList
VolSurfaceList * volSurfaceList(DetElement &det)
Definition: Surface.cpp:528
dd4hep::rec::VolSurfaceBase::outerMaterial
virtual const IMaterial & outerMaterial() const
Access to the material in direction of the normal.
Definition: Surface.cpp:74
dd4hep::rec::VolConeImpl::VolConeImpl
VolConeImpl()
default c'tor
Definition: Surface.h:412
dd4hep::DetElement
Handle class describing a detector element.
Definition: DetElement.h:188
dd4hep::rec::VolSurfaceBase::insideBounds
virtual bool insideBounds(const Vector3D &point, double epsilon=1e-4) const
Checks if the given point lies within the surface.
Definition: Surface.cpp:208
dd4hep::Volume
Handle class holding a placed volume (also called physical volume)
Definition: Volumes.h:378
dd4hep::rec::VolConeImpl::globalToLocal
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:411
dd4hep::rec::VolSurface::id
virtual long64 id() const
The id of this surface - always 0 for VolSurfaces.
Definition: Surface.cpp:247
dd4hep::rec::VolSurfaceBase::localToGlobal
virtual Vector3D localToGlobal(const Vector2D &point) const
Definition: Surface.cpp:66
dd4hep::rec::VolSurfaceBase::globalToLocal
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:50
dd4hep::rec::VolSurface
Definition: Surface.h:181
dd4hep::rec::VolSurfaceBase::setV
virtual void setV(const Vector3D &v)
setter for daughter classes
Definition: Surface.cpp:38
dd4hep::rec::VolConeImpl::normal
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:405
dd4hep::Alignment::worldTransformation
const TGeoHMatrix & worldTransformation() const
Create cached matrix to transform to world coordinates.
Definition: Alignments.cpp:68
dd4hep::rec::Surface::type
virtual const SurfaceType & type() const
Definition: Surface.cpp:613
dd4hep::rec::long64
long long int long64
Definition: ISurface.h:26
dd4hep::rec::Vector3D::unit
Vector3D unit() const
Definition: Vector3D.h:199
dd4hep::rec::VolSurface::distance
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:261
dd4hep::rec::Surface::distance
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:701
dd4hep::rec::VolSurface::origin
virtual const Vector3D & origin() const
Definition: Surface.cpp:252
dd4hep::rec::VolSurfaceList
Definition: Surface.h:300
dd4hep::rec::Vector2D::u
double u() const
Definition: Vector2D.h:34
dd4hep::rec::Surface::_id
long64 _id
Definition: Surface.h:506
dd4hep::rec::VolSurface::length_along_u
virtual double length_along_u() const
Definition: Surface.cpp:259
dd4hep::rec::CylinderSurface::radius
virtual double radius() const
the radius of the cylinder (rho of the origin vector)
Definition: Surface.cpp:1263
dd4hep::rec::Surface::Surface
Surface()=delete
default c'tor etc. removed
dd4hep::rec::Surface::length_along_u
virtual double length_along_u() const
Definition: Surface.cpp:621
dd4hep::rec::VolConeImpl::_zt1
double _zt1
Definition: Surface.h:406
dd4hep::rec::VolSurfaceBase::outerThickness
virtual double outerThickness() const
Definition: Surface.cpp:76
dd4hep::rec::SurfaceType::ParallelToZ
@ ParallelToZ
Definition: ISurface.h:149
dd4hep::rec::MaterialManager
Definition: MaterialManager.h:41
dd4hep::rec::VolConeImpl::u
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:398
dd4hep::rec::SurfaceList::~SurfaceList
virtual ~SurfaceList()
d'tor deletes all owned surfaces
Definition: Surface.cpp:519
dd4hep::rec::Surface::volumeOrigin
virtual Vector3D volumeOrigin() const
Definition: Surface.cpp:689
dd4hep::Alignment
Main handle class to hold an alignment object.
Definition: Alignments.h:115
dd4hep::DetElement::nominal
Alignment nominal() const
Access to the constant ideal (nominal) alignment information.
Definition: DetElement.cpp:184
dd4hep::rec::Vector2D
Definition: Vector2D.h:24
dd4hep::rec::Surface::getLines
virtual std::vector< std::pair< Vector3D, Vector3D > > getLines(unsigned nMax=100)
Definition: Surface.cpp:833
dd4hep::rec::VolConeImpl::_zt0
double _zt0
Definition: Surface.h:405
dd4hep::rec::VolConeImpl::_ztip
double _ztip
Definition: Surface.h:404
dd4hep::rec::VolSurfaceBase::origin
virtual const Vector3D & origin() const
Definition: Surface.cpp:48
dd4hep::rec::VolSurfaceBase::type
virtual const SurfaceType & type() const
Definition: Surface.cpp:44
dd4hep::rec::VolSurfaceBase::innerMaterial
virtual const IMaterial & innerMaterial() const
Access to the material in opposite direction of the normal.
Definition: Surface.cpp:73
dd4hep::rec::SurfaceType
Definition: ISurface.h:140
dd4hep::rec::VolSurface::setOuterMaterial
void setOuterMaterial(const IMaterial &mat)
set the outer Materal
Definition: Surface.h:282
dd4hep::rec::Surface::volume
Volume volume() const
The volume that has the surface attached.
Definition: Surface.h:537
dd4hep::rec::Vector3D::cross
Vector3D cross(const Vector3D &v) const
Definition: Vector3D.h:191
dd4hep::rec::Vector3D::const_array
const double * const_array() const
direct access to data as const double*
Definition: Vector3D.h:211
dd4hep::rec::ConeSurface::z0
virtual double z0() const
the start z of the cone
Definition: Surface.cpp:1287
dd4hep::rec::MaterialManager::materialsBetween
const MaterialVec & materialsBetween(const Vector3D &p0, const Vector3D &p1, double epsilon=1e-4)
Definition: MaterialManager.cpp:40
dd4hep::rec::Surface::_volSurf
VolSurface _volSurf
Definition: Surface.h:503
dd4hep::rec::VolSurfaceBase
Definition: Surface.h:43
dd4hep::rec::VolSurface::insideBounds
virtual bool insideBounds(const Vector3D &point, double epsilon=1e-4) const
Checks if the given point lies within the surface.
Definition: Surface.cpp:262
dd4hep::rec::ConeSurface::radius0
virtual double radius0() const
the start radius of the cone
Definition: Surface.cpp:1271
dd4hep::rec::VolSurface::innerThickness
virtual double innerThickness() const
Definition: Surface.cpp:257
Memory.h
dd4hep::rec::Vector3D::y
double y() const
Definition: Vector3D.h:106
dd4hep::rec::VolSurface::setInnerMaterial
void setInnerMaterial(const IMaterial &mat)
set the innerMaterial
Definition: Surface.h:279
dd4hep::rec::Surface::innerMaterial
virtual const IMaterial & innerMaterial() const
Access to the material in opposite direction of the normal.
Definition: Surface.cpp:626
dd4hep::rec::VolSurfaceBase::_n
Vector3D _n
Definition: Surface.h:51
dd4hep::rec::SurfaceType::Cylinder
@ Cylinder
Definition: ISurface.h:145
dd4hep::rec::Surface::normal
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Access to the normal direction at the given point.
Definition: Surface.cpp:617
dd4hep::rec::VolSurfaceBase::setU
virtual void setU(const Vector3D &u)
setter for daughter classes
Definition: Surface.cpp:37
dd4hep::rec::VolSurfaceBase::getLines
virtual std::vector< std::pair< Vector3D, Vector3D > > getLines(unsigned nMax=100)
Definition: Surface.cpp:238
dd4hep::rec::Vector3D::spherical
static Spherical spherical()
Definition: Vector3D.h:279
MaterialManager.h
dist
double dist(const Position &p0, const Position &p1)
Definition: test_cellid_position_converter.cpp:46
dd4hep::Handle::ptr
T * ptr() const
Access to the held object.
Definition: Handle.h:153
dd4hep::rec::VolCylinderImpl::u
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:297
dd4hep::rec::Surface::_v
Vector3D _v
Definition: Surface.h:510
dd4hep::rec::IMaterial::Z
virtual double Z() const =0
averaged proton number
dd4hep
Namespace for the AIDA detector description toolkit.
Definition: AlignmentsCalib.h:28
dd4hep::rec::VolSurfaceBase::length_along_u
virtual double length_along_u() const
Definition: Surface.cpp:79
dd4hep::rec::VolSurface::globalToLocal
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:253
dd4hep::rec::CylinderSurface::normal
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:1237
dd4hep::rec::Surface::origin
virtual const Vector3D & origin() const
Definition: Surface.cpp:618
det
DetElement::Object * det
Definition: AlignmentsCalculator.cpp:66
dd4hep::rec::ConeSurface::radius1
virtual double radius1() const
the end radius of the cone
Definition: Surface.cpp:1279
dd4hep::PlacedVolume::volume
Volume volume() const
Logical volume of this placement.
Definition: Volumes.cpp:452
dd4hep::PlacedVolume::data
Object * data() const
Check if placement is properly instrumented.
Definition: Volumes.cpp:431
dd4hep::rec::Surface::v
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:616
dd4hep::rec::VolSurfaceBase::setOrigin
virtual void setOrigin(const Vector3D &o)
setter for daughter classes
Definition: Surface.cpp:40
dd4hep::rec::SurfaceType::setProperty
void setProperty(unsigned prop, bool val=true)
set the given peorperty
Definition: ISurface.h:196
dd4hep::rec::Surface::insideBounds
virtual bool insideBounds(const Vector3D &point, double epsilon=1.e-4) const
Checks if the given point lies within the surface.
Definition: Surface.cpp:710
dd4hep::rec::VolConeImpl::distance
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:442
dd4hep::rec::CylinderSurface::u
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:1220
dd4hep::rec::IMaterial
Definition: IMaterial.h:28
dd4hep::rec::VolConeImpl::v
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:392
dd4hep::rec::CylinderSurface::globalToLocal
virtual Vector2D globalToLocal(const Vector3D &point) const
Definition: Surface.cpp:1245
dd4hep::rec::Surface::_o
Vector3D _o
Definition: Surface.h:512
dd4hep::rec::VolCylinderImpl::distance
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:335
dd4hep::rec::VolPlaneImpl::distance
virtual double distance(const Vector3D &point) const
Definition: Surface.cpp:272
dd4hep::rec::CylinderSurface::center
virtual Vector3D center() const
the center of the cylinder
Definition: Surface.cpp:1265
dd4hep::rec::Surface::outerThickness
virtual double outerThickness() const
Definition: Surface.cpp:620
dd4hep::rec::CylinderSurface::v
virtual Vector3D v(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:1229
dd4hep::rec::VolConeImpl::getLines
virtual std::vector< std::pair< Vector3D, Vector3D > > getLines(unsigned nMax=100)
create outer bounding lines for the given symmetry of the polyhedron
Definition: Surface.cpp:467
dd4hep::rec::Surface::u
virtual Vector3D u(const Vector3D &point=Vector3D()) const
Definition: Surface.cpp:615
dd4hep::rec::SurfaceType::Cone
@ Cone
Definition: ISurface.h:153
dd4hep::rec::VolSurface::normal
virtual Vector3D normal(const Vector3D &point=Vector3D()) const
Access to the normal direction at the given point.
Definition: Surface.cpp:251
dd4hep::rec::VolSurface::outerMaterial
virtual const IMaterial & outerMaterial() const
Access to the material in direction of the normal.
Definition: Surface.cpp:256
dd4hep::dd4hep_ptr
Out version of the std auto_ptr implementation base either on auto_ptr or unique_ptr.
Definition: Memory.h:46
dd4hep::rec::Vector3D::theta
double theta() const
Definition: Vector3D.h:179