00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "G4FCylindricalSurface.hh"
00037 #include "G4PhysicalConstants.hh"
00038 #include "G4Sort.hh"
00039
00040
00041 G4FCylindricalSurface::G4FCylindricalSurface()
00042 : radius(0.), length(1.)
00043 {
00044 }
00045
00046
00047 G4FCylindricalSurface::~G4FCylindricalSurface()
00048 {
00049 }
00050
00051
00052 G4FCylindricalSurface::G4FCylindricalSurface( const G4Point3D& o,
00053 const G4Vector3D& a,
00054 G4double r,
00055 G4double l
00056 )
00057 {
00058
00059
00060 G4Vector3D dir(1,1,1);
00061 Position.Init(dir, a, o);
00062
00063 origin = o;
00064 radius = r;
00065
00066
00067 if ( l >= 0.0 )
00068 length = l;
00069 else
00070 {
00071 std::ostringstream message;
00072 message << "Negative length." << G4endl
00073 << "Default length of 0.0 is used.";
00074 G4Exception("G4FCylindricalSurface::G4FCylindricalSurface()",
00075 "GeomSolids1001", JustWarning, message);
00076
00077 length = 0.0;
00078 }
00079
00080
00081 if ( r >= 0.0 )
00082 radius = r;
00083 else
00084 {
00085 std::ostringstream message;
00086 message << "Negative radius." << G4endl
00087 << "Default value of 0.0 is used.";
00088 G4Exception("G4FCylindricalSurface::G4FCylindricalSurface()",
00089 "GeomSolids1001", JustWarning, message);
00090
00091 radius = 0.0;
00092 }
00093 }
00094
00095
00096 const char* G4FCylindricalSurface::NameOf() const
00097 {
00098 return "G4FCylindricalSurface";
00099 }
00100
00101
00102 void G4FCylindricalSurface::PrintOn( std::ostream& os ) const
00103 {
00104 os << "G4FCylindricalSurface with origin: " << origin << "\t"
00105 << "and axis: " << Position.GetAxis() << "\n"
00106 << "\t radius: " << radius << "\t and length: "
00107 << length << "\n";
00108 }
00109
00110
00111 G4double G4FCylindricalSurface::Area() const
00112 {
00113 return ( 2.0 * pi * radius * length );
00114 }
00115
00116
00117
00118
00119 void G4FCylindricalSurface::CalcBBox()
00120 {
00121
00122
00123
00124
00125 G4Point3D Max = G4Point3D(-PINFINITY);
00126 G4Point3D Min = G4Point3D( PINFINITY);
00127
00128 G4Point3D Tmp;
00129 G4Point3D Origin = Position.GetLocation();
00130 G4Point3D EndOrigin = G4Point3D( Origin + (length*Position.GetAxis()) );
00131 G4Point3D Radius(radius, radius, 0);
00132
00133
00134 G4Point3D Tolerance(kCarTolerance, kCarTolerance, kCarTolerance);
00135 G4Point3D BoxMin(Origin-Tolerance);
00136 G4Point3D BoxMax(Origin+Tolerance);
00137
00138 bbox = new G4BoundingBox3D();
00139 bbox->Init(BoxMin, BoxMax);
00140
00141
00142 Tmp = (Origin - Radius);
00143 bbox->Extend(Tmp);
00144
00145 Tmp = Origin + Radius;
00146 bbox->Extend(Tmp);
00147
00148 Tmp = EndOrigin - Radius;
00149 bbox->Extend(Tmp);
00150
00151 Tmp = EndOrigin + Radius;
00152 bbox->Extend(Tmp);
00153 }
00154
00155
00156 G4int G4FCylindricalSurface::Intersect( const G4Ray& ry )
00157 {
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 distance = kInfinity;
00171 closest_hit = PINFINITY;
00172
00173
00174 G4Point3D x = ry.GetStart();
00175 G4Vector3D dhat = ry.GetDir();
00176
00177
00178 G4Vector3D ahat = Position.GetAxis();
00179
00180
00181 G4double sol[2];
00182 sol[0]=-1.0;
00183 sol[1]=-1.0;
00184
00185
00186 G4Vector3D gamma = G4Vector3D( x - Position.GetLocation() );
00187
00188 G4double ga = gamma * ahat;
00189 G4double da = dhat * ahat;
00190
00191 G4double A = da * da - dhat * dhat;
00192 G4double B = 2 * ( -gamma * dhat + ga * da );
00193 G4double C = -gamma * gamma + ga * ga + radius * radius ;
00194
00195 G4double radical = B * B - 4.0 * A * C;
00196
00197 if ( radical < 0.0 )
00198
00199 return 0;
00200 else
00201 {
00202 G4double root = std::sqrt( radical );
00203 sol[0] = ( - B + root ) / ( 2. * A );
00204 sol[1] = ( - B - root ) / ( 2. * A );
00205 }
00206
00207
00208
00209 G4Point3D p0 = G4Point3D( x + sol[0]*dhat );
00210 G4Point3D p1 = G4Point3D( x + sol[1]*dhat );
00211
00212 if( !GetBBox()->Inside(p0) )
00213 sol[0] = kInfinity;
00214
00215 if( !GetBBox()->Inside(p1) )
00216 sol[1] = kInfinity;
00217
00218
00219
00220 G4int nbinter = 0;
00221 distance = kInfinity;
00222
00223 for ( G4int i = 0; i < 2; i++ )
00224 {
00225 if(sol[i] < kInfinity) {
00226 if ( sol[i] >= kCarTolerance*0.5 ) {
00227 nbinter ++;
00228
00229
00230 if( distance > sol[i]*sol[i]) {
00231 distance = sol[i]*sol[i];
00232 }
00233 }
00234 }
00235 }
00236
00237 return nbinter;
00238 }
00239
00240
00241 G4double G4FCylindricalSurface::HowNear( const G4Vector3D& x ) const
00242 {
00243
00244
00245
00246 G4double hownear;
00247
00248 G4Vector3D upcorner = G4Vector3D ( radius, 0 , origin.z()+length);
00249 G4Vector3D downcorner = G4Vector3D ( radius, 0 , origin.z());
00250 G4Vector3D xd;
00251
00252 xd = G4Vector3D ( std::sqrt ( x.x()*x.x() + x.y()*x.y() ) , 0 , x.z() );
00253
00254
00255 G4double Zinter = (xd.z()) ;
00256
00257 if ( ((Zinter >= downcorner.z()) && (Zinter <=upcorner.z())) ) {
00258 hownear = std::fabs( radius - xd.x() );
00259 } else {
00260 hownear = std::min ( (xd-upcorner).mag() , (xd-downcorner).mag() );
00261 }
00262
00263 return hownear;
00264 }
00265
00266 G4int G4FCylindricalSurface::WithinBoundary( const G4Vector3D& x ) const
00267 {
00268
00269
00270 if ( std::fabs( ( x - Position.GetLocation()) * Position.GetAxis() )
00271 <= 0.5 * length )
00272 return 1;
00273 else
00274 return 0;
00275 }
00276
00277
00278 G4double G4FCylindricalSurface::Scale() const
00279 {
00280
00281
00282
00283 if ( radius == 0.0 )
00284 return length;
00285 else
00286 return radius;
00287 }
00288
00289
00290 G4Vector3D G4FCylindricalSurface::SurfaceNormal( const G4Point3D& p ) const
00291 {
00292
00293
00294
00295 G4Vector3D n = G4Vector3D( ( p - Position.GetLocation() ) -
00296 ( ( p - Position.GetLocation()) *
00297 Position.GetAxis() ) * Position.GetAxis() );
00298 G4double nmag = n.mag();
00299
00300 if ( nmag != 0.0 )
00301 n = n * (1/nmag);
00302
00303 if( !sameSense )
00304 n = -n;
00305
00306 return n;
00307 }
00308
00309 G4int G4FCylindricalSurface::Inside ( const G4Vector3D& x ) const
00310 {
00311
00312
00313
00314
00315 if ( HowNear( x ) >= -0.5*kCarTolerance )
00316 return 1;
00317 else
00318 return 0;
00319 }
00320
00321
00322 void G4FCylindricalSurface::resize( G4double r, G4double l )
00323 {
00324
00325
00326 if ( r >= 0.0 )
00327 radius = r;
00328 else
00329 {
00330 std::ostringstream message;
00331 message << "Negative radius." << G4endl
00332 << "Original value of " << radius << " is retained.";
00333 G4Exception("G4FCylindricalSurface::resize()",
00334 "GeomSolids1001", JustWarning, message);
00335 }
00336
00337
00338 if ( l > 0.0 )
00339 length = l;
00340 else
00341 {
00342 std::ostringstream message;
00343 message << "Negative or zero length." << G4endl
00344 << "Original value of " << length << " is retained.";
00345 G4Exception("G4FCylindricalSurface::resize()",
00346 "GeomSolids1001", JustWarning, message);
00347 }
00348 }