Changeset 1315 for trunk/source/geometry/solids/CSG/src
- Timestamp:
- Jun 18, 2010, 11:42:07 AM (14 years ago)
- Location:
- trunk/source/geometry/solids/CSG/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/geometry/solids/CSG/src/G4Box.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4Box.cc,v 1.4 4 2006/10/19 15:33:37gcosmo Exp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4Box.cc,v 1.49 2010/05/25 10:14:41 gcosmo Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // … … 32 32 // Implementation for G4Box class 33 33 // 34 // 24.06.98 - V. 34 // 24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v) 35 35 // 20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v) 36 36 // 07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0 … … 67 67 if ( (pX > 2*kCarTolerance) 68 68 && (pY > 2*kCarTolerance) 69 && (pZ > 2*kCarTolerance) ) 69 && (pZ > 2*kCarTolerance) ) // limit to thickness of surfaces 70 70 { 71 71 fDx = pX ; … … 105 105 void G4Box::SetXHalfLength(G4double dx) 106 106 { 107 if(dx > 2*kCarTolerance) 107 if(dx > 2*kCarTolerance) // limit to thickness of surfaces 108 { 108 109 fDx = dx; 110 } 109 111 else 110 112 { … … 122 124 void G4Box::SetYHalfLength(G4double dy) 123 125 { 124 if(dy > 2*kCarTolerance) 126 if(dy > 2*kCarTolerance) // limit to thickness of surfaces 127 { 125 128 fDy = dy; 129 } 126 130 else 127 131 { … … 139 143 void G4Box::SetZHalfLength(G4double dz) 140 144 { 141 if(dz > 2*kCarTolerance) 145 if(dz > 2*kCarTolerance) // limit to thickness of surfaces 146 { 142 147 fDz = dz; 148 } 143 149 else 144 150 { … … 153 159 fpPolyhedron = 0; 154 160 } 155 156 157 161 158 162 //////////////////////////////////////////////////////////////////////// … … 193 197 if (pVoxelLimit.IsXLimited()) 194 198 { 195 if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance||196 xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false ;199 if ((xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) || 200 (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance)) { return false ; } 197 201 else 198 202 { 199 if (xMin < pVoxelLimit.GetMinXExtent()) 200 { 201 xMin = pVoxelLimit.GetMinXExtent() ; 202 } 203 if (xMax > pVoxelLimit.GetMaxXExtent()) 204 { 205 xMax = pVoxelLimit.GetMaxXExtent() ; 206 } 203 xMin = std::max(xMin, pVoxelLimit.GetMinXExtent()); 204 xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent()); 207 205 } 208 206 } … … 213 211 if (pVoxelLimit.IsYLimited()) 214 212 { 215 if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance||216 yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false ;213 if ((yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) || 214 (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance)) { return false ; } 217 215 else 218 216 { 219 if (yMin < pVoxelLimit.GetMinYExtent()) 220 { 221 yMin = pVoxelLimit.GetMinYExtent() ; 222 } 223 if (yMax > pVoxelLimit.GetMaxYExtent()) 224 { 225 yMax = pVoxelLimit.GetMaxYExtent() ; 226 } 217 yMin = std::max(yMin, pVoxelLimit.GetMinYExtent()); 218 yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent()); 227 219 } 228 220 } … … 233 225 if (pVoxelLimit.IsZLimited()) 234 226 { 235 if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance||236 zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false ;227 if ((zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) || 228 (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance)) { return false ; } 237 229 else 238 230 { 239 if (zMin < pVoxelLimit.GetMinZExtent()) 240 { 241 zMin = pVoxelLimit.GetMinZExtent() ; 242 } 243 if (zMax > pVoxelLimit.GetMaxZExtent()) 244 { 245 zMax = pVoxelLimit.GetMaxZExtent() ; 246 } 231 zMin = std::max(zMin, pVoxelLimit.GetMinZExtent()); 232 zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent()); 247 233 } 248 234 } … … 286 272 if (pVoxelLimit.IsLimited(pAxis) == false) 287 273 { 288 if ( pMin != kInfinity || pMax != -kInfinity)274 if ( (pMin != kInfinity) || (pMax != -kInfinity) ) 289 275 { 290 276 existsAfterClip = true ; … … 292 278 // Add 2*tolerance to avoid precision troubles 293 279 294 pMin 295 pMax 280 pMin -= kCarTolerance; 281 pMax += kCarTolerance; 296 282 } 297 283 } … … 303 289 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); 304 290 305 if ( pMin != kInfinity || pMax != -kInfinity)291 if ( (pMin != kInfinity) || (pMax != -kInfinity) ) 306 292 { 307 293 existsAfterClip = true ; … … 356 342 EInside G4Box::Inside(const G4ThreeVector& p) const 357 343 { 344 static const G4double delta=0.5*kCarTolerance; 358 345 EInside in = kOutside ; 359 360 if ( std::fabs(p.x()) <= fDx - kCarTolerance*0.5 ) 361 { 362 if (std::fabs(p.y()) <= fDy - kCarTolerance*0.5 ) 363 { 364 if (std::fabs(p.z()) <= fDz - kCarTolerance*0.5 ) in = kInside ; 365 else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ; 366 } 367 else if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 ) 368 { 369 if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ; 370 } 371 } 372 else if (std::fabs(p.x()) <= fDx + kCarTolerance*0.5 ) 373 { 374 if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 ) 375 { 376 if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5) in = kSurface ; 346 G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z())); 347 348 if ( q.x() <= (fDx - delta) ) 349 { 350 if (q.y() <= (fDy - delta) ) 351 { 352 if ( q.z() <= (fDz - delta) ) { in = kInside ; } 353 else if ( q.z() <= (fDz + delta) ) { in = kSurface ; } 354 } 355 else if ( q.y() <= (fDy + delta) ) 356 { 357 if ( q.z() <= (fDz + delta) ) { in = kSurface ; } 358 } 359 } 360 else if ( q.x() <= (fDx + delta) ) 361 { 362 if ( q.y() <= (fDy + delta) ) 363 { 364 if ( q.z() <= (fDz + delta) ) { in = kSurface ; } 377 365 } 378 366 } … … 389 377 { 390 378 G4double distx, disty, distz ; 391 G4ThreeVector norm 379 G4ThreeVector norm(0.,0.,0.); 392 380 393 381 // Calculate distances as if in 1st octant … … 400 388 // normals 401 389 402 const G4double delta = 0.5*kCarTolerance;390 static const G4double delta = 0.5*kCarTolerance; 403 391 const G4ThreeVector nX = G4ThreeVector( 1.0, 0,0 ); 404 392 const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0 ); … … 415 403 { 416 404 noSurfaces ++; 417 if ( p.x() >= 0.){ // on +X surface 418 normX= nX ; // G4ThreeVector( 1.0, 0., 0. ); 419 }else{ 420 normX= nmX; // G4ThreeVector(-1.0, 0., 0. ); 421 } 405 if ( p.x() >= 0. ) { normX= nX ; } // on +X surface : (1,0,0) 406 else { normX= nmX; } // (-1,0,0) 422 407 sumnorm= normX; 423 408 } … … 426 411 { 427 412 noSurfaces ++; 428 if ( p.y() >= 0.){ // on +Y surface 429 normY= nY; 430 }else{ 431 normY = nmY; 432 } 413 if ( p.y() >= 0. ) { normY= nY; } // on +Y surface 414 else { normY= nmY; } 433 415 sumnorm += normY; 434 416 } … … 437 419 { 438 420 noSurfaces ++; 439 if ( p.z() >= 0.){ // on +Z surface 440 normZ= nZ; 441 }else{ 442 normZ = nmZ; 443 } 444 sumnorm += normZ; 445 } 446 447 // sumnorm= normX + normY + normZ; 448 const G4double invSqrt2 = 1.0 / std::sqrt( 2.0); 449 const G4double invSqrt3 = 1.0 / std::sqrt( 3.0); 450 451 norm= G4ThreeVector( 0., 0., 0.); 421 if ( p.z() >= 0. ) { normZ= nZ; } // on +Z surface 422 else { normZ= nmZ; } 423 sumnorm += normZ; 424 } 425 426 static const G4double invSqrt2 = 1.0 / std::sqrt(2.0); 427 static const G4double invSqrt3 = 1.0 / std::sqrt(3.0); 428 452 429 if( noSurfaces > 0 ) 453 430 { 454 if( noSurfaces == 1 ){ 431 if( noSurfaces == 1 ) 432 { 455 433 norm= sumnorm; 456 }else{ 434 } 435 else 436 { 457 437 // norm = sumnorm . unit(); 458 if( noSurfaces == 2 ) { 438 if( noSurfaces == 2 ) 439 { 459 440 // 2 surfaces -> on edge 460 441 norm = invSqrt2 * sumnorm; 461 } else { 442 } 443 else 444 { 462 445 // 3 surfaces (on corner) 463 446 norm = invSqrt3 * sumnorm; 464 447 } 465 448 } 466 }else{ 449 } 450 else 451 { 467 452 #ifdef G4CSGDEBUG 468 453 G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning, … … 483 468 { 484 469 G4double distx, disty, distz ; 485 G4ThreeVector norm 470 G4ThreeVector norm(0.,0.,0.); 486 471 487 472 // Calculate distances as if in 1st octant … … 495 480 if ( distx <= distz ) // Closest to X 496 481 { 497 if ( p.x() < 0 ) norm = G4ThreeVector(-1.0,0,0) ;498 else norm = G4ThreeVector( 1.0,0,0) ;482 if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; } 483 else { norm = G4ThreeVector( 1.0,0,0) ; } 499 484 } 500 485 else // Closest to Z 501 486 { 502 if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;503 else norm = G4ThreeVector(0,0, 1.0) ;487 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; } 488 else { norm = G4ThreeVector(0,0, 1.0) ; } 504 489 } 505 490 } … … 508 493 if ( disty <= distz ) // Closest to Y 509 494 { 510 if ( p.y() < 0 ) norm = G4ThreeVector(0,-1.0,0) ;511 else norm = G4ThreeVector(0, 1.0,0) ;495 if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; } 496 else { norm = G4ThreeVector(0, 1.0,0) ; } 512 497 } 513 498 else // Closest to Z 514 499 { 515 if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;516 else norm = G4ThreeVector(0,0, 1.0) ;500 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; } 501 else { norm = G4ThreeVector(0,0, 1.0) ; } 517 502 } 518 503 } … … 541 526 // shape. 542 527 543 G4double G4Box::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const 528 G4double G4Box::DistanceToIn(const G4ThreeVector& p, 529 const G4ThreeVector& v) const 544 530 { 545 531 G4double safx, safy, safz ; … … 549 535 G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ; 550 536 537 static const G4double delta = 0.5*kCarTolerance; 538 551 539 safx = std::fabs(p.x()) - fDx ; // minimum distance to x surface of shape 552 540 safy = std::fabs(p.y()) - fDy ; … … 558 546 // travel is in a direction away from the shape. 559 547 560 if ( ((p.x()*v.x() >= 0.0) && safx > -kCarTolerance*0.5)561 || ((p.y()*v.y() >= 0.0) && safy > -kCarTolerance*0.5)562 || ((p.z()*v.z() >= 0.0) && safz > -kCarTolerance*0.5) )548 if ( ((p.x()*v.x() >= 0.0) && (safx > -delta)) 549 || ((p.y()*v.y() >= 0.0) && (safy > -delta)) 550 || ((p.z()*v.z() >= 0.0) && (safz > -delta)) ) 563 551 { 564 552 return kInfinity ; // travel away or parallel within tolerance … … 568 556 // X Planes 569 557 570 if ( v.x() )558 if ( v.x() ) // != 0 571 559 { 572 560 stmp = 1.0/std::fabs(v.x()) ; … … 579 567 else 580 568 { 581 if (v.x() > 0) sOut = (fDx - p.x())*stmp ;582 if (v.x() < 0) sOut = (fDx + p.x())*stmp ;569 if (v.x() < 0) { sOut = (fDx + p.x())*stmp ; } 570 else { sOut = (fDx - p.x())*stmp ; } 583 571 } 584 572 } … … 586 574 // Y Planes 587 575 588 if ( v.y() )576 if ( v.y() ) // != 0 589 577 { 590 578 stmp = 1.0/std::fabs(v.y()) ; … … 595 583 smaxy = (fDy+std::fabs(p.y()))*stmp ; 596 584 597 if (sminy > smin) smin=sminy ;598 if (smaxy < smax) smax=smaxy ;599 600 if (smin >= smax-kCarTolerance*0.5)585 if (sminy > smin) { smin=sminy ; } 586 if (smaxy < smax) { smax=smaxy ; } 587 588 if (smin >= (smax-delta)) 601 589 { 602 590 return kInfinity ; // touch XY corner … … 605 593 else 606 594 { 607 if (v.y() > 0) sOuty = (fDy - p.y())*stmp ;608 if (v.y() < 0) sOuty = (fDy + p.y())*stmp ;609 if( sOuty < sOut ) sOut = sOuty ;595 if (v.y() < 0) { sOuty = (fDy + p.y())*stmp ; } 596 else { sOuty = (fDy - p.y())*stmp ; } 597 if( sOuty < sOut ) { sOut = sOuty ; } 610 598 } 611 599 } … … 613 601 // Z planes 614 602 615 if ( v.z() ) 603 if ( v.z() ) // != 0 616 604 { 617 605 stmp = 1.0/std::fabs(v.z()) ; 618 606 619 if ( safz >= 0.0 )607 if ( safz >= 0.0 ) 620 608 { 621 609 sminz = safz*stmp ; 622 610 smaxz = (fDz+std::fabs(p.z()))*stmp ; 623 611 624 if (sminz > smin) smin = sminz ;625 if (smaxz < smax) smax = smaxz ;626 627 if (smin >= smax-kCarTolerance*0.5)612 if (sminz > smin) { smin = sminz ; } 613 if (smaxz < smax) { smax = smaxz ; } 614 615 if (smin >= (smax-delta)) 628 616 { 629 617 return kInfinity ; // touch ZX or ZY corners … … 632 620 else 633 621 { 634 if (v.z() > 0) sOutz = (fDz - p.z())*stmp ;635 if (v.z() < 0) sOutz = (fDz + p.z())*stmp ;636 if( sOutz < sOut ) sOut = sOutz ;637 } 638 } 639 640 if ( sOut <= smin + 0.5*kCarTolerance) // travel over edge622 if (v.z() < 0) { sOutz = (fDz + p.z())*stmp ; } 623 else { sOutz = (fDz - p.z())*stmp ; } 624 if( sOutz < sOut ) { sOut = sOutz ; } 625 } 626 } 627 628 if (sOut <= (smin + delta)) // travel over edge 641 629 { 642 630 return kInfinity ; 643 631 } 644 if (smin < 0.5*kCarTolerance) smin = 0.0 ;632 if (smin < delta) { smin = 0.0 ; } 645 633 646 634 return smin ; … … 662 650 safez = std::fabs(p.z()) - fDz ; 663 651 664 if (safex > safe) safe = safex ;665 if (safey > safe) safe = safey ;666 if (safez > safe) safe = safez ;652 if (safex > safe) { safe = safex ; } 653 if (safey > safe) { safe = safey ; } 654 if (safez > safe) { safe = safez ; } 667 655 668 656 return safe ; … … 671 659 ///////////////////////////////////////////////////////////////////////// 672 660 // 673 // Calc luate distance to surface of box from inside661 // Calculate distance to surface of box from inside 674 662 // by calculating distances to box's x/y/z planes. 675 663 // Smallest distance is exact distance to exiting. … … 682 670 { 683 671 ESide side = kUndefined ; 684 G4double pdist,stmp,snxt; 685 686 if (calcNorm) *validNorm = true ; // All normals are valid 687 688 if (v.x() > 0) // X planes 672 G4double pdist,stmp,snxt=kInfinity; 673 674 static const G4double delta = 0.5*kCarTolerance; 675 676 if (calcNorm) { *validNorm = true ; } // All normals are valid 677 678 if (v.x() > 0) // X planes 689 679 { 690 680 pdist = fDx - p.x() ; 691 681 692 if (pdist > kCarTolerance*0.5)682 if (pdist > delta) 693 683 { 694 684 snxt = pdist/v.x() ; … … 697 687 else 698 688 { 699 if (calcNorm) *n = G4ThreeVector(1,0,0) ;700 return 701 } 702 } 703 else if (v.x() < 0) 689 if (calcNorm) { *n = G4ThreeVector(1,0,0) ; } 690 return snxt = 0 ; 691 } 692 } 693 else if (v.x() < 0) 704 694 { 705 695 pdist = fDx + p.x() ; 706 696 707 if (pdist > kCarTolerance*0.5)697 if (pdist > delta) 708 698 { 709 699 snxt = -pdist/v.x() ; … … 712 702 else 713 703 { 714 if (calcNorm) *n = G4ThreeVector(-1,0,0) ;704 if (calcNorm) { *n = G4ThreeVector(-1,0,0) ; } 715 705 return snxt = 0 ; 716 706 } 717 707 } 718 else snxt = kInfinity ; 719 720 if ( v.y() > 0 ) // Y planes 721 { 722 pdist=fDy-p.y(); 723 724 if (pdist>kCarTolerance*0.5) 725 { 726 stmp=pdist/v.y(); 727 728 if (stmp<snxt) 729 { 730 snxt=stmp; 731 side=kPY; 732 } 733 } 734 else 735 { 736 if (calcNorm) *n = G4ThreeVector(0,1,0) ; 737 return snxt = 0 ; 738 } 739 } 740 else if ( v.y() < 0 ) 741 { 742 pdist = fDy + p.y() ; 743 744 if (pdist > kCarTolerance*0.5) 745 { 746 stmp=-pdist/v.y(); 747 748 if (stmp<snxt) 749 { 750 snxt=stmp; 751 side=kMY; 752 } 753 } 754 else 755 { 756 if (calcNorm) *n = G4ThreeVector(0,-1,0) ; 757 return snxt = 0 ; 758 } 759 } 760 if (v.z()>0) // Z planes 761 { 762 pdist=fDz-p.z(); 763 764 if (pdist > kCarTolerance*0.5) 765 { 766 stmp=pdist/v.z(); 708 709 if (v.y() > 0) // Y planes 710 { 711 pdist = fDy-p.y(); 712 713 if (pdist > delta) 714 { 715 stmp = pdist/v.y(); 767 716 768 717 if (stmp < snxt) 769 718 { 770 snxt =stmp;771 side =kPZ;719 snxt = stmp; 720 side = kPY; 772 721 } 773 722 } 774 723 else 775 724 { 776 if (calcNorm) *n = G4ThreeVector(0,0,1) ;777 return 778 } 779 } 780 else if (v. z()<0)781 { 782 pdist = fD z + p.z() ;783 784 if (pdist > kCarTolerance*0.5)785 { 786 stmp =-pdist/v.z();787 788 if ( stmp < snxt)725 if (calcNorm) { *n = G4ThreeVector(0,1,0) ; } 726 return snxt = 0 ; 727 } 728 } 729 else if (v.y() < 0) 730 { 731 pdist = fDy + p.y() ; 732 733 if (pdist > delta) 734 { 735 stmp = -pdist/v.y(); 736 737 if ( stmp < snxt ) 789 738 { 790 snxt =stmp;791 side =kMZ;739 snxt = stmp; 740 side = kMY; 792 741 } 793 742 } 794 743 else 795 744 { 796 if (calcNorm) *n = G4ThreeVector(0,0,-1) ; 797 return snxt = 0 ; 798 } 799 } 745 if (calcNorm) { *n = G4ThreeVector(0,-1,0) ; } 746 return snxt = 0 ; 747 } 748 } 749 750 if (v.z() > 0) // Z planes 751 { 752 pdist = fDz-p.z(); 753 754 if ( pdist > delta ) 755 { 756 stmp = pdist/v.z(); 757 758 if ( stmp < snxt ) 759 { 760 snxt = stmp; 761 side = kPZ; 762 } 763 } 764 else 765 { 766 if (calcNorm) { *n = G4ThreeVector(0,0,1) ; } 767 return snxt = 0 ; 768 } 769 } 770 else if (v.z() < 0) 771 { 772 pdist = fDz + p.z(); 773 774 if ( pdist > delta ) 775 { 776 stmp = -pdist/v.z(); 777 778 if ( stmp < snxt ) 779 { 780 snxt = stmp; 781 side = kMZ; 782 } 783 } 784 else 785 { 786 if (calcNorm) { *n = G4ThreeVector(0,0,-1) ; } 787 return snxt = 0 ; 788 } 789 } 790 800 791 if (calcNorm) 801 792 { … … 875 866 // shortest Dist to any boundary now MIN(safx1,safx2,safy1..) 876 867 877 if (safx2 < safx1) safe = safx2 ;878 else safe = safx1 ;879 if (safy1 < safe) safe = safy1 ;880 if (safy2 < safe) safe = safy2 ;881 if (safz1 < safe) safe = safz1 ;882 if (safz2 < safe) safe = safz2 ;883 884 if (safe < 0) safe = 0 ;868 if (safx2 < safx1) { safe = safx2; } 869 else { safe = safx1; } 870 if (safy1 < safe) { safe = safy1; } 871 if (safy2 < safe) { safe = safy2; } 872 if (safz1 < safe) { safe = safz1; } 873 if (safz2 < safe) { safe = safz2; } 874 875 if (safe < 0) { safe = 0 ; } 885 876 return safe ; 886 877 } … … 979 970 py = -fDy +2*fDy*G4UniformRand(); 980 971 981 if(G4UniformRand() > 0.5) pz = fDz;982 else pz = -fDz;972 if(G4UniformRand() > 0.5) { pz = fDz; } 973 else { pz = -fDz; } 983 974 } 984 975 else if ( ( select - Sxy ) < Sxz ) … … 987 978 pz = -fDz +2*fDz*G4UniformRand(); 988 979 989 if(G4UniformRand() > 0.5) py = fDy;990 else py = -fDy;980 if(G4UniformRand() > 0.5) { py = fDy; } 981 else { py = -fDy; } 991 982 } 992 983 else … … 995 986 pz = -fDz +2*fDz*G4UniformRand(); 996 987 997 if(G4UniformRand() > 0.5) px = fDx;998 else px = -fDx;988 if(G4UniformRand() > 0.5) { px = fDx; } 989 else { px = -fDx; } 999 990 } 1000 991 return G4ThreeVector(px,py,pz); -
trunk/source/geometry/solids/CSG/src/G4Orb.cc
r1228 r1315 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4Orb.cc,v 1.3 0 2009/11/30 10:20:38 gcosmoExp $27 // GEANT4 tag $Name: geant4-09-0 3$26 // $Id: G4Orb.cc,v 1.31 2009/12/04 15:39:56 grichine Exp $ 27 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 28 28 // 29 29 // class G4Orb … … 296 296 297 297 298 rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z() 298 rad2 = p.x()*p.x()+p.y()*p.y()+p.z()*p.z(); 299 299 300 300 G4double rad = std::sqrt(rad2); … … 304 304 // sets `in' 305 305 306 tolRMax = fRmax - fRmaxTolerance*0.5 306 tolRMax = fRmax - fRmaxTolerance*0.5; 307 307 308 if ( rad <= tolRMax ) { in = kInside 308 if ( rad <= tolRMax ) { in = kInside; } 309 309 else 310 310 { 311 tolRMax = fRmax + fRmaxTolerance*0.5 312 if ( rad <= tolRMax ) { in = kSurface 313 else { in = kOutside 311 tolRMax = fRmax + fRmaxTolerance*0.5; 312 if ( rad <= tolRMax ) { in = kSurface; } 313 else { in = kOutside; } 314 314 } 315 315 return in; … … 357 357 const G4ThreeVector& v ) const 358 358 { 359 G4double snxt = kInfinity 360 361 G4double rad 2, pDotV3d, tolORMax2, tolIRMax2;362 G4double c, d2, s = kInfinity 359 G4double snxt = kInfinity; // snxt = default return value 360 361 G4double rad, pDotV3d; // , tolORMax2, tolIRMax2; 362 G4double c, d2, s = kInfinity; 363 363 364 364 const G4double dRmax = 100.*fRmax; … … 366 366 // General Precalcs 367 367 368 rad 2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z();369 pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() 368 rad = std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()); 369 pDotV3d = p.x()*v.x() + p.y()*v.y() + p.z()*v.z(); 370 370 371 371 // Radial Precalcs 372 372 373 tolORMax2 = (fRmax+fRmaxTolerance*0.5)*(fRmax+fRmaxTolerance*0.5);374 tolIRMax2 = (fRmax-fRmaxTolerance*0.5)*(fRmax-fRmaxTolerance*0.5);373 // tolORMax2 = (fRmax+fRmaxTolerance*0.5)*(fRmax+fRmaxTolerance*0.5); 374 // tolIRMax2 = (fRmax-fRmaxTolerance*0.5)*(fRmax-fRmaxTolerance*0.5); 375 375 376 376 // Outer spherical shell intersection … … 388 388 // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2)) 389 389 390 391 G4double rad = std::sqrt(rad2);392 390 c = (rad - fRmax)*(rad + fRmax); 393 391 394 if ( c > fRmaxTolerance*fRmax ) 395 { 396 // If outside tolerant boundary of outer G4Orb 397 // [ should be std::sqrt(rad2) - fRmax > fRmaxTolerance*0.5 ] 398 399 d2 = pDotV3d*pDotV3d - c ; 400 401 if ( d2 >= 0 ) 402 { 403 s = -pDotV3d - std::sqrt(d2) ; 404 if ( s >= 0 ) 405 { 406 if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on 407 { // 64 bits systems. Split long distances and recompute 408 G4double fTerm = s-std::fmod(s,dRmax); 409 s = fTerm + DistanceToIn(p+fTerm*v,v); 410 } 411 return snxt = s; 412 } 413 } 414 else // No intersection with G4Orb 415 { 416 return snxt = kInfinity; 417 } 418 } 419 else 420 { 421 if ( c > -fRmaxTolerance*fRmax ) // on surface 422 { 423 d2 = pDotV3d*pDotV3d - c ; 424 if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) ) 392 if( rad > fRmax-fRmaxTolerance*0.5 ) // not inside in terms of Inside(p) 393 { 394 if ( c > fRmaxTolerance*fRmax ) 395 { 396 // If outside tolerant boundary of outer G4Orb in terms of c 397 // [ should be std::sqrt(rad2) - fRmax > fRmaxTolerance*0.5 ] 398 399 d2 = pDotV3d*pDotV3d - c; 400 401 if ( d2 >= 0 ) 402 { 403 s = -pDotV3d - std::sqrt(d2); 404 if ( s >= 0 ) 405 { 406 if ( s > dRmax ) // Avoid rounding errors due to precision issues seen on 407 { // 64 bits systems. Split long distances and recompute 408 G4double fTerm = s - std::fmod(s,dRmax); 409 s = fTerm + DistanceToIn(p+fTerm*v,v); 410 } 411 return snxt = s; 412 } 413 } 414 else // No intersection with G4Orb 425 415 { 426 416 return snxt = kInfinity; 427 417 } 428 else 429 { 430 return snxt = 0.; 431 } 432 } 418 } 419 else // not outside in terms of c 420 { 421 if ( c > -fRmaxTolerance*fRmax ) // on surface 422 { 423 d2 = pDotV3d*pDotV3d - c; 424 if ( (d2 < fRmaxTolerance*fRmax) || (pDotV3d >= 0) ) 425 { 426 return snxt = kInfinity; 427 } 428 else 429 { 430 return snxt = 0.; 431 } 432 } 433 } 434 } 433 435 #ifdef G4CSGDEBUG 434 435 436 else // inside ??? 437 { 436 438 G4Exception("G4Orb::DistanceToIn(p,v)", "Notification", 437 439 JustWarning, "Point p is inside !?"); 438 440 } 439 441 #endif 440 } 442 441 443 return snxt; 442 444 } … … 520 522 if(calcNorm) 521 523 { 522 *validNorm = true 523 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) 524 *validNorm = true; 525 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax); 524 526 } 525 527 return snxt = 0; … … 528 530 { 529 531 snxt = -pDotV3d + std::sqrt(d2); // second root since inside Rmax 530 side = kRMax 532 side = kRMax; 531 533 } 532 534 } … … 596 598 if( Inside(p) == kOutside ) 597 599 { 598 G4cout.precision(16) 599 G4cout << G4endl 600 G4cout.precision(16); 601 G4cout << G4endl; 600 602 DumpInfo(); 601 G4cout << "Position:" << G4endl << G4endl 602 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl 603 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl 604 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl 603 G4cout << "Position:" << G4endl << G4endl; 604 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl; 605 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl; 606 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl; 605 607 G4Exception("G4Orb::DistanceToOut(p)", "Notification", JustWarning, 606 608 "Point p is outside !?" );
Note: See TracChangeset
for help on using the changeset viewer.