00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // $Id$ 00028 // 00029 // ---------------------------------------------------------------------- 00030 // GEANT 4 class source file 00031 // 00032 // G4CylindricalSurface.cc 00033 // 00034 // ---------------------------------------------------------------------- 00035 00036 #include "G4CylindricalSurface.hh" 00037 #include "G4Sort.hh" 00038 #include "G4Globals.hh" 00039 00040 G4CylindricalSurface::G4CylindricalSurface() : G4Surface() 00041 { 00042 // default constructor 00043 // default axis is ( 1.0, 0.0, 0.0 ), default radius is 1.0 00044 00045 axis = G4Vector3D( 1.0, 0.0, 0.0 ); 00046 radius = 1.0; 00047 } 00048 00049 00050 G4CylindricalSurface::G4CylindricalSurface( const G4Vector3D& o, 00051 const G4Vector3D& a, 00052 G4double r ) 00053 : G4Surface() 00054 { 00055 // Normal constructor 00056 // require axis to be a unit vector 00057 00058 G4double amag = a.mag(); 00059 if ( amag != 0.0 ) 00060 { 00061 axis = a * (1/ amag); // this makes the axis a unit vector 00062 } 00063 else 00064 { 00065 std::ostringstream message; 00066 message << "Axis has zero length." << G4endl 00067 << "Default axis ( 1.0, 0.0, 0.0 ) is used."; 00068 G4Exception("G4CylindricalSurface::G4CylindricalSurface()", 00069 "GeomSolids1001", JustWarning, message); 00070 00071 axis = G4Vector3D( 1.0, 0.0, 0.0 ); 00072 } 00073 00074 // Require radius to be non-negative 00075 // 00076 if ( r >= 0.0 ) 00077 { 00078 radius = r; 00079 } 00080 else 00081 { 00082 std::ostringstream message; 00083 message << "Negative radius." << G4endl 00084 << "Default radius of 1.0 is used."; 00085 G4Exception("G4CylindricalSurface::G4CylindricalSurface()", 00086 "GeomSolids1001", JustWarning, message); 00087 00088 radius = 1.0; 00089 } 00090 00091 origin =o; 00092 } 00093 00094 00095 G4CylindricalSurface::~G4CylindricalSurface() 00096 { 00097 } 00098 00099 00100 G4CylindricalSurface::G4CylindricalSurface( const G4CylindricalSurface& c ) 00101 : G4Surface(), axis(c.axis), radius(c.radius) 00102 { 00103 } 00104 00105 00106 G4CylindricalSurface& 00107 G4CylindricalSurface::operator=( const G4CylindricalSurface& c ) 00108 { 00109 if (&c == this) { return *this; } 00110 axis = c.axis; 00111 radius = c.radius; 00112 00113 return *this; 00114 } 00115 00116 00117 const char* G4CylindricalSurface::NameOf() const 00118 { 00119 return "G4CylindricalSurface"; 00120 } 00121 00122 void G4CylindricalSurface::PrintOn( std::ostream& os ) const 00123 { 00124 os << "G4CylindricalSurface surface with origin: " << origin << "\t" 00125 << "radius: " << radius << "\tand axis " << axis << "\n"; 00126 } 00127 00128 G4int G4CylindricalSurface::Intersect(const G4Ray& ry) 00129 { 00130 // Distance along a Ray (straight line with G4ThreeVec) to leave or enter 00131 // a G4CylindricalSurface. The input variable which_way should be set 00132 // to +1 to indicate leaving a G4CylindricalSurface, -1 to indicate 00133 // entering a G4CylindricalSurface. 00134 // p is the point of intersection of the Ray with the G4CylindricalSurface. 00135 // If the G4Vector3D of the Ray is opposite to that of the Normal to 00136 // the G4CylindricalSurface at the intersection point, it will not leave 00137 // the G4CylindricalSurface. 00138 // Similarly, if the G4Vector3D of the Ray is along that of the Normal 00139 // to the G4CylindricalSurface at the intersection point, it will not enter 00140 // the G4CylindricalSurface. 00141 // This method is called by all finite shapes sub-classed to 00142 // G4CylindricalSurface. 00143 // Use the virtual function table to check if the intersection point 00144 // is within the boundary of the finite shape. 00145 // A negative result means no intersection. 00146 // If no valid intersection point is found, set the distance 00147 // and intersection point to large numbers. 00148 00149 G4int which_way=1; 00150 00151 if(!Inside(ry.GetStart())) { which_way = -1; } 00152 00153 distance = kInfinity; 00154 G4Vector3D lv ( kInfinity, kInfinity, kInfinity ); 00155 00156 closest_hit = lv; 00157 00158 // Origin and G4Vector3D unit vector of Ray. 00159 // 00160 G4Vector3D x = ry.GetStart(); 00161 G4Vector3D dhat = ry.GetDir(); 00162 00163 // Axis unit vector of the G4CylindricalSurface. 00164 // 00165 G4Vector3D ahat = GetAxis(); 00166 G4int isoln = 0, 00167 maxsoln = 2; 00168 00169 // Array of solutions in distance along the Ray. 00170 // 00171 G4double sol[2]; 00172 sol[0] = -1.0; 00173 sol[1] = -1.0 ; 00174 00175 // Calculate the two solutions (quadratic equation) 00176 // 00177 G4Vector3D d = x - GetOrigin(); 00178 G4double radiu = GetRadius(); 00179 00180 G4double dsq = d * d; 00181 G4double da = d * ahat; 00182 G4double dasq = da * da; 00183 G4double rsq = radiu * radiu; 00184 G4double qsq = dsq - dasq; 00185 G4double dira = dhat * ahat; 00186 G4double a = 1.0 - dira * dira; 00187 00188 if ( a <= 0.0 ) { return 0; } 00189 00190 G4double b = 2. * ( d * dhat - da * dira ); 00191 G4double c = rsq - qsq; 00192 G4double radical = b * b + 4. * a * c; 00193 00194 if ( radical < 0.0 ) { return 0; } 00195 00196 G4double root = std::sqrt( radical ); 00197 sol[0] = ( - b + root ) / ( 2. * a ); 00198 sol[1] = ( - b - root ) / ( 2. * a ); 00199 00200 // Order the possible solutions by increasing distance along the Ray 00201 // 00202 sort_double( sol, isoln, maxsoln-1 ); 00203 00204 // Now loop over each positive solution, keeping the first one (smallest 00205 // distance along the Ray) which is within the boundary of the sub-shape 00206 // and which also has the correct G4Vector3D with respect to the Normal to 00207 // the G4CylindricalSurface at the intersection point 00208 // 00209 for ( isoln = 0; isoln < maxsoln; isoln++ ) 00210 { 00211 if ( sol[isoln] >= kCarTolerance*0.5 ) 00212 { 00213 if ( sol[isoln] >= kInfinity ) // quit if too large 00214 { return 0; } 00215 00216 distance = sol[isoln]; 00217 closest_hit = ry.GetPoint( distance ); 00218 G4double tmp = dhat * (Normal( closest_hit )); 00219 00220 if ((tmp * which_way) >= 0.0 ) 00221 { 00222 if ( WithinBoundary( closest_hit ) == 1 ) 00223 { 00224 distance = distance*distance; 00225 } 00226 } 00227 return 1; 00228 } 00229 } 00230 00231 // Get here only if there was no solution within the boundary, Reset 00232 // distance and intersection point to large numbers 00233 // 00234 distance = kInfinity; 00235 closest_hit = lv; 00236 00237 return 0; 00238 } 00239 00240 00241 G4double G4CylindricalSurface::HowNear( const G4Vector3D& x ) const 00242 { 00243 // Distance from the point x to the infinite G4CylindricalSurface. 00244 // The distance will be positive if the point is Inside the 00245 // G4CylindricalSurface, negative if the point is outside. 00246 // Note that this may not be correct for a bounded cylindrical object 00247 // subclassed to G4CylindricalSurface. 00248 00249 G4Vector3D d = x - origin; 00250 G4double dA = d * axis; 00251 G4double rds = std::sqrt( d.mag2() - dA*dA ); 00252 G4double hownear = std::fabs( radius - rds ); 00253 00254 return hownear; 00255 } 00256 00257 00258 /* 00259 G4double G4CylindricalSurface:: 00260 distanceAlongRay( G4int which_way, const G4Ray* ry, G4Vector3D& p ) const 00261 { 00262 // Distance along a Ray (straight line with G4Vector3D) to leave or enter 00263 // a G4CylindricalSurface. The input variable which_way should be set to +1 00264 // to indicate leaving a G4CylindricalSurface, -1 to indicate entering a 00265 // G4CylindricalSurface. 00266 // p is the point of intersection of the Ray with the G4CylindricalSurface. 00267 // If the G4Vector3D of the Ray is opposite to that of the Normal to 00268 // the G4CylindricalSurface at the intersection point, it will not leave the 00269 // G4CylindricalSurface. 00270 // Similarly, if the G4Vector3D of the Ray is along that of the Normal 00271 // to the G4CylindricalSurface at the intersection point, it will not enter 00272 // the G4CylindricalSurface. 00273 // This method is called by all finite shapes sub-classed to 00274 // G4CylindricalSurface. 00275 // Use the virtual function table to check if the intersection point 00276 // is within the boundary of the finite shape. 00277 // A negative result means no intersection. 00278 // If no valid intersection point is found, set the distance 00279 // and intersection point to large numbers. 00280 00281 G4double Dist = kInfinity; 00282 G4Vector3D lv ( kInfinity, kInfinity, kInfinity ); 00283 p = lv; 00284 00285 // Origin and G4Vector3D unit vector of Ray. 00286 // 00287 G4Vector3D x = ry->Position(); 00288 G4Vector3D dhat = ry->Direction( 0.0 ); 00289 00290 // Axis unit vector of the G4CylindricalSurface. 00291 // 00292 G4Vector3D ahat = GetAxis(); 00293 G4int isoln = 0, maxsoln = 2; 00294 00295 // Array of solutions in distance along the Ray 00296 // 00297 G4double s[2];s[0] = -1.0; s[1]= -1.0 ; 00298 00299 // Calculate the two solutions (quadratic equation) 00300 // 00301 G4Vector3D d = x - GetOrigin(); 00302 G4double radius = GetRadius(); 00303 00304 // Quit with no intersection if radius of the G4CylindricalSurface is zero 00305 // 00306 if ( radius <= 0.0 ) { return Dist; } 00307 00308 G4double dsq = d * d; 00309 G4double da = d * ahat; 00310 G4double dasq = da * da; 00311 G4double rsq = radius * radius; 00312 G4double qsq = dsq - dasq; 00313 G4double dira = dhat * ahat; 00314 G4double a = 1.0 - dira * dira; 00315 if ( a <= 0.0 ) { return Dist; } 00316 00317 G4double b = 2. * ( d * dhat - da * dira ); 00318 G4double c = rsq - qsq; 00319 G4double radical = b * b + 4. * a * c; 00320 if ( radical < 0.0 ) { return Dist; } 00321 00322 G4double root = std::sqrt( radical ); 00323 s[0] = ( - b + root ) / ( 2. * a ); 00324 s[1] = ( - b - root ) / ( 2. * a ); 00325 00326 // Order the possible solutions by increasing distance along the Ray 00327 // 00328 G4Sort_double( s, isoln, maxsoln-1 ); 00329 00330 // Now loop over each positive solution, keeping the first one (smallest 00331 // distance along the Ray) which is within the boundary of the sub-shape 00332 // and which also has the correct G4Vector3D with respect to the Normal to 00333 // the G4CylindricalSurface at the intersection point 00334 // 00335 for ( isoln = 0; isoln < maxsoln; isoln++ ) 00336 { 00337 if ( s[isoln] >= 0.0 ) 00338 { 00339 if ( s[isoln] >= kInfinity ) // quit if too large 00340 { return Dist; } 00341 Dist = s[isoln]; 00342 p = ry->Position( Dist ); 00343 if ( ( ( dhat * Normal( p ) * which_way ) >= 0.0 ) 00344 && ( WithinBoundary( p ) == 1 ) ) { return Dist; } 00345 } 00346 } 00347 00348 // Get here only if there was no solution within the boundary, Reset 00349 // distance and intersection point to large numbers 00350 00351 p = lv; 00352 return kInfinity; 00353 } 00354 00355 00356 G4double G4CylindricalSurface:: 00357 distanceAlongHelix( G4int which_way, const Helix* hx, G4Vector3D& p ) const 00358 { 00359 // Distance along a Helix to leave or enter a G4CylindricalSurface. 00360 // The input variable which_way should be set to +1 to 00361 // indicate leaving a G4CylindricalSurface, -1 to indicate entering a 00362 // G4CylindricalSurface. 00363 // p is the point of intersection of the Helix with the G4CylindricalSurface. 00364 // If the G4Vector3D of the Helix is opposite to that of the Normal to 00365 // the G4CylindricalSurface at the intersection point, it will not leave the 00366 // G4CylindricalSurface. 00367 // Similarly, if the G4Vector3D of the Helix is along that of the Normal 00368 // to the G4CylindricalSurface at the intersection point, it will not enter 00369 // the G4CylindricalSurface. 00370 // This method is called by all finite shapes sub-classed to 00371 // G4CylindricalSurface. 00372 // Use the virtual function table to check if the intersection point 00373 // is within the boundary of the finite shape. 00374 // If no valid intersection point is found, set the distance 00375 // and intersection point to large numbers. 00376 // Possible negative distance solutions are discarded. 00377 00378 G4double Dist = kInfinity; 00379 G4Vector3D lv ( kInfinity, kInfinity, kInfinity ); 00380 G4Vector3D zerovec; // zero G4Vector3D 00381 00382 p = lv; 00383 G4int isoln = 0, maxsoln = 4; 00384 00385 // Array of solutions in turning angle 00386 // 00387 G4double s[4];s[0]=-1.0;s[1]= -1.0;s[2]= -1.0;s[3]= -1.0; 00388 00389 // Flag set to 1 if exact solution is found 00390 // 00391 G4int exact = 0; 00392 00393 // Helix parameters 00394 // 00395 G4double rh = hx->GetRadius(); // radius of Helix 00396 G4Vector3D ah = hx->GetAxis(); // axis of Helix 00397 G4Vector3D oh = hx->position(); // origin of Helix 00398 G4Vector3D dh = hx->direction( 0.0 ); // initial G4Vector3D of Helix 00399 G4Vector3D prp = hx->getPerp(); // perpendicular vector 00400 G4double prpmag = prp.Magnitude(); 00401 G4double rhp = rh / prpmag; 00402 00403 // G4CylindricalSurface parameters 00404 // 00405 G4double rc = GetRadius(); // radius of G4CylindricalSurface 00406 if ( rc == 0.0 ) { return Dist; } // quit if zero radius 00407 00408 G4Vector3D oc = GetOrigin(); // origin of G4CylindricalSurface 00409 G4Vector3D ac = GetAxis(); // axis of G4CylindricalSurface 00410 00411 // Calculate quantities of use later on 00412 // 00413 G4Vector3D alpha = rhp * prp; 00414 G4Vector3D beta = rhp * dh; 00415 G4Vector3D gamma = oh - oc; 00416 00417 // Declare variables used later on in several places 00418 // 00419 G4double rcd2 = 0.0, alpha2 = 0.0; 00420 G4double A = 0.0, B = 0.0, C = 0.0, F = 0.0, G = 0.0, H = 0.0; 00421 G4double CoverB = 0.0, radical = 0.0, root = 0.0, s1 = 0.0, s2 = 0.0; 00422 G4Vector3D ghat; 00423 00424 // Set flag for special cases 00425 // 00426 G4int special_case = 0; // 0 means general case 00427 00428 // Test to see if axes of Helix and G4CylindricalSurface are parallel, 00429 // in which case there are exact solutions 00430 // 00431 if ( ( std::fabs( ah.AngleBetween(ac) ) < kAngTolerance ) 00432 || ( std::fabs( ah.AngleBetween(ac) - pi ) < kAngTolerance ) ) 00433 { 00434 special_case = 1; 00435 00436 // If, in addition, gamma is a zero vector or is parallel to the 00437 // G4CylindricalSurface axis, this simplifies the previous case 00438 // 00439 if ( gamma == zerovec ) 00440 { 00441 special_case = 3; 00442 ghat = gamma; 00443 } 00444 else 00445 { 00446 ghat = gamma / gamma.Magnitude(); 00447 if ( ( std::fabs( ghat.AngleBetween(ac) ) < kAngTolerance ) 00448 || ( std::fabs( ghat.AngleBetween(ac) - pi ) < kAngTolerance ) ) 00449 { 00450 special_case = 3; 00451 } 00452 } 00453 00454 // Test to see if, in addition to the axes of the Helix and 00455 // G4CylindricalSurface being parallel, the axis of the surface is 00456 // perpendicular to the initial G4Vector3D of the Helix 00457 // 00458 if ( std::fabs( ( ac * dh ) ) < kAngTolerance ) 00459 { 00460 // And, if, in addition to all this, the difference in origins of the 00461 // Helix and G4CylindricalSurface is perpendicular to the initial 00462 // G4Vector3D of the Helix, there is a separate special case 00463 // 00464 if ( std::fabs( ( ghat * dh ) ) < kAngTolerance ) 00465 { 00466 special_case = 4; 00467 } 00468 } 00469 } // end of section with axes of Helix and G4CylindricalSurface parallel 00470 00471 // Another peculiar case occurs if the axis of the G4CylindricalSurface and 00472 // the initial G4Vector3D of the Helix line up and their origins are the same. 00473 // This will require a higher order approximation than the general case 00474 // 00475 if ( ( ( std::fabs( dh.AngleBetween(ac) ) < kAngTolerance ) 00476 || ( std::fabs( dh.AngleBetween(ac) - pi ) < kAngTolerance ) ) 00477 && ( gamma == zerovec ) ) 00478 { 00479 special_case = 2; 00480 } 00481 00482 // Now all the special cases have been tagged, so solutions are found 00483 // for each case. Exact solutions are indicated by setting exact = 1. 00484 // [For some reason switch doesn't work here, so do series of if's.] 00485 00486 if ( special_case == 0 ) // approximate quadratic solutions 00487 { 00488 A = beta * beta - ( beta * ac ) * ( beta * ac ) 00489 + gamma * alpha - ( gamma * ac ) * ( alpha * ac ); 00490 B = 2.0 * gamma * beta 00491 - 2.0 * ( gamma * ac ) * ( beta * ac ); 00492 C = gamma * gamma 00493 - ( gamma * ac ) * ( gamma * ac ) - rc * rc; 00494 00495 if ( std::fabs( A ) < kCarTolerance ) // no quadratic term 00496 { 00497 if ( B == 0.0 ) // no intersection, quit 00498 { 00499 return Dist; 00500 } 00501 else // B != 0 00502 { 00503 s[0] = -C / B; 00504 } 00505 } 00506 else // A != 0, general quadratic solution 00507 { 00508 radical = B * B - 4.0 * A * C; 00509 if ( radical < 0.0 ) // no solution, quit 00510 { 00511 return Dist; 00512 } 00513 root = std::sqrt( radical ); 00514 s[0] = ( -B + root ) / ( 2.0 * A ); 00515 s[1] = ( -B - root ) / ( 2.0 * A ); 00516 if ( rh < 0.0 ) 00517 { 00518 s[0] = -s[0]; 00519 s[1] = -s[1]; 00520 } 00521 s[2] = s[0] + 2.0 * pi; 00522 s[3] = s[1] + 2.0 * pi; 00523 } 00524 } 00525 else if ( special_case == 1 ) // exact solutions 00526 { 00527 exact = 1; 00528 H = 2.0 * ( alpha * alpha + gamma * alpha ); 00529 F = gamma * gamma - ( ( gamma * ac ) * ( gamma * ac ) ) - rc * rc + H; 00530 G = 2.0 * rhp * ( gamma * dh - ( gamma * ac ) * ( ac * dh ) ); 00531 A = G * G + H * H; 00532 B = -2.0 * F * H; 00533 C = F * F - G * G; 00534 00535 if ( std::fabs( A ) < kCarTolerance ) // no quadratic term 00536 { 00537 if ( B == 0.0 ) // no intersection, quit 00538 { 00539 return Dist; 00540 } 00541 else // B != 0 00542 { 00543 CoverB = -C / B; 00544 if ( std::fabs( CoverB ) > 1.0 ) 00545 { 00546 return Dist; 00547 } 00548 s[0] = std::acos( CoverB ); 00549 } 00550 } 00551 } 00552 else // A != 0, general quadratic solution 00553 { 00554 // Try a different method of calculation using F, G, and H to avoid 00555 // precision problems. 00556 00557 // radical = B * B - 4.0 * A * C; 00558 // if ( radical < 0.0 ) 00559 // { 00560 if ( std::fabs( H ) > kCarTolerance ) 00561 { 00562 G4double r1 = G / H; 00563 G4double r2 = F / H; 00564 G4double radsq = 1.0 + r1*r1 - r2*r2; 00565 if ( radsq < 0.0 ) { return Dist; } 00566 root = G * std::sqrt( radsq ); 00567 00568 G4double denominator = H * ( 1.0 + r1*r1 ); 00569 s1 = ( F + root ) / denominator; 00570 s2 = ( F - root ) / denominator; 00571 } 00572 else 00573 { 00574 return Dist; 00575 } 00576 // } // end radical < 0 condition 00577 // else 00578 // { 00579 // root = std::sqrt( radical ); 00580 // s1 = ( -B + root ) / ( 2.0 * A ); 00581 // s2 = ( -B - root ) / ( 2.0 * A ); 00582 // } 00583 if ( std::fabs( s1 ) <= 1.0 ) 00584 { 00585 s[0] = std::acos( s1 ); 00586 s[2] = 2.0 * pi - s[0]; 00587 } 00588 if ( std::fabs( s2 ) <= 1.0 ) 00589 { 00590 s[1] = std::acos( s2 ); 00591 s[3] = 2.0 * pi - s[1]; 00592 } 00593 00594 // Must take only solutions which satisfy original unsquared equation: 00595 // Gsin(s) - Hcos(s) + F = 0. Take best solution of pair and set false 00596 // solutions to -1. Only do this if the result is significantly different 00597 // from zero. 00598 00599 G4double temp1 = 0.0, temp2 = 0.0; 00600 G4double rsign = 1.0; 00601 if ( rh < 0.0 ) rsign = -1.0; 00602 if ( s[0] > 0.0 ) 00603 { 00604 temp1 = G * rsign * std::sin( s[0] ) - H * std::cos( s[0] ) + F; 00605 temp2 = G * rsign * std::sin( s[2] ) - H * std::cos( s[2] ) + F; 00606 if ( std::fabs( temp1 ) > std::fabs( temp2 ) ) 00607 { 00608 if ( std::fabs( temp1 ) > kAngTolerance ) { s[0] = -1.0; } 00609 else if ( std::fabs( temp2 ) > kAngTolerance ) { s[2] = -1.0; } 00610 } 00611 } 00612 if ( s[1] > 0.0 ) 00613 { 00614 temp1 = G * rsign * std::sin( s[1] ) - H * std::cos( s[1] ) + F; 00615 temp2 = G * rsign * std::sin( s[3] ) - H * std::cos( s[3] ) + F; 00616 if ( std::fabs( temp1 ) > std::fabs( temp2 ) ) 00617 { 00618 if ( std::fabs( temp1 ) > kAngTolerance ) { s[1] = -1.0; } 00619 else if ( std::fabs( temp2 ) > kAngTolerance ) { s[3] = -1.0; } 00620 } 00621 } 00622 } 00623 else if ( special_case == 2 ) // approximate solution 00624 { 00625 G4Vector3D e = ah.cross( ac ); 00626 G4double re = std::fabs( rhp ) * e.Magnitude(); 00627 s[0] = std::sqrt( 2.0 * rc / re ); 00628 } 00629 else if ( special_case == 3 ) // exact solutions 00630 { 00631 exact = 1; 00632 alpha2 = alpha * alpha; 00633 rcd2 = rhp * rhp * ( 1.0 - ( (ac*dh) * (ac*dh) ) ); 00634 A = alpha2 - rcd2; 00635 B = - 2.0 * alpha2; 00636 C = alpha2 + rcd2 - rc*rc; 00637 if ( std::fabs( A ) < kCarTolerance ) // no quadratic term 00638 { 00639 if ( B == 0.0 ) { return Dist; } // no intersection, quit 00640 else 00641 { 00642 CoverB = -C / B; 00643 if ( std::fabs( CoverB ) > 1.0 ) { return Dist; } 00644 s[0] = std::acos( CoverB ); 00645 } 00646 } 00647 else // A != 0, general quadratic solution 00648 { 00649 radical = B * B - 4.0 * A * C; 00650 if ( radical < 0.0 ) { return Dist; } 00651 root = std::sqrt( radical ); 00652 s1 = ( -B + root ) / ( 2.0 * A ); 00653 s2 = ( -B - root ) / ( 2.0 * A ); 00654 if ( std::fabs( s1 ) <= 1.0 ) { s[0] = std::acos( s1 ); } 00655 if ( std::fabs( s2 ) <= 1.0 ) { s[1] = std::acos( s2 ); } 00656 } 00657 } 00658 else if ( special_case == 4 ) // exact solution 00659 { 00660 exact = 1; 00661 F = gamma * gamma - ( ( gamma * ac ) * ( gamma * ac ) ) - rc * rc; 00662 G = 2.0 * ( rhp * rhp + gamma * alpha ); 00663 if ( G == 0.0 ) { return Dist; } // no intersection, quit 00664 00665 G4double cs = 1.0 + ( F / G ); 00666 if ( std::fabs( cs ) > 1.0 ) { return Dist; } // no intersection, quit 00667 s[0] = std::acos( cs ); 00668 } 00669 else // shouldn't get here 00670 { 00671 return Dist; 00672 } 00673 00674 // ************************************************************************** 00675 // 00676 // Order the possible solutions by increasing turning angle 00677 // 00678 G4Sort_double( s, isoln, maxsoln-1 ); 00679 00680 // Now loop over each positive solution, keeping the first one (smallest 00681 // distance along the Helix) which is within the boundary of the sub-shape 00682 // 00683 for ( isoln = 0; isoln < maxsoln; isoln++ ) 00684 { 00685 if ( s[isoln] >= 0.0 ) 00686 { 00687 // Calculate distance along Helix and position and G4Vector3D vectors 00688 // 00689 Dist = s[isoln] * std::fabs( rhp ); 00690 p = hx->position( Dist ); 00691 G4Vector3D d = hx->direction( Dist ); 00692 if ( exact == 0 ) // only for approximate solutions 00693 { 00694 // Now do approximation to get remaining distance to correct this 00695 // solution iterate it until the accuracy is below the user-set 00696 // surface precision. 00697 00698 G4double delta = 0.0; 00699 G4double delta0 = kInfinity; 00700 G4int dummy = 1; 00701 G4int iter = 0; 00702 G4int in0 = Inside( hx->position ( 0.0 ) ); 00703 G4int in1 = Inside( p ); 00704 G4double sc = Scale(); 00705 while ( dummy ) 00706 { 00707 iter++; 00708 00709 // Terminate loop after 50 iterations and Reset distance to large 00710 // number, indicating no intersection with G4CylindricalSurface. 00711 // This generally occurs if Helix curls too tightly to Intersect it. 00712 // 00713 if ( iter > 50 ) 00714 { 00715 Dist = kInfinity; 00716 p = lv; 00717 break; 00718 } 00719 00720 // Find distance from the current point along the above-calculated 00721 // G4Vector3D using a Ray. The G4Vector3D of the Ray and the Sign of 00722 // the distance are determined by whether the starting point of the 00723 // Helix is Inside or outside of the G4CylindricalSurface. 00724 // 00725 in1 = Inside( p ); 00726 if ( in1 ) // current point Inside 00727 { 00728 if ( in0 ) // starting point Inside 00729 { 00730 Ray* r = new Ray( p, d ); 00731 delta = distanceAlongRay( 1, r, p ); 00732 delete r; 00733 } 00734 else // starting point outside 00735 { 00736 Ray* r = new Ray( p, -d ); 00737 delta = -distanceAlongRay( 1, r, p ); 00738 delete r; 00739 } 00740 } 00741 else // current point outside 00742 { 00743 if ( in0 ) // starting point Inside 00744 { 00745 Ray* r = new Ray( p, -d ); 00746 delta = -distanceAlongRay( -1, r, p ); 00747 delete r; 00748 } 00749 else // starting point outside 00750 { 00751 Ray* r = new Ray( p, d ); 00752 delta = distanceAlongRay( -1, r, p ); 00753 delete r; 00754 } 00755 } 00756 00757 // Test if distance is less than the surface precision, 00758 // if so Terminate loop. 00759 // 00760 if ( std::fabs( delta / sc ) <= SURFACE_PRECISION ) 00761 { break; } 00762 00763 // If delta has not changed sufficiently from the previous 00764 // iteration, skip out of this loop. 00765 // 00766 if ( std::fabs( ( delta - delta0 ) / sc ) <= SURFACE_PRECISION ) 00767 { break; } 00768 00769 // If delta has increased in absolute value from the previous 00770 // iteration either the Helix doesn't Intersect the surface or the 00771 // approximate solution is too far from the real solution. 00772 // Try groping for a solution. If not found, Reset distance to large 00773 // number, indicating no intersection with the G4CylindricalSurface. 00774 00775 if ( std::fabs( delta ) > std::fabs( delta0 ) ) 00776 { 00777 Dist = std::fabs( rhp ) * gropeAlongHelix( hx ); 00778 if ( Dist < 0.0 ) 00779 { 00780 Dist = kInfinity; 00781 p = lv; 00782 } 00783 else 00784 { 00785 p = hx->position( Dist ); 00786 } 00787 break; 00788 } 00789 00790 // Set old delta to new one 00791 // 00792 delta0 = delta; 00793 00794 // Add distance to G4CylindricalSurface to distance along Helix 00795 // 00796 Dist += delta; 00797 00798 // Negative distance along Helix means Helix doesn't Intersect 00799 // G4CylindricalSurface. 00800 // Reset distance to large number, indicating no intersection with 00801 // G4CylindricalSurface. 00802 // 00803 if ( Dist < 0.0 ) 00804 { 00805 Dist = kInfinity; 00806 p = lv; 00807 break; 00808 } 00809 00810 // Recalculate point along Helix and the G4Vector3D 00811 // 00812 p = hx->position( Dist ); 00813 d = hx->direction( Dist ); 00814 } // end of while loop 00815 } // end of exact == 0 condition 00816 00817 // Now have best value of distance along Helix and position for this 00818 // solution, so test if it is within the boundary of the sub-shape 00819 // and require that it point in the correct G4Vector3D with respect to 00820 // the Normal to the G4CylindricalSurface. 00821 00822 if ( ( Dist < kInfinity ) && 00823 ( ( hx->direction( Dist ) * Normal( p ) * which_way ) >= 0.0 ) && 00824 ( WithinBoundary( p ) == 1 ) ) { return Dist; } 00825 00826 } // end of if s[isoln] >= 0.0 condition 00827 } // end of for loop over solutions 00828 00829 // If one gets here, there is no solution, so set distance along Helix 00830 // and position to large numbers 00831 00832 Dist = kInfinity; 00833 p = lv; 00834 00835 return Dist; 00836 } 00837 */ 00838 00839 00840 G4Vector3D G4CylindricalSurface::Normal( const G4Vector3D& p ) const 00841 { 00842 // return the Normal unit vector to the G4CylindricalSurface 00843 // at a point p on (or nearly on) the G4CylindricalSurface 00844 00845 G4Vector3D n = ( p - origin ) - ( ( p - origin ) * axis ) * axis; 00846 G4double nmag = n.mag(); 00847 00848 if ( nmag != 0.0 ) { n = n * (1/nmag); } 00849 00850 return n; 00851 } 00852 00853 00854 G4Vector3D G4CylindricalSurface::SurfaceNormal( const G4Point3D& p ) const 00855 { 00856 // return the Normal unit vector to the G4CylindricalSurface at a point 00857 // p on (or nearly on) the G4CylindricalSurface 00858 00859 G4Vector3D n = ( p - origin ) - ( ( p - origin ) * axis ) * axis; 00860 G4double nmag = n.mag(); 00861 00862 if ( nmag != 0.0 ) { n = n * (1/nmag); } 00863 00864 return n; 00865 } 00866 00867 00868 G4int G4CylindricalSurface::Inside ( const G4Vector3D& x ) const 00869 { 00870 // Return 0 if point x is outside G4CylindricalSurface, 1 if Inside. 00871 // Outside means that the distance to the G4CylindricalSurface would 00872 // be negative. Use the HowNear function to calculate this distance. 00873 00874 if ( HowNear( x ) >= -0.5*kCarTolerance ) 00875 { return 1; } 00876 else 00877 { return 0; } 00878 } 00879 00880 00881 G4int G4CylindricalSurface::WithinBoundary( const G4Vector3D& x ) const 00882 { 00883 // return 1 if point x is on the G4CylindricalSurface, otherwise return zero 00884 // base this on the surface precision factor 00885 00886 if ( std::fabs( HowNear( x ) / Scale() ) <= SURFACE_PRECISION ) 00887 { return 1; } 00888 else 00889 { return 0; } 00890 } 00891 00892 00893 G4double G4CylindricalSurface::Scale() const 00894 { 00895 // Returns the radius of a G4CylindricalSurface unless it is zero, in which 00896 // case returns the arbitrary number 1.0. 00897 // This is ok since derived finite-sized classes will overwrite this. 00898 // Used for Scale-invariant tests of surface thickness. 00899 00900 if ( radius == 0.0 ) 00901 { return 1.0; } 00902 else 00903 { return radius; } 00904 } 00905 00906 /* 00907 void G4CylindricalSurface:: 00908 rotate( G4double alpha, G4double beta, G4double gamma, 00909 G4ThreeMat& m, G4int inverse ) 00910 { 00911 // Rotate G4CylindricalSurface first about global x-axis by angle alpha, 00912 // second about global y-axis by angle beta, 00913 // and third about global z-axis by angle gamma 00914 // by creating and using G4ThreeMat objects in Surface::rotate 00915 // angles are assumed to be given in radians 00916 // if inverse is non-zero, the order of rotations is reversed 00917 // the axis is rotated here, the origin is rotated by calling 00918 // Surface::rotate 00919 00920 G4Surface::rotate( alpha, beta, gamma, m, inverse ); 00921 axis = m * axis; 00922 } 00923 00924 void G4CylindricalSurface:: 00925 rotate( G4double alpha, G4double beta, G4double gamma, G4int inverse ) 00926 { 00927 // Rotate G4CylindricalSurface first about global x-axis by angle alpha, 00928 // second about global y-axis by angle beta, 00929 // and third about global z-axis by angle gamma 00930 // by creating and using G4ThreeMat objects in Surface::rotate 00931 // angles are assumed to be given in radians 00932 // if inverse is non-zero, the order of rotations is reversed 00933 // the axis is rotated here, the origin is rotated by calling 00934 // Surface::rotate 00935 00936 G4ThreeMat m; 00937 G4Surface::rotate( alpha, beta, gamma, m, inverse ); 00938 axis = m * axis; 00939 } 00940 */ 00941 00942 void G4CylindricalSurface::SetRadius( G4double r ) 00943 { 00944 // Reset the radius of the G4CylindricalSurface 00945 // Require radius to be non-negative 00946 00947 if ( r >= 0.0 ) { radius = r; } 00948 else // Use old value (do not change radius) if out of the range, 00949 { // but print warning message 00950 std::ostringstream message; 00951 message << "Negative radius." << G4endl 00952 << "Default radius of " << radius << " is used."; 00953 G4Exception("G4CylindricalSurface::SetRadius()", 00954 "GeomSolids1001", JustWarning, message); 00955 } 00956 } 00957 00958 /* 00959 G4double G4CylindricalSurface::gropeAlongHelix( const Helix* hx ) const 00960 { 00961 // Grope for a solution of a Helix intersecting a G4CylindricalSurface. 00962 // This function returns the turning angle (in radians) where the 00963 // intersection occurs with only positive values allowed, or -1.0 if 00964 // no intersection is found. 00965 // The idea is to start at the beginning of the Helix, then take steps 00966 // of some fraction of a turn. If at the end of a Step, the current position 00967 // along the Helix and the previous position are on opposite sides of the 00968 // G4CylindricalSurface, then the solution must lie somewhere in between. 00969 00970 G4int one_over_f = 8; // one over fraction of a turn to go in each Step 00971 G4double turn_angle = 0.0; 00972 G4double dist_along = 0.0; 00973 G4double d_new; 00974 G4double fk = 1.0 / G4double( one_over_f ); 00975 G4double scal = Scale(); 00976 G4double d_old = HowNear( hx->position( dist_along ) ); 00977 G4double rh = hx->GetRadius(); // radius of Helix 00978 G4Vector3D prp = hx->getPerp(); // perpendicular vector 00979 G4double prpmag = prp.Magnitude(); 00980 G4double rhp = rh / prpmag; 00981 G4int max_iter = one_over_f * HELIX_MAX_TURNS; 00982 00983 // Take up to a user-settable number of turns along the Helix, 00984 // groping for an intersection point. 00985 // 00986 for ( G4int k = 1; k < max_iter; k++ ) 00987 { 00988 turn_angle = 2.0 * pi * k / one_over_f; 00989 dist_along = turn_angle * std::fabs( rhp ); 00990 d_new = HowNear( hx->position( dist_along ) ); 00991 if ( ( d_old < 0.0 && d_new > 0.0 ) || 00992 ( d_old > 0.0 && d_new < 0.0 ) ) 00993 { 00994 d_old = d_new; 00995 00996 // Old and new points are on opposite sides of the G4CylindricalSurface, 00997 // therefore a solution lies in between, use a binary search to pin the 00998 // point down to the surface precision, but don't do more than 50 00999 // iterations. 01000 01001 G4int itr = 0; 01002 while ( std::fabs( d_new / scal ) > SURFACE_PRECISION ) 01003 { 01004 itr++; 01005 if ( itr > 50 ) { return turn_angle; } 01006 turn_angle -= fk * pi; 01007 dist_along = turn_angle * std::fabs( rhp ); 01008 d_new = HowNear( hx->position( dist_along ) ); 01009 if ( ( d_old < 0.0 && d_new > 0.0 ) || ( d_old > 0.0 && d_new < 0.0 ) ) 01010 { fk *= -0.5; } 01011 else 01012 { fk *= 0.5; } 01013 d_old = d_new; 01014 } // end of while loop 01015 return turn_angle; // this is the best solution 01016 } // end of if condition 01017 } // end of for loop 01018 01019 // Get here only if no solution is found, so return -1.0 to indicate that. 01020 // 01021 return -1.0; 01022 } 01023 */