DDTest/src/test_surface.cc

0001 #include "DD4hep/DDTest.h"
0002
0003 #include "DD4hep/Detector.h"
0004 #include "DD4hep/DD4hepUnits.h"
0005 #include "DD4hep/Volumes.h"
0006 #include "DD4hep/DetElement.h"
0007
0008 #include "DDRec/Surface.h"
0009 #include "STR.h"
0010
0011 #include <exception>
0012 #include <iostream>
0013 #include <assert.h>
0014 #include <cmath>
0015
0016 using namespace std ;
0017 using namespace dd4hep ;
0018 using namespace dd4hep::rec ;
0019 using namespace detail;
0020
0021 // this should be the first line in your test
0022 static DDTest test( "surface" ) ;
0023 //=============================================================================
0024
0025 int main(int argc, char** argv ){
0026
0027   test.log( "test units" );
0028
0029   if( argc < 2 ) {
0030     std::cout << " usage:  test_surface units.xml " << std::endl ;
0031     exit(1) ;
0032   }
0033
0034   try{
0035
0036     // ----- write your tests in here -------------------------------------
0037
0038     Detector& description = Detector::getInstance();
0039
0040     description.fromCompact( argv[1] );
0041
0042     // --- test a planar surface
0043     double thick  = 0.005 ;
0044     double width  = 1.0  ;
0045     double length = 10.0 ;
0046
0047     Material    mat    = description.material( "Silicon" );
0048     Box         box   ( thick/2.,  width/2.,  length/2. );
0049     Volume      vol   ( "test_box", box, mat);
0050
0051     Vector3D u( 0. , 1. , 0. ) ;
0052     Vector3D v( 0. , 0. , 1. ) ;
0053     Vector3D n( 1. , 0. , 0. ) ;
0054     Vector3D o( 0. , 0. , 0. ) ;
0055
0056
0057     VolPlane surf( vol , SurfaceType( SurfaceType::Sensitive ), thick/2, thick/2 , u,v,n,o ) ;
0058
0059     // test inside bounds for some points:
0060
0061     test( surf.insideBounds(  Vector3D(  0. , 23. , 42.  )  ) , false , " insideBounds Vector3D(   0. , 23. , 42. )   " ) ;
0062
0063     test( surf.insideBounds(  Vector3D(  0,  .23 ,  .42  )  ) , true , " insideBounds Vector3D(    0,  .23 ,  .42  )   " ) ;
0064
0065     test( surf.insideBounds(  Vector3D(  0.00003 ,  .23 ,  .42  )  ) , true , " insideBounds Vector3D(   0.00003 ,  .23 ,  .42  )   " ) ;
0066
0067
0068     //=============== test global to local ===================
0069
0070     Vector3D point = o + 34.3 * u - 42.7 * v ;
0071
0072     Vector2D lp = surf.globalToLocal( point ) ;
0073     //    std::cout << " --- local coordinates of " << point << " : (" << lp[0] << "," << lp[1] << ")" << std::endl ;
0074     test(  STR( lp[0] ) == STR( 34.3 ) , true , " local u coordinate is 34.4 "  ) ;
0075     test(  STR( lp[1] ) == STR( -42.7 ) , true , " local v coordinate is -42.7 "  ) ;
0076
0077     Vector3D pointPrime = surf.localToGlobal( lp ) ;
0078     test(  pointPrime.isEqual( point ) , true , " point after global to local to global is the same " ) ;
0079
0080     // ----- test with rotated coordinates
0081     Vector3D ur,vr ;
0082     ur.fill( 0. ,  cos( 3.5/180.*M_PI  ) ,  sin( 3.5/180.*M_PI  ) ) ;
0083     vr.fill( 0. , -sin( 3.5/180.*M_PI  ) ,  cos( 3.5/180.*M_PI  ) ) ;
0084     VolPlane surfR( vol , SurfaceType( SurfaceType::Sensitive ), thick/2, thick/2 , ur,vr,n,o ) ;
0085
0086     Vector3D pointR = o + 34.3 * ur - 42.7 * vr ;
0087
0088     lp = surfR.globalToLocal( pointR ) ;
0089     //    std::cout << " --- local coordinates of " << pointR << " : (" << lp[0] << "," << lp[1] << ")" << std::endl ;
0090     test(  STR( lp[0] ) == STR( 34.3 ) , true , " local u coordinate is 34.4 "  ) ;
0091     test(  STR( lp[1] ) == STR( -42.7 ) , true , " local v coordinate is -42.7 "  ) ;
0092
0093     Vector3D pointPrimeR = surfR.localToGlobal( lp ) ;
0094     test(  pointPrimeR.isEqual( pointR ) , true , " point after global to local to global is the same " ) ;
0095
0096     // ----- test with non-orthogonal rotated coordinates
0097     Vector3D vr2 ;
0098     vr2.fill( 0. , -sin( 35./180.*M_PI  ) ,  cos( 35./180.*M_PI  ) ) ;
0099     VolPlane surfR2( vol , SurfaceType( SurfaceType::Sensitive ), thick/2, thick/2 , ur,vr2,n,o ) ;
0100
0101     Vector3D pointR2 = o + 34.3 * ur - 42.7 * vr2 ;
0102
0103     lp = surfR2.globalToLocal( pointR2 ) ;
0104     //    std::cout << " --- local coordinates of " << pointR << " : (" << lp[0] << "," << lp[1] << ")" << std::endl ;
0105     test(  STR( lp[0] ) == STR( 34.3 ) , true , " local u coordinate is 34.4 "  ) ;
0106     test(  STR( lp[1] ) == STR( -42.7 ) , true , " local v coordinate is -42.7 "  ) ;
0107
0108     Vector3D pointPrimeR2 = surfR2.localToGlobal( lp ) ;
0109     test(  pointPrimeR2.isEqual( pointR2 ) , true , " point after global to local to global is the same " ) ;
0110
0111
0112     //==================================================
0113
0114
0115
0116     // --- test MaterialData
0117     MaterialData sm( mat ) ;
0118
0119     // material properies of Si :
0120     test( STR( sm.A() )  , STR( 28.0855 ) , "   MaterialData.A() == 28.0855 " ) ;
0121
0122     test( STR( sm.Z() )  , STR( 14 ) , "   MaterialData.Z() == 14 " ) ;
0123
0124     test( STR( sm.density() )  , STR( 2.33 ) , "   MaterialData.density() == 2.33 " ) ;
0125
0126     test( STR( sm.radiationLength() / dd4hep::mm )  , STR( 93.6612 ) , "   MaterialData.radiationLength() == 93.4961 * mm " ) ;
0127
0128     test( STR( sm.interactionLength() / dd4hep::mm )  , STR( 457.533 ) , "   MaterialData.interactionLength() == 457.532 * mm " ) ;
0129
0130
0131
0132     // test surface type:
0133
0134     test( surf.type().isSensitive() , true , " surface is sensitive " )  ;
0135
0136     test( surf.type().isPlane() , true , " surface is Plane " )  ;
0137
0138     surf.type().checkParallelToZ(  surf ) ;
0139
0140     test( surf.type().isZPlane() , true , " surface is ZPlane " )  ;
0141
0142     test( surf.type().isCylinder() , false , " surface is no Cylinder " )  ;
0143
0144
0145     // --------------------------------------------------------------------
0146
0147     // test a cylindrical surface
0148
0149     thick  = 1.0 ;
0150     length = 100.0 ;
0151     double radius = 42.0  ;
0152
0153     mat    = description.material( "Air" );
0154     Tube    tube  ( radius - thick/2., radius + thick/2. , length/2. );
0155     vol = Volume( "test_tube", tube , mat);
0156
0157     Vector3D o_radius = Vector3D( radius , 0. , 0.  ) ;
0158
0159     VolCylinder surfT( vol , SurfaceType( SurfaceType::Sensitive ), thick/2, thick/2 , o_radius  ) ;
0160
0161     std::cout << " *** cylinder surface : " << surfT << std::endl ;
0162
0163     test( surfT.insideBounds(  Vector3D(  radius * sin(0.75) , radius * cos( 0.75 ) , 49.  )) , true
0164       , " insideBounds Vector3D(  radius * sin(0.75) , radius * cos( 0.75 ) , 49. )  " ) ;
0165
0166     test( surfT.insideBounds(  Vector3D(  radius * sin(0.75) , radius * cos( 0.75 ) , 50.01  )) , false
0167       , " insideBounds Vector3D(  radius * sin(0.75) , radius * cos( 0.75 ) , 50.01 )  " ) ;
0168
0169     test( surfT.insideBounds(  Vector3D(  (radius+0.001) * sin(0.75) , (radius+0.001) * cos( 0.75 ) , 49.  )) , false
0170       , " insideBounds Vector3D(  (radius+0.001) * sin(0.75) , (radius+0.001) * cos( 0.75 ) , 49. )  " ) ;
0171
0172
0173     test( surfT.insideBounds(  Vector3D(  (radius+0.00005) * sin(0.75) , (radius+0.00005) * cos( 0.75 ) , 49.  )) , true
0174       , " insideBounds Vector3D(  (radius+0.00005) * sin(0.75) , (radius+0.00005) * cos( 0.75 ) , 49. )  " ) ;
0175
0176
0177     Vector3D y( 0. , radius , 42 ) ;
0178
0179     Vector3D yn = surfT.normal( y ) ;
0180
0181     bool dummy =   yn.isEqual( Vector3D( 0. , 1. , 0 ) ) ;
0182     test( dummy , true , "  normal at (0.,radius,42) is  Vector3D( 0. , 1. , 0 ) " ) ;
0183
0184     if( ! dummy )
0185       std::cout << " ** yn = " << yn << std::endl ;
0186
0187
0188     Vector3D yv = surfT.u( y ) ;
0189     dummy = yv.isEqual( Vector3D( -1. , 0. , 0 ) ) ;
0190     test( dummy , true , "  u at (0.,radius,42) is  Vector3D( -1. , 0. , 0 ) " ) ;
0191     if( ! dummy )
0192       std::cout << " ** yv = " << yv << std::endl ;
0193
0194
0195     //=============== test global to local ===================
0196
0197     Vector3D pointC( radius , -42.7/radius , 34.3  , Vector3D::cylindrical  )  ;
0198
0199     Vector2D lpC = surfT.globalToLocal( pointC ) ;
0200
0201     // std::cout << " --- local coordinates of " << pointC << " : (" << lpC[0] << "," << lpC[1] << ")" << std::endl ;
0202
0203     test(  STR( lpC[0] ) == STR( -42.7 ) , true , " local u coordinate is -42.7 "  ) ;
0204     test(  STR( lpC[1] ) == STR( 34.3 ) , true  , " local v coordinate is 34.4 "  ) ;
0205
0206     Vector3D pointPrimeC = surfT.localToGlobal( lpC ) ;
0207
0208     // std::cout << " --- global coordinates of point after local to global ( " << pointC.rho() << ", " << pointC.phi() << ", " << pointC.z()
0209     //            << " ) : (" << lpC[0] << "," << lpC[1] << ")" << std::endl ;
0210
0211     test(  pointPrimeC.isEqual( pointC ) , true , " point after global to local to global is the same " ) ;
0212
0213     //========================================================
0214
0215
0216
0217
0218
0219     // test surface type:
0220
0221     test( surfT.type().isSensitive() , true , " surface is sensitive " )  ;
0222
0223     test( surfT.type().isPlane() , false , " surface is no Plane " )  ;
0224
0225
0226     surfT.type().checkParallelToZ(  surfT ) ;
0227
0228     test( surfT.type().isZCylinder() , true , " surface is ZCylinder " )  ;
0229
0230
0231    // --------------------------------------------------------------------
0232
0233
0234   } catch( exception &e ){
0235     //} catch( ... ){
0236
0237     test.log( e.what() );
0238     test.error( "exception occurred" );
0239   }
0240
0241   return 0;
0242 }
0243
0244 //=============================================================================