23 #include "TGeoMatrix.h"
24 #include "TGeoShape.h"
25 #include "TRotation.h"
32 using namespace detail ;
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 ;
63 return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
68 Vector3D g = origin() + point[0] * u() + point[1] *
v() ;
82 const Vector3D& u_val = this->u( o ) ;
93 if( volume()->GetShape()->Contains( o.
const_array() ) ){
95 dist_p = volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o.
const_array() ) ,
97 dist_m = volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o.
const_array() ) ,
98 const_cast<double*
> ( um.
array() ) ) ;
109 dist_p = volume()->GetShape()->DistFromOutside(
const_cast<double*
> ( o.
const_array() ) ,
111 dist_m = volume()->GetShape()->DistFromOutside(
const_cast<double*
> ( o.
const_array() ) ,
112 const_cast<double*
> ( um.
array() ) ) ;
122 Vector3D o_1 = this->origin() + dist_p * u_val ;
123 Vector3D o_2 = this->origin() + dist_m * um ;
125 dist_p += volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o_1.
const_array() ) ,
126 const_cast<double*
> ( u_val.const_array() ) ) ;
128 dist_m += volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o_2.
const_array() ) ,
129 const_cast<double*
> ( um.array() ) ) ;
137 return dist_p + dist_m ;
144 const Vector3D& o = this->origin() ;
156 if( volume()->GetShape()->Contains( o.
const_array() ) ){
158 dist_p = volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o.
const_array() ) ,
160 dist_m = volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o.
const_array() ) ,
161 const_cast<double*
> ( vm.
array() ) ) ;
172 dist_p = volume()->GetShape()->DistFromOutside(
const_cast<double*
> ( o.
const_array() ) ,
174 dist_m = volume()->GetShape()->DistFromOutside(
const_cast<double*
> ( o.
const_array() ) ,
175 const_cast<double*
> ( vm.
array() ) ) ;
185 Vector3D o_1 = this->origin() + dist_p * v_val ;
186 Vector3D o_2 = this->origin() + dist_m * vm ;
188 dist_p += volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o_1.
const_array() ) ,
189 const_cast<double*
> ( v_val.const_array() ) ) ;
191 dist_m += volume()->GetShape()->DistFromInside(
const_cast<double*
> ( o_2.
const_array() ) ,
192 const_cast<double*
> ( vm.array() ) ) ;
200 return dist_p + dist_m ;
212 bool inShape = ( type().isUnbounded() ? true : volume()->GetShape()->Contains( point.
const_array() ) ) ;
214 double dist = std::abs ( distance( point ) ) ;
216 std::cout <<
" ** Surface::insideBound( " << point <<
" ) - distance = " <<
dist
217 <<
" origin = " << origin() <<
" normal = " << normal()
218 <<
" p * n = " << point * normal()
219 <<
" isInShape : " << inShape << std::endl ;
224 if( type().isUnbounded() ){
226 return std::abs ( distance( point ) ) <
epsilon ;
230 return ( std::abs ( distance( point ) ) <
epsilon && volume()->GetShape()->Contains( point.
const_array() ) ) ;
240 std::vector< std::pair<Vector3D, Vector3D> > lines ;
263 return _surf->insideBounds( point,
epsilon ) ;
266 return _surf->getLines(nMax) ;
273 return ( point - origin() ) * normal() ;
278 double thickness_inner ,
double thickness_outer,
Vector3D o ) :
315 while( phi < -
M_PI ) phi += 2.*
M_PI ;
316 while( phi >
M_PI ) phi -= 2.*
M_PI ;
324 double z = point.
v() +
origin().
z() ;
327 while( phi < -
M_PI ) phi += 2.*
M_PI ;
328 while( phi >
M_PI ) phi -= 2.*
M_PI ;
342 double thickness_inner ,
double thickness_outer,
Vector3D v_val,
Vector3D o_val ) :
346 Vector3D o_rphi( o_val.
x() , o_val.
y() , 0. ) ;
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 ;
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() ) ;
359 double theta = v_val.
theta() ;
371 _ztip = o_val.
z() - tipoffset ;
373 double dist_p = vol->GetShape()->DistFromInside(
const_cast<double*
> ( o_val.
const_array() ) ,
376 double dist_m = vol->GetShape()->DistFromInside(
const_cast<double*
> ( o_val.
const_array() ) ,
377 const_cast<double*
> ( vm.
array() ) ) ;
379 double costh = std::cos( theta) ;
380 _zt0 = tipoffset - dist_m * costh ;
381 _zt1 = tipoffset + dist_p * costh ;
402 return av.
cross( n ) ;
416 while( phi < -
M_PI ) phi += 2.*
M_PI ;
417 while( phi >
M_PI ) phi -= 2.*
M_PI ;
434 while( phi < -
M_PI ) phi += 2.*
M_PI ;
435 while( phi >
M_PI ) phi -= 2.*
M_PI ;
460 double zp = point.
z() -
_ztip ;
462 return r * std::cos(
_v.
theta() ) ;
469 std::vector< std::pair<Vector3D, Vector3D> > lines ;
471 lines.reserve( nMax ) ;
473 double theta =
v().
theta() ;
476 Vector3D zv( 0. , 0. , half_length ) ;
478 double dr = half_length * tan( theta ) ;
484 unsigned n = nMax / 4 ;
485 double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
487 for(
unsigned i = 0 ; i < n ; ++i ) {
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. ) ;
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. ) ;
500 lines.emplace_back( pl0, pl1 );
501 lines.emplace_back( pl1, pl2 );
502 lines.emplace_back( pl2, pl3 );
503 lines.emplace_back( pl3, pl0 );
522 std::for_each(begin(), end(), detail::deleteObject<ISurface>);
542 volList.emplace_back( pv ) ;
561 const TGeoNode* node = pv.
ptr();
567 throw std::runtime_error(
"*** findVolume: Invalid placement: - node pointer Null ! " + std::string( pv.
name() ) );
573 for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
575 TGeoNode* daughter = node->GetDaughter(idau);
578 if ( !placement.
data() ) {
579 throw std::runtime_error(
"*** findVolume: Invalid not instrumented placement:"+std::string(daughter->GetName())
580 +
" [Internal error -- bad detector constructor]");
585 if(
findVolume( pv_dau , theVol , volList ) ) {
606 _wtM() , _id( 0) , _type( _volSurf.type() ) {
630 if( ! ( mat.
Z() > 0 ) ) {
640 materials[0].first ) ;
649 if( ! ( mat.
Z() > 0 ) ) {
659 materials[0].first ) ;
672 double uv =
u() *
v() ;
675 double uup =
u() * uprime ;
676 double vvp =
v() * vprime ;
678 return Vector2D( p*uprime / uup , p*vprime / vvp ) ;
704 _wtM->MasterToLocal( point , pa ) ;
713 _wtM->MasterToLocal( point , pa ) ;
722 std::list< PlacedVolume > pVList ;
728 std::stringstream sst ; sst <<
" ***** ERROR: Volume " << theVol.
name() <<
" not found for DetElement " <<
_det.
name() <<
" with surface " ;
729 throw std::runtime_error( sst.str() ) ;
738 for( std::list<PlacedVolume>::iterator it= pVList.begin(), n = pVList.end() ; it != n ; ++it ){
740 TGeoMatrix* m = pv->GetMatrix();
741 std::cout <<
" +++ matrix for placed volume : " << std::endl ;
749 std::unique_ptr<TGeoHMatrix> wtI(
new TGeoHMatrix( wm ) ) ;
753 for( std::list<PlacedVolume>::iterator it = ++( pVList.begin() ) , n = pVList.end() ; it != n ; ++it ){
756 TGeoMatrix* m = pvol->GetMatrix();
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 ???
770 _wtM = std::move(wtI);
774 _wtM = std::move(wtI);
780 double ua[3], va[3], na[3], oa[3] ;
836 const static double epsilon = 1e-6 ;
838 std::vector< std::pair<Vector3D, Vector3D> > lines ;
842 const std::vector< std::pair<Vector3D, Vector3D> >& local_lines =
_volSurf.
getLines() ;
844 if( local_lines.size() > 0 ) {
845 unsigned n=local_lines.size() ;
848 for(
unsigned i=0;i<n;++i){
851 _wtM->LocalToMaster( local_lines[i].first , av.
array() ) ;
852 _wtM->LocalToMaster( local_lines[i].second , bv.
array() ) ;
854 lines.emplace_back( av, bv );
869 const TGeoShape* shape = vol->GetShape() ;
872 if(
type().isPlane() ) {
874 if( shape->IsA() == TGeoBBox::Class() ) {
876 TGeoBBox* box = ( TGeoBBox* ) shape ;
878 Vector3D boxDim( box->GetDX() , box->GetDY() , box->GetDZ() ) ;
881 bool isYZ = std::fabs( ln.
x() - 1.0 ) <
epsilon ;
882 bool isXZ = std::fabs( ln.
y() - 1.0 ) <
epsilon ;
883 bool isXY = std::fabs( ln.
z() - 1.0 ) <
epsilon ;
886 if( isYZ || isXZ || isXY ) {
897 ubl.
fill( 1., 0., 0. ) ;
898 vbl.
fill( 0., 0., 1. ) ;
904 ubl.
fill( 1., 0., 0. ) ;
905 vbl.
fill( 0., 1., 0. ) ;
912 _wtM->LocalToMasterVect( ubl , ub.
array() ) ;
913 _wtM->LocalToMasterVect( vbl , vb.
array() ) ;
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 );
925 }
else if( shape->InheritsFrom(
"TGeoTube") || shape->InheritsFrom(
"TGeoCone") ) {
929 if(
type().isZDisk() ) {
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() ;
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() ;
960 double r0 = rmin1 + ( rmin2 - rmin1 ) / ( 2. * zhalf ) * ( zhalf + lo.
z() ) ;
961 double r1 = rmax1 + ( rmax2 - rmax1 ) / ( 2. * zhalf ) * ( zhalf + lo.
z() ) ;
964 unsigned n = nMax / 4 ;
965 double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
967 for(
unsigned i = 0 ; i < n ; ++i ) {
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. ) ;
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. ) ;
985 _wtM->LocalToMaster( pl0, pg0.
array() ) ;
986 _wtM->LocalToMaster( pl1, pg1.
array() ) ;
987 _wtM->LocalToMaster( pl2, pg2.
array() ) ;
988 _wtM->LocalToMaster( pl3, pg3.
array() ) ;
990 lines.emplace_back( pg0, pg1 );
991 lines.emplace_back( pg2, pg3 );
996 n = 4 ; dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
998 for(
unsigned i = 0 ; i < n ; ++i ) {
1000 Vector3D rv0( r0*sin( i * dPhi ) , r0*cos( i * dPhi ) , 0. ) ;
1001 Vector3D rv1( r1*sin( i * dPhi ) , r1*cos( i * dPhi ) , 0. ) ;
1008 _wtM->LocalToMaster( pl0, pg0.
array() ) ;
1009 _wtM->LocalToMaster( pl1, pg1.
array() ) ;
1011 lines.emplace_back(pg0, pg1);
1019 else if(shape->IsA() == TGeoTrap::Class()) {
1020 TGeoTrap* trapezoid = ( TGeoTrap* ) shape;
1022 double dx1 = trapezoid->GetBl1();
1023 double dx2 = trapezoid->GetTl1();
1024 double dz = trapezoid->GetH1();
1034 _wtM->LocalToMasterVect( ubl , ub.
array() ) ;
1035 _wtM->LocalToMasterVect( vbl , vb.
array() ) ;
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);
1048 else if(shape->IsA() == TGeoTrd1::Class()){
1049 TGeoTrd1* trapezoid = ( TGeoTrd1* ) shape;
1051 double dx1 = trapezoid->GetDx1();
1052 double dx2 = trapezoid->GetDx2();
1054 double dz = trapezoid->GetDz();
1065 _wtM->LocalToMasterVect( ubl , ub.
array() ) ;
1066 _wtM->LocalToMasterVect( vbl , vb.
array() ) ;
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);
1079 else if(shape->IsA() == TGeoTrd2::Class()){
1080 TGeoTrd2* trapezoid = ( TGeoTrd2* ) shape;
1082 double dx1 = trapezoid->GetDx1();
1083 double dx2 = trapezoid->GetDx2();
1086 double dz = trapezoid->GetDz();
1097 _wtM->LocalToMasterVect( ubl , ub.
array() ) ;
1098 _wtM->LocalToMasterVect( vbl , vb.
array() ) ;
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);
1119 lines.reserve( nMax ) ;
1121 double dAlpha = 2.* ROOT::Math::Pi() / double( nMax ) ;
1123 TVector3 norm( ln.
x() , ln.
y() , ln.
z() ) ;
1128 for(
unsigned i=0 ; i< nMax ; ++i ){
1130 double alpha = double(i) * dAlpha ;
1132 TVector3 vec( lu.
x() , lu.
y() , lu.
z() ) ;
1135 rot.Rotate( alpha , norm );
1137 TVector3 vecR = rot * vec ;
1140 luRot.
fill( vecR ) ;
1142 double dist = shape->DistFromInside(
const_cast<double*
> (lo.
const_array()) ,
const_cast<double*
> (luRot.
const_array()) , 3, 0.1 ) ;
1149 _wtM->LocalToMaster( lp , gp.
array() ) ;
1160 lines.emplace_back(previous, gp);
1166 lines.emplace_back(previous, first);
1169 }
else if(
type().isCylinder() ) {
1171 if( shape->InheritsFrom(
"TGeoTube") || shape->InheritsFrom(
"TGeoCone") ) {
1173 lines.reserve( nMax ) ;
1175 TGeoBBox* tube = ( TGeoBBox* ) shape ;
1177 double zHalf = tube->GetDZ() ;
1181 double r = lo.
rho() ;
1184 unsigned n = nMax / 4 ;
1185 double dPhi = 2.* ROOT::Math::Pi() / double( n ) ;
1187 for(
unsigned i = 0 ; i < n ; ++i ) {
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. ) ;
1201 _wtM->LocalToMaster( pl0, pg0.
array() ) ;
1202 _wtM->LocalToMaster( pl1, pg1.
array() ) ;
1203 _wtM->LocalToMaster( pl2, pg2.
array() ) ;
1204 _wtM->LocalToMaster( pl3, pg3.
array() ) ;
1206 lines.emplace_back( pg0, pg1 );
1207 lines.emplace_back( pg1, pg2 );
1208 lines.emplace_back( pg2, pg3 );
1209 lines.emplace_back( pg3, pg0 );
1223 _wtM->MasterToLocal( point , lp.
array() ) ;
1225 _wtM->LocalToMasterVect( lu , u_val.
array() ) ;
1231 _wtM->MasterToLocal( point , lp.
array() ) ;
1233 _wtM->LocalToMasterVect( lv , v_val.
array() ) ;
1239 _wtM->MasterToLocal( point , lp.
array() ) ;
1241 _wtM->LocalToMasterVect( ln , n.
array() ) ;
1248 _wtM->MasterToLocal( point , lp.
array() ) ;
1276 return origin().
rho() - 0.5 * l * tan( theta ) ;
1284 return origin().
rho() + 0.5 * l * tan( theta ) ;
1292 return origin().
z() - 0.5 * l ;
1300 return origin().
z() + 0.5 * l ;