Changeset 850 for trunk/source/geometry/solids/CSG/src
- Timestamp:
- Sep 10, 2008, 5:40:37 PM (16 years ago)
- Location:
- trunk/source/geometry/solids/CSG/src
- Files:
-
- 10 edited
-
G4Box.cc (modified) (1 diff)
-
G4CSGSolid.cc (modified) (1 diff)
-
G4Cons.cc (modified) (1 diff)
-
G4Orb.cc (modified) (1 diff)
-
G4Para.cc (modified) (1 diff)
-
G4Sphere.cc (modified) (17 diffs)
-
G4Torus.cc (modified) (1 diff)
-
G4Trap.cc (modified) (1 diff)
-
G4Trd.cc (modified) (1 diff)
-
G4Tubs.cc (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/geometry/solids/CSG/src/G4Box.cc
r831 r850 26 26 // 27 27 // $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // -
trunk/source/geometry/solids/CSG/src/G4CSGSolid.cc
r831 r850 26 26 // 27 27 // $Id: G4CSGSolid.cc,v 1.13 2006/10/19 15:33:37 gcosmo Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // -------------------------------------------------------------------- -
trunk/source/geometry/solids/CSG/src/G4Cons.cc
r831 r850 25 25 // 26 26 // 27 // $Id: G4Cons.cc,v 1.5 4.4.1 2008/04/23 09:05:23gcosmo Exp $28 // GEANT4 tag $Name: geant4-09-01-patch-02$27 // $Id: G4Cons.cc,v 1.56 2008/02/20 08:56:16 gcosmo Exp $ 28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // -
trunk/source/geometry/solids/CSG/src/G4Orb.cc
r831 r850 25 25 // 26 26 // $Id: G4Orb.cc,v 1.24 2007/05/18 07:38:01 gcosmo Exp $ 27 // GEANT4 tag $Name: $27 // GEANT4 tag $Name: HEAD $ 28 28 // 29 29 // class G4Orb -
trunk/source/geometry/solids/CSG/src/G4Para.cc
r831 r850 26 26 // 27 27 // $Id: G4Para.cc,v 1.39 2006/10/19 15:33:37 gcosmo Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // class G4Para -
trunk/source/geometry/solids/CSG/src/G4Sphere.cc
r831 r850 25 25 // 26 26 // 27 // $Id: G4Sphere.cc,v 1. 57.2.1 2008/04/23 09:05:23 gcosmoExp $28 // GEANT4 tag $Name: geant4-09-01-patch-02$27 // $Id: G4Sphere.cc,v 1.68 2008/07/07 09:35:16 grichine Exp $ 28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // class G4Sphere … … 34 34 // History: 35 35 // 36 // 12.06.08 V.Grichine: fix for theta intersections in DistanceToOut(p,v,...) 36 37 // 22.07.05 O.Link : Added check for intersection with double cone 37 38 // 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal … … 47 48 // 25.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), phi intersections 48 49 // 12.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), theta intersections 49 // 09.10.98 V.Grichine: modifications in Distance 50 // 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...) 50 51 // 17.09.96 V.Grichine: final modifications to commit 51 52 // 28.03.94 P.Kent: old C++ code converted to tolerant geometry … … 608 609 distETheta = std::fabs(pTheta-fSTheta-fDTheta); 609 610 610 nTs = G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi), 611 -std::cos(fSTheta)*std::sin(pPhi), 612 std::sin(fSTheta) ); 613 nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi), 614 std::cos(fSTheta+fDTheta)*std::sin(pPhi), 615 -std::sin(fSTheta+fDTheta) ); 611 nTs = G4ThreeVector(-std::cos(fSTheta)*p.x()/rho, // *std::cos(pPhi), 612 -std::cos(fSTheta)*p.y()/rho, // *std::sin(pPhi), 613 std::sin(fSTheta) ); 614 615 nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*p.x()/rho, // *std::cos(pPhi), 616 std::cos(fSTheta+fDTheta)*p.y()/rho, // *std::sin(pPhi), 617 -std::sin(fSTheta+fDTheta) ); 616 618 } 617 619 else if( !fRmin ) 618 620 { 619 if ( fSTheta ) distSTheta = 0.; 620 if ( fSTheta + fDTheta < pi ) distETheta = 0.; 621 if ( fSTheta ) 622 { 623 distSTheta = 0.; 624 nTs = G4ThreeVector(0.,0.,-1.); 625 } 626 if ( fSTheta + fDTheta < pi ) // distETheta = 0.; 627 { 628 distETheta = 0.; 629 nTe = G4ThreeVector(0.,0.,1.); 630 } 621 631 } 622 632 } … … 1423 1433 d = std::sqrt(d2) ; 1424 1434 s = -b - d ; // First root 1425 1426 if ( s < 0 ) 1427 { 1428 s=-b+d; // Second root 1435 zi = p.z() + s*v.z(); 1436 1437 if ( s < 0 || zi*(fSTheta - halfpi) > 0 ) 1438 { 1439 s = -b+d; // Second root 1429 1440 } 1430 1441 if (s >= 0 && s < snxt) 1431 1442 { 1432 xi = p.x() + s*v.x() 1433 yi = p.y() + s*v.y() 1434 zi = p.z() + s*v.z() 1435 rhoi2 = xi*xi + yi*yi 1436 radi2 = rhoi2 + zi*zi 1443 xi = p.x() + s*v.x(); 1444 yi = p.y() + s*v.y(); 1445 zi = p.z() + s*v.z(); 1446 rhoi2 = xi*xi + yi*yi; 1447 radi2 = rhoi2 + zi*zi; 1437 1448 if ( (radi2 <= tolORMax2) 1438 1449 && (radi2 >= tolORMin2) … … 1501 1512 } 1502 1513 } 1503 else if (pTheta > tolETheta) 1504 { // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0) 1505 // Inside (theta>etheta+tol) e theta cone 1514 else if ( pTheta > tolETheta ) 1515 { 1516 // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0) 1517 // Inside (theta > etheta+tol) e-theta cone 1506 1518 // First root of etheta cone, second if first root `imaginary' 1507 1519 … … 1517 1529 d = std::sqrt(d2) ; 1518 1530 s = -b - d ; // First root 1519 if (s < 0) 1531 zi = p.z() + s*v.z(); 1532 1533 if (s < 0 || zi*(fSTheta + fDTheta - halfpi) > 0) 1520 1534 { 1521 1535 s = -b + d ; // second root … … 1959 1973 1960 1974 G4bool segTheta; // Theta flag and precals 1961 G4double tanSTheta=0.,tanETheta , rhoSecTheta;1975 G4double tanSTheta=0.,tanETheta=0., rhoSecTheta; 1962 1976 G4double tanSTheta2=0.,tanETheta2=0.; 1963 G4double dist2STheta, dist2ETheta;1977 G4double dist2STheta, dist2ETheta, distTheta; 1964 1978 G4double d2,s; 1965 1979 1966 1980 // General Precalcs 1967 1981 1968 rho2 =p.x()*p.x()+p.y()*p.y();1969 rad2 =rho2+p.z()*p.z();1982 rho2 = p.x()*p.x()+p.y()*p.y(); 1983 rad2 = rho2+p.z()*p.z(); 1970 1984 // G4double rad=std::sqrt(rad2); 1971 1985 1972 pTheta =std::atan2(std::sqrt(rho2),p.z());1973 1974 pDotV2d =p.x()*v.x()+p.y()*v.y();1975 pDotV3d =pDotV2d+p.z()*v.z();1986 pTheta = std::atan2(std::sqrt(rho2),p.z()); 1987 1988 pDotV2d = p.x()*v.x()+p.y()*v.y(); 1989 pDotV3d = pDotV2d+p.z()*v.z(); 1976 1990 1977 1991 // Set phi divided flag and precalcs 1978 1992 1979 if( fDPhi<twopi)1993 if( fDPhi < twopi ) 1980 1994 { 1981 1995 segPhi=true; … … 1996 2010 // Theta precalcs 1997 2011 1998 if (fDTheta < pi) 1999 { 2000 segTheta=true; 2001 tolSTheta=fSTheta-kAngTolerance*0.5; 2002 tolETheta=fSTheta+fDTheta+kAngTolerance*0.5; 2003 } 2004 else 2005 { 2006 segTheta=false; 2007 } 2012 if ( fDTheta < pi ) 2013 { 2014 segTheta = true; 2015 tolSTheta = fSTheta - kAngTolerance*0.5; 2016 tolETheta = fSTheta + fDTheta + kAngTolerance*0.5; 2017 } 2018 else segTheta = false; 2019 2008 2020 2009 2021 // Radial Intersections from G4Sphere::DistanceToIn … … 2024 2036 // 2025 2037 // const G4double fractionTolerance = 1.0e-12; 2026 const G4double flexRadMaxTolerance = // kRadTolerance; 2038 2039 const G4double flexRadMaxTolerance = // kRadTolerance; 2027 2040 std::max(kRadTolerance, fEpsilon * fRmax); 2028 2041 2029 2042 const G4double Rmax_plus = fRmax + flexRadMaxTolerance*0.5; 2043 2030 2044 const G4double flexRadMinTolerance = std::max(kRadTolerance, 2031 2045 fEpsilon * fRmin); 2046 2032 2047 const G4double Rmin_minus= (fRmin > 0) ? fRmin-flexRadMinTolerance*0.5 : 0 ; 2033 2048 … … 2063 2078 else 2064 2079 { 2065 snxt =-pDotV3d+std::sqrt(d2); // second root since inside Rmax2066 side = kRMax ;2080 snxt = -pDotV3d+std::sqrt(d2); // second root since inside Rmax 2081 side = kRMax ; 2067 2082 } 2068 2083 } … … 2077 2092 d2 = pDotV3d*pDotV3d - c; 2078 2093 2079 if ( c >- flexRadMinTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin2094 if ( c >- flexRadMinTolerance*fRmin ) // 2.0 * (0.5*kRadTolerance) * fRmin 2080 2095 { 2081 2096 if( c < flexRadMinTolerance*fRmin && 2082 2097 d2 >= flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin 2083 2098 { 2084 if(calcNorm) 2085 { 2086 *validNorm = false ; // Rmin surface is concave 2087 } 2088 return snxt = 0 ; 2099 if(calcNorm) *validNorm = false ; // Rmin surface is concave 2100 return snxt = 0 ; 2089 2101 } 2090 2102 else 2091 2103 { 2092 if (d2 >= 0) 2093 { 2094 s = -pDotV3d-std::sqrt(d2) ; 2095 if (s>=0) // Always intersect Rmin first 2104 if ( d2 >= 0. ) 2105 { 2106 s = -pDotV3d-std::sqrt(d2); 2107 2108 if ( s >= 0. ) // Always intersect Rmin first 2096 2109 { 2097 2110 snxt = s ; … … 2129 2142 // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0 2130 2143 // 2144 2145 /* //////////////////////////////////////////////////////// 2146 2131 2147 tanSTheta=std::tan(fSTheta); 2132 2148 tanSTheta2=tanSTheta*tanSTheta; … … 2288 2304 s = kInfinity ; // wrong cone 2289 2305 } 2290 if (s < stheta)2291 {2292 stheta = s ;2293 sidetheta = kETheta;2294 }2306 } 2307 if (s < stheta) 2308 { 2309 stheta = s ; 2310 sidetheta = kETheta ; 2295 2311 } 2296 2312 } 2297 2313 } 2298 2314 } 2299 } 2315 */ //////////////////////////////////////////////////////////// 2316 2317 if(fSTheta) // intersection with first cons 2318 { 2319 2320 tanSTheta = std::tan(fSTheta); 2321 2322 if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0 2323 { 2324 if( v.z() > 0. ) 2325 { 2326 if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 ) 2327 { 2328 if(calcNorm) 2329 { 2330 *validNorm = true; 2331 *n = G4ThreeVector(0.,0.,1.); 2332 } 2333 return snxt = 0 ; 2334 } 2335 // s = -p.z()/v.z(); 2336 stheta = -p.z()/v.z(); 2337 sidetheta = kSTheta; 2338 } 2339 } 2340 else // kons is not plane 2341 { 2342 tanSTheta2 = tanSTheta*tanSTheta; 2343 t1 = 1-v.z()*v.z()*(1+tanSTheta2); 2344 t2 = pDotV2d-p.z()*v.z()*tanSTheta2; // ~vDotN if p on cons 2345 dist2STheta = rho2-p.z()*p.z()*tanSTheta2; // t3 2346 2347 // distTheta = std::sqrt(std::fabs(dist2STheta/(1+tanSTheta2))); 2348 distTheta = std::sqrt(rho2)-p.z()*tanSTheta; 2349 2350 if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons 2351 { 2352 if( v.z() > 0. ) 2353 { 2354 if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface 2355 { 2356 if( fSTheta < halfpi && p.z() > 0. ) 2357 { 2358 if( calcNorm ) *validNorm = false; 2359 return snxt = 0.; 2360 } 2361 else if( fSTheta > halfpi && p.z() <= 0) 2362 { 2363 if( calcNorm ) 2364 { 2365 *validNorm = true; 2366 if (rho2) 2367 { 2368 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2)); 2369 2370 *n = G4ThreeVector( p.x()/rhoSecTheta, 2371 p.y()/rhoSecTheta, 2372 std::sin(fSTheta) ); 2373 } 2374 else *n = G4ThreeVector(0.,0.,1.); 2375 } 2376 return snxt = 0.; 2377 } 2378 } 2379 // s = -0.5*dist2STheta/t2; 2380 2381 stheta = -0.5*dist2STheta/t2; 2382 sidetheta = kSTheta; 2383 } 2384 } 2385 else // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve 2386 { 2387 if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface 2388 { 2389 if( fSTheta > halfpi && t2 >= 0. ) // leave 2390 { 2391 if( calcNorm ) 2392 { 2393 *validNorm = true; 2394 if (rho2) 2395 { 2396 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2)); 2397 2398 *n = G4ThreeVector( p.x()/rhoSecTheta, 2399 p.y()/rhoSecTheta, 2400 std::sin(fSTheta) ); 2401 } 2402 else *n = G4ThreeVector(0.,0.,1.); 2403 } 2404 return snxt = 0.; 2405 } 2406 else if( fSTheta < halfpi && t2 < 0. && p.z() >=0. ) // leave 2407 { 2408 if( calcNorm ) *validNorm = false; 2409 return snxt = 0.; 2410 } 2411 } 2412 b = t2/t1; 2413 c = dist2STheta/t1; 2414 d2 = b*b - c ; 2415 2416 if ( d2 >= 0. ) 2417 { 2418 d = std::sqrt(d2); 2419 2420 if( fSTheta > halfpi ) 2421 { 2422 s = -b - d; // First root 2423 2424 if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) || 2425 s < 0. || 2426 ( s > 0. && p.z() + s*v.z() > 0.) ) 2427 { 2428 s = -b + d ; // 2nd root 2429 } 2430 if( s > flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.) 2431 { 2432 stheta = s; 2433 sidetheta = kSTheta; 2434 } 2435 } 2436 else // sTheta < pi/2, concave surface, no normal 2437 { 2438 s = -b - d; // First root 2439 2440 if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) || 2441 s < 0. || 2442 ( s > 0. && p.z() + s*v.z() < 0.) ) 2443 { 2444 s = -b + d ; // 2nd root 2445 } 2446 if( s > flexRadMaxTolerance*0.5 && p.z() + s*v.z() >= 0.) 2447 { 2448 stheta = s; 2449 sidetheta = kSTheta; 2450 } 2451 } 2452 } 2453 } 2454 } 2455 } 2456 if (fSTheta + fDTheta < pi) // intersection with second cons 2457 { 2458 2459 tanETheta = std::tan(fSTheta+fDTheta); 2460 2461 if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0 2462 { 2463 if( v.z() < 0. ) 2464 { 2465 if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 ) 2466 { 2467 if(calcNorm) 2468 { 2469 *validNorm = true; 2470 *n = G4ThreeVector(0.,0.,-1.); 2471 } 2472 return snxt = 0 ; 2473 } 2474 s = -p.z()/v.z(); 2475 2476 if( s < stheta) 2477 { 2478 stheta = s; 2479 sidetheta = kETheta; 2480 } 2481 } 2482 } 2483 else // kons is not plane 2484 { 2485 tanETheta2 = tanETheta*tanETheta; 2486 t1 = 1-v.z()*v.z()*(1+tanETheta2); 2487 t2 = pDotV2d-p.z()*v.z()*tanETheta2; // ~vDotN if p on cons 2488 dist2ETheta = rho2-p.z()*p.z()*tanETheta2; // t3 2489 2490 // distTheta = std::sqrt(std::fabs(dist2ETheta/(1+tanETheta2))); 2491 distTheta = std::sqrt(rho2)-p.z()*tanETheta; 2492 2493 if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons 2494 { 2495 if( v.z() < 0. ) 2496 { 2497 if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // p on surface 2498 { 2499 if( fSTheta+fDTheta > halfpi && p.z() < 0. ) 2500 { 2501 if( calcNorm ) *validNorm = false; 2502 return snxt = 0.; 2503 } 2504 else if( fSTheta+fDTheta < halfpi && p.z() >= 0) 2505 { 2506 if( calcNorm ) 2507 { 2508 *validNorm = true; 2509 if (rho2) 2510 { 2511 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2)); 2512 2513 *n = G4ThreeVector( p.x()/rhoSecTheta, 2514 p.y()/rhoSecTheta, 2515 -std::sin(fSTheta+fDTheta) ); 2516 } 2517 else *n = G4ThreeVector(0.,0.,-1.); 2518 } 2519 return snxt = 0.; 2520 } 2521 } 2522 s = -0.5*dist2ETheta/t2; 2523 2524 if( s < stheta) 2525 { 2526 stheta = s; 2527 sidetheta = kETheta; 2528 } 2529 } 2530 } 2531 else // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve 2532 { 2533 if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface 2534 { 2535 if( fSTheta+fDTheta < halfpi && t2 >= 0. ) // leave 2536 { 2537 if( calcNorm ) 2538 { 2539 *validNorm = true; 2540 if (rho2) 2541 { 2542 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2)); 2543 2544 *n = G4ThreeVector( p.x()/rhoSecTheta, 2545 p.y()/rhoSecTheta, 2546 -std::sin(fSTheta+fDTheta) ); 2547 } 2548 else *n = G4ThreeVector(0.,0.,-1.); 2549 } 2550 return snxt = 0.; 2551 } 2552 else if( fSTheta+fDTheta > halfpi && t2 < 0. && p.z() <=0. ) // leave 2553 { 2554 if( calcNorm ) *validNorm = false; 2555 return snxt = 0.; 2556 } 2557 } 2558 b = t2/t1; 2559 c = dist2ETheta/t1; 2560 d2 = b*b - c ; 2561 2562 if ( d2 >= 0. ) 2563 { 2564 d = std::sqrt(d2); 2565 2566 if( fSTheta+fDTheta < halfpi ) 2567 { 2568 s = -b - d; // First root 2569 2570 if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 < 0.) || 2571 s < 0. ) 2572 { 2573 s = -b + d ; // 2nd root 2574 } 2575 if( s > flexRadMaxTolerance*0.5 ) 2576 { 2577 if( s < stheta ) 2578 { 2579 stheta = s; 2580 sidetheta = kETheta; 2581 } 2582 } 2583 } 2584 else // sTheta+fDTheta > pi/2, concave surface, no normal 2585 { 2586 s = -b - d; // First root 2587 2588 if( (std::fabs(s) < flexRadMaxTolerance*0.5 && t2 >= 0.) || 2589 s < 0. || 2590 ( s > 0. && p.z() + s*v.z() > 0.) ) 2591 { 2592 s = -b + d ; // 2nd root 2593 } 2594 if( s > flexRadMaxTolerance*0.5 && p.z() + s*v.z() <= 0.) 2595 { 2596 if( s < stheta ) 2597 { 2598 stheta = s; 2599 sidetheta = kETheta; 2600 } 2601 } 2602 } 2603 } 2604 } 2605 } 2606 } 2607 2608 } // end theta intersections 2300 2609 2301 2610 // Phi Intersection … … 2578 2887 *validNorm=true; 2579 2888 break; 2889 2580 2890 case kRMin: 2581 2891 *validNorm=false; // Rmin is concave 2582 2892 break; 2893 2583 2894 case kSPhi: 2584 if ( fDPhi<=pi) // Normal to Phi-2895 if ( fDPhi <= pi ) // Normal to Phi- 2585 2896 { 2586 2897 *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); … … 2589 2900 else *validNorm=false; 2590 2901 break ; 2902 2591 2903 case kEPhi: 2592 if ( fDPhi<=pi) // Normal to Phi+2904 if ( fDPhi <= pi ) // Normal to Phi+ 2593 2905 { 2594 2906 *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); … … 2597 2909 else *validNorm=false; 2598 2910 break; 2911 2599 2912 case kSTheta: 2600 if( fSTheta == pi*0.5)2601 { 2602 *n=G4ThreeVector(0 ,0,1);2913 if( fSTheta == halfpi ) 2914 { 2915 *n=G4ThreeVector(0.,0.,1.); 2603 2916 *validNorm=true; 2604 2917 } 2605 else if ( fSTheta > pi ) 2918 else if ( fSTheta > halfpi ) 2919 { 2920 xi = p.x() + snxt*v.x(); 2921 yi = p.y() + snxt*v.y(); 2922 rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2)); 2923 *n = G4ThreeVector( xi/rhoSecTheta, // N- 2924 yi/rhoSecTheta, 2925 -tanSTheta/std::sqrt(1+tanSTheta2)); 2926 *validNorm=true; 2927 } 2928 else *validNorm=false; // Concave STheta cone 2929 break; 2930 2931 case kETheta: 2932 if( ( fSTheta + fDTheta ) == halfpi ) 2933 { 2934 *n = G4ThreeVector(0.,0.,-1.); 2935 *validNorm = true; 2936 } 2937 else if ( ( fSTheta + fDTheta ) < halfpi) 2606 2938 { 2607 2939 xi=p.x()+snxt*v.x(); 2608 2940 yi=p.y()+snxt*v.y(); 2609 rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2)) ; 2610 *n = G4ThreeVector(-xi/rhoSecTheta, // N- 2611 -yi/rhoSecTheta, 2612 tanSTheta/std::sqrt(1+tanSTheta2)) ; 2613 *validNorm=true; 2614 } 2615 else *validNorm=false; // Concave STheta cone 2616 break; 2617 case kETheta: 2618 if( ( fSTheta + fDTheta ) == pi*0.5 ) 2619 { 2620 *n = G4ThreeVector(0,0,-1); 2621 *validNorm = true ; 2622 } 2623 else if ( ( fSTheta + fDTheta ) < pi ) 2624 { 2625 xi=p.x()+snxt*v.x(); 2626 yi=p.y()+snxt*v.y(); 2627 rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2)) ; 2941 rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2)); 2628 2942 *n = G4ThreeVector( xi/rhoSecTheta, // N+ 2629 2943 yi/rhoSecTheta, 2630 -tan STheta/std::sqrt(1+tanSTheta2) );2944 -tanETheta/std::sqrt(1+tanETheta2) ); 2631 2945 *validNorm=true; 2632 2946 } 2633 2947 else *validNorm=false; // Concave ETheta cone 2634 2948 break; 2949 2635 2950 default: 2636 2951 G4cout.precision(16); -
trunk/source/geometry/solids/CSG/src/G4Torus.cc
r831 r850 26 26 // 27 27 // $Id: G4Torus.cc,v 1.63 2007/10/02 09:34:17 gcosmo Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // -
trunk/source/geometry/solids/CSG/src/G4Trap.cc
r831 r850 25 25 // 26 26 // 27 // $Id: G4Trap.cc,v 1.4 2.4.1 2008/04/23 09:55:26gcosmo Exp $28 // GEANT4 tag $Name: geant4-09-01-patch-02$27 // $Id: G4Trap.cc,v 1.45 2008/04/23 09:49:57 gcosmo Exp $ 28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // class G4Trap -
trunk/source/geometry/solids/CSG/src/G4Trd.cc
r831 r850 26 26 // 27 27 // $Id: G4Trd.cc,v 1.34 2006/10/19 15:33:38 gcosmo Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // -
trunk/source/geometry/solids/CSG/src/G4Tubs.cc
r831 r850 25 25 // 26 26 // 27 // $Id: G4Tubs.cc,v 1.6 6 2007/11/23 09:07:43 tnikitinExp $28 // GEANT4 tag $Name: $27 // $Id: G4Tubs.cc,v 1.68 2008/06/23 13:37:39 gcosmo Exp $ 28 // GEANT4 tag $Name: HEAD $ 29 29 // 30 30 // … … 983 983 iden = std::sqrt(t3) ; 984 984 cosPsi = inum/iden ; 985 if (cosPsi >= cosHDPhiIT) { return 0.0; } 985 if (cosPsi >= cosHDPhiIT) 986 { 987 // In the old version, the small negative tangent for the point 988 // on surface was not taken in account, and returning 0.0 ... 989 // New version: check the tangent for the point on surface and 990 // if no intersection, return kInfinity, if intersection instead 991 // return s. 992 // 993 c = t3-fRMax*fRMax; 994 if ( c<=0.0 ) 995 { 996 return 0.0; 997 } 998 else 999 { 1000 c = c/t1 ; 1001 d = b*b-c; 1002 if ( d>=0.0 ) 1003 { 1004 snxt = c/(-b+std::sqrt(d)); // using safe solution 1005 // for quadratic equation 1006 if ( snxt<kCarTolerance*0.5 ) { snxt=0; } 1007 return snxt ; 1008 } 1009 else 1010 { 1011 return kInfinity; 1012 } 1013 } 1014 } 986 1015 } 987 1016 else 988 { 989 return 0.0 ; 990 } 991 } 992 } 1017 { 1018 // In the old version, the small negative tangent for the point 1019 // on surface was not taken in account, and returning 0.0 ... 1020 // New version: check the tangent for the point on surface and 1021 // if no intersection, return kInfinity, if intersection instead 1022 // return s. 1023 // 1024 c = t3 - fRMax*fRMax; 1025 if ( c<=0.0 ) 1026 { 1027 return 0.0; 1028 } 1029 else 1030 { 1031 c = c/t1 ; 1032 d = b*b-c; 1033 if ( d>=0.0 ) 1034 { 1035 snxt= c/(-b+std::sqrt(d)); // using safe solution 1036 // for quadratic equation 1037 if ( snxt<kCarTolerance*0.5 ) { snxt=0; } 1038 return snxt ; 1039 } 1040 else 1041 { 1042 return kInfinity; 1043 } 1044 } 1045 } // end if (seg) 1046 } // end if (t3>tolIRMin2) 1047 } // end if (Inside Outer Radius) 993 1048 if ( fRMin ) // Try inner cylinder intersection 994 1049 { … … 1459 1514 // 1460 1515 if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)) 1461 1516 { sidephi = kSPhi; 1462 1517 if (((fSPhi-0.5*kAngTolerance)<=vphi) 1463 1518 &&((ePhi+0.5*kAngTolerance)>=vphi))
Note: See TracChangeset
for help on using the changeset viewer.
