Changeset 921 for trunk/source/geometry/solids/CSG/src
- Timestamp:
- Feb 16, 2009, 10:14:30 AM (16 years ago)
- Location:
- trunk/source/geometry/solids/CSG/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/geometry/solids/CSG/src/G4Cons.cc
r850 r921 25 25 // 26 26 // 27 // $Id: G4Cons.cc,v 1. 56 2008/02/20 08:56:16gcosmo Exp $28 // GEANT4 tag $Name: HEAD$27 // $Id: G4Cons.cc,v 1.60 2008/11/06 15:26:53 gcosmo Exp $ 28 // GEANT4 tag $Name: geant4-09-02-cand-01 $ 29 29 // 30 30 // … … 87 87 88 88 if ( pDz > 0 ) 89 fDz = pDz ; 89 { 90 fDz = pDz; 91 } 90 92 else 91 93 { … … 99 101 // Check radii 100 102 101 if ( pRmin1 < pRmax1 && pRmin2 < pRmax2 && pRmin1 >= 0 && pRmin2 >= 0)103 if ( (pRmin1<pRmax1) && (pRmin2<pRmax2) && (pRmin1>=0) && (pRmin2>=0) ) 102 104 { 103 105 … … 106 108 fRmin2 = pRmin2 ; 107 109 fRmax2 = pRmax2 ; 108 if( (pRmin1 == 0.0 && pRmin2 > 0.0) ) fRmin1 = 1e3*kRadTolerance ;109 if( (pRmin2 == 0.0 && pRmin1 > 0.0) ) fRmin2 = 1e3*kRadTolerance ;110 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; } 111 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; } 110 112 } 111 113 else … … 119 121 } 120 122 121 // Check angles 122 123 if ( pDPhi >= twopi ) 123 fPhiFullCone = true; 124 if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles 124 125 { 125 126 fDPhi=twopi; … … 128 129 else 129 130 { 130 if ( pDPhi > 0 ) fDPhi = pDPhi ; 131 fPhiFullCone = false; 132 if ( pDPhi > 0 ) 133 { 134 fDPhi = pDPhi; 135 } 131 136 else 132 137 { … … 135 140 << pDPhi << G4endl; 136 141 G4Exception("G4Cons::G4Cons()", "InvalidSetup", 137 FatalException, "Invalid pDPhi.") ; 138 } 139 140 // Ensure pSPhi in 0-2PI or -2PI-0 range if shape crosses 0 141 142 if ( pSPhi < 0 ) fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi) ; 143 else fSPhi = std::fmod(pSPhi,twopi) ; 144 145 if (fSPhi + fDPhi > twopi) fSPhi -= twopi ; 146 } 142 FatalException, "Invalid dphi."); 143 } 144 145 // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0 146 147 if ( pSPhi < 0 ) 148 { 149 fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi); 150 } 151 else 152 { 153 fSPhi = std::fmod(pSPhi,twopi) ; 154 } 155 if ( fSPhi+fDPhi > twopi ) 156 { 157 fSPhi -= twopi ; 158 } 159 } 160 InitializeTrigonometry(); 147 161 } 148 162 … … 173 187 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ; 174 188 EInside in; 175 176 if (std::fabs(p.z()) > fDz + kCarTolerance*0.5 ) return in = kOutside; 177 else if(std::fabs(p.z()) >= fDz - kCarTolerance*0.5 ) in = kSurface; 178 else in = kInside; 189 static const G4double halfCarTolerance=kCarTolerance*0.5; 190 static const G4double halfRadTolerance=kRadTolerance*0.5; 191 static const G4double halfAngTolerance=kAngTolerance*0.5; 192 193 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; } 194 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; } 195 else { in = kInside; } 179 196 180 197 r2 = p.x()*p.x() + p.y()*p.y() ; … … 184 201 // rh2 = rh*rh; 185 202 186 tolRMin = rl - kRadTolerance*0.5;187 if ( tolRMin < 0 ) tolRMin = 0 ;188 tolRMax = rh + kRadTolerance*0.5;189 190 if ( r2 < tolRMin*tolRMin || r2 > tolRMax*tolRMax ) return in = kOutside;191 192 if (rl) tolRMin = rl + kRadTolerance*0.5 ;193 else tolRMin = 0.0 ;194 tolRMax = rh - kRadTolerance*0.5;203 tolRMin = rl - halfRadTolerance; 204 if ( tolRMin < 0 ) { tolRMin = 0; } 205 tolRMax = rh + halfRadTolerance; 206 207 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; } 208 209 if (rl) { tolRMin = rl + halfRadTolerance; } 210 else { tolRMin = 0.0; } 211 tolRMax = rh - halfRadTolerance; 195 212 196 213 if (in == kInside) // else it's kSurface already 197 214 { 198 if (r2 < tolRMin*tolRMin || r2 >= tolRMax*tolRMax) in = kSurface; 199 // if (r2 <= tolRMin*tolRMin || r2-rh2 >= -rh*kRadTolerance) in = kSurface; 200 } 201 if ( ( fDPhi < twopi - kAngTolerance ) && 202 ( (p.x() != 0.0 ) || (p.y() != 0.0) ) ) 215 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; } 216 } 217 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) ) 203 218 { 204 219 pPhi = std::atan2(p.y(),p.x()) ; 205 220 206 if ( pPhi < fSPhi - kAngTolerance*0.5 ) pPhi += twopi ;207 else if ( pPhi > fSPhi + fDPhi + kAngTolerance*0.5 ) pPhi -= twopi;221 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; } 222 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; } 208 223 209 if ( (pPhi < fSPhi - kAngTolerance*0.5) ||210 (pPhi > fSPhi + fDPhi + kAngTolerance*0.5) ) return in = kOutside;224 if ( (pPhi < fSPhi - halfAngTolerance) || 225 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; } 211 226 212 227 else if (in == kInside) // else it's kSurface anyway already 213 228 { 214 if ( (pPhi < fSPhi + kAngTolerance*0.5) ||215 (pPhi > fSPhi + fDPhi - kAngTolerance*0.5) ) in = kSurface ;216 } 217 } 218 else if ( fDPhi < twopi - kAngTolerance ) in = kSurface ;229 if ( (pPhi < fSPhi + halfAngTolerance) || 230 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; } 231 } 232 } 233 else if ( !fPhiFullCone ) { in = kSurface; } 219 234 220 235 return in ; … … 244 259 G4double& pMax ) const 245 260 { 246 if ( !pTransform.IsRotated() && 247 fDPhi == twopi && fRmin1 == 0 && fRmin2 == 0)261 if ( !pTransform.IsRotated() && (fDPhi == twopi) 262 && (fRmin1 == 0) && (fRmin2 == 0) ) 248 263 { 249 264 // Special case handling for unrotated solid cones … … 265 280 if (pVoxelLimit.IsZLimited()) 266 281 { 267 if( zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance||268 zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance)282 if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) || 283 (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) ) 269 284 { 270 285 return false ; … … 290 305 if (pVoxelLimit.IsXLimited()) 291 306 { 292 if ( xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance||293 xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance)307 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) || 308 (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) ) 294 309 { 295 310 return false ; … … 315 330 if (pVoxelLimit.IsYLimited()) 316 331 { 317 if ( yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance||318 yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance)332 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) || 333 (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) ) 319 334 { 320 335 return false ; … … 338 353 yoff2 = yMax - yoffset ; 339 354 340 if ( yoff1 >= 0 && yoff2 >= 0) // Y limits cross max/min x => no change341 { 355 if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x 356 { // => no change 342 357 pMin = xMin ; 343 358 pMax = xMax ; … … 362 377 xoff2 = xMax - xoffset ; 363 378 364 if ( xoff1 >= 0 && xoff2 >= 0 ) // X limits cross max/min y => no change365 { 379 if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y 380 { // => no change 366 381 pMin = yMin ; 367 382 pMax = yMax ; … … 409 424 for ( i = 0 ; i < noEntries ; i += 4 ) 410 425 { 411 ClipCrossSection(vertices, i,pVoxelLimit,pAxis,pMin,pMax) ;426 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ; 412 427 } 413 428 for ( i = 0 ; i < noBetweenSections4 ; i += 4 ) 414 429 { 415 ClipBetweenSections(vertices, i,pVoxelLimit,pAxis,pMin,pMax) ;430 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ; 416 431 } 417 if ( pMin != kInfinity || pMax != -kInfinity)432 if ( (pMin != kInfinity) || (pMax != -kInfinity) ) 418 433 { 419 434 existsAfterClip = true ; … … 462 477 G4double tanRMin, secRMin, pRMin, widRMin; 463 478 G4double tanRMax, secRMax, pRMax, widRMax; 464 G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance; 479 480 static const G4double delta = 0.5*kCarTolerance; 481 static const G4double dAngle = 0.5*kAngTolerance; 465 482 466 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1. 0);483 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.); 467 484 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe; 468 485 … … 482 499 distRMax = std::fabs(pRMax - widRMax)/secRMax; 483 500 484 if ( fDPhi < twopi) // && rho )// Protected against (0,0,z)501 if (!fPhiFullCone) // Protected against (0,0,z) 485 502 { 486 503 if ( rho ) … … 488 505 pPhi = std::atan2(p.y(),p.x()); 489 506 490 if (pPhi < fSPhi-delta) pPhi += twopi;491 else if (pPhi > fSPhi+fDPhi+delta) pPhi -= twopi;507 if (pPhi < fSPhi-delta) { pPhi += twopi; } 508 else if (pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; } 492 509 493 510 distSPhi = std::fabs( pPhi - fSPhi ); 494 distEPhi = std::fabs( pPhi - fSPhi - fDPhi);511 distEPhi = std::fabs( pPhi - fSPhi - fDPhi ); 495 512 } 496 513 else if( !(fRmin1) || !(fRmin2) ) … … 499 516 distEPhi = 0.; 500 517 } 501 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi),0);502 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi),0);518 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0); 519 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0); 503 520 } 504 521 if ( rho > delta ) 505 522 { 506 nR = G4ThreeVector(p.x()/rho/secRMax,p.y()/rho/secRMax,-tanRMax/secRMax); 507 if (fRmin1 || fRmin2) nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin); 523 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax); 524 if (fRmin1 || fRmin2) 525 { 526 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin); 527 } 508 528 } 509 529 … … 513 533 sumnorm += nR; 514 534 } 515 if( (fRmin1 || fRmin2) && distRMin <= delta)535 if( (fRmin1 || fRmin2) && (distRMin <= delta) ) 516 536 { 517 537 noSurfaces ++; 518 538 sumnorm += nr; 519 539 } 520 if( fDPhi < twopi)540 if( !fPhiFullCone ) 521 541 { 522 542 if (distSPhi <= dAngle) … … 534 554 { 535 555 noSurfaces ++; 536 if ( p.z() >= 0.) sumnorm += nZ;537 else sumnorm -= nZ;556 if ( p.z() >= 0.) { sumnorm += nZ; } 557 else { sumnorm -= nZ; } 538 558 } 539 559 if ( noSurfaces == 0 ) … … 545 565 norm = ApproxSurfaceNormal(p); 546 566 } 547 else if ( noSurfaces == 1 ) norm = sumnorm; 548 else norm = sumnorm.unit(); 567 else if ( noSurfaces == 1 ) { norm = sumnorm; } 568 else { norm = sumnorm.unit(); } 569 549 570 return norm ; 550 571 } 551 572 552 //////////////////////////////////////////////////////////////////////////// //////573 //////////////////////////////////////////////////////////////////////////// 553 574 // 554 575 // Algorithm for SurfaceNormal() following the original specification … … 605 626 } 606 627 } 607 if ( fDPhi < twopi&& rho ) // Protected against (0,0,z)628 if ( !fPhiFullCone && rho ) // Protected against (0,0,z) 608 629 { 609 630 phi = std::atan2(p.y(),p.x()) ; 610 631 611 if (phi < 0) phi += twopi ;612 613 if (fSPhi < 0) distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;614 else distSPhi = std::fabs(phi - fSPhi)*rho ;632 if (phi < 0) { phi += twopi; } 633 634 if (fSPhi < 0) { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; } 635 else { distSPhi = std::fabs(phi - fSPhi)*rho; } 615 636 616 637 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ; … … 620 641 if (distSPhi < distEPhi) 621 642 { 622 if (distSPhi < distMin) side = kNSPhi ;643 if (distSPhi < distMin) { side = kNSPhi; } 623 644 } 624 645 else 625 646 { 626 if (distEPhi < distMin) side = kNEPhi ;647 if (distEPhi < distMin) { side = kNEPhi; } 627 648 } 628 649 } … … 631 652 case kNRMin: // Inner radius 632 653 rho *= secRMin ; 633 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho,tanRMin/secRMin) ;654 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ; 634 655 break ; 635 656 case kNRMax: // Outer radius 636 657 rho *= secRMax ; 637 norm = G4ThreeVector(p.x()/rho, p.y()/rho,-tanRMax/secRMax) ;658 norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ; 638 659 break ; 639 660 case kNZ: // +/- dz 640 if (p.z() > 0) norm = G4ThreeVector(0,0,1) ;641 else norm = G4ThreeVector(0,0,-1) ;661 if (p.z() > 0) { norm = G4ThreeVector(0,0,1); } 662 else { norm = G4ThreeVector(0,0,-1); } 642 663 break ; 643 664 case kNSPhi: 644 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi),0) ;665 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ; 645 666 break ; 646 667 case kNEPhi: 647 norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi),0) ;668 norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ; 648 669 break ; 649 670 default: … … 677 698 // 678 699 // NOTE: 679 // - Precalculations for phi trigonometry are Done `just in time'680 700 // - `if valid' implies tolerant checking of intersection points 681 701 // - z, phi intersection from Tubs … … 686 706 G4double snxt = kInfinity ; // snxt = default return value 687 707 688 G4bool seg ; // true if segmented in phi 689 G4double hDPhi,hDPhiOT,hDPhiIT,cosHDPhiOT=0.,cosHDPhiIT=0. ; 690 // half dphi + outer tolerance 691 G4double cPhi,sinCPhi=0.,cosCPhi=0. ; // central phi 708 static const G4double halfCarTolerance=kCarTolerance*0.5; 709 static const G4double halfRadTolerance=kRadTolerance*0.5; 692 710 693 711 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; // Data for cones … … 704 722 G4double nt1,nt2,nt3 ; 705 723 G4double Comp ; 706 G4double cosSPhi,sinSPhi ; // Trig for phi start intersect707 G4double ePhi,cosEPhi,sinEPhi ; // for phi end intersect708 709 //710 // Set phi divided flag and precalcs711 //712 if (fDPhi < twopi)713 {714 seg = true ;715 hDPhi = 0.5*fDPhi ; // half delta phi716 cPhi = fSPhi + hDPhi ; ;717 hDPhiOT = hDPhi + 0.5*kAngTolerance ; // outers tol' half delta phi718 hDPhiIT = hDPhi - 0.5*kAngTolerance ;719 sinCPhi = std::sin(cPhi) ;720 cosCPhi = std::cos(cPhi) ;721 cosHDPhiOT = std::cos(hDPhiOT) ;722 cosHDPhiIT = std::cos(hDPhiIT) ;723 }724 else seg = false ;725 724 726 725 // Cone Precalcs … … 730 729 rMinAv = (fRmin1 + fRmin2)*0.5 ; 731 730 732 if (rMinAv > kRadTolerance*0.5)733 { 734 rMinOAv = rMinAv - kRadTolerance*0.5;735 rMinIAv = rMinAv + kRadTolerance*0.5;731 if (rMinAv > halfRadTolerance) 732 { 733 rMinOAv = rMinAv - halfRadTolerance ; 734 rMinIAv = rMinAv + halfRadTolerance ; 736 735 } 737 736 else … … 743 742 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 744 743 rMaxAv = (fRmax1 + fRmax2)*0.5 ; 745 rMaxOAv = rMaxAv + kRadTolerance*0.5;744 rMaxOAv = rMaxAv + halfRadTolerance ; 746 745 747 746 // Intersection with z-surfaces 748 747 749 tolIDz = fDz - kCarTolerance*0.5;750 tolODz = fDz + kCarTolerance*0.5;748 tolIDz = fDz - halfCarTolerance ; 749 tolODz = fDz + halfCarTolerance ; 751 750 752 751 if (std::fabs(p.z()) >= tolIDz) … … 756 755 s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance 757 756 758 if( s < 0.0 ) s = 0.0 ;// negative dist -> zero757 if( s < 0.0 ) { s = 0.0; } // negative dist -> zero 759 758 760 759 xi = p.x() + s*v.x() ; // Intersection coords 761 760 yi = p.y() + s*v.y() ; 762 rhoi2 = xi*xi + yi*yi 761 rhoi2 = xi*xi + yi*yi ; 763 762 764 763 // Check validity of intersection … … 767 766 if (v.z() > 0) 768 767 { 769 tolORMin = fRmin1 - 0.5*kRadTolerance*secRMin ;770 tolIRMin = fRmin1 + 0.5*kRadTolerance*secRMin ;771 tolIRMax = fRmax1 - 0.5*kRadTolerance*secRMin ;772 tolORMax2 = (fRmax1 + 0.5*kRadTolerance*secRMax)*773 (fRmax1 + 0.5*kRadTolerance*secRMax) ;768 tolORMin = fRmin1 - halfRadTolerance*secRMin ; 769 tolIRMin = fRmin1 + halfRadTolerance*secRMin ; 770 tolIRMax = fRmax1 - halfRadTolerance*secRMin ; 771 tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)* 772 (fRmax1 + halfRadTolerance*secRMax) ; 774 773 } 775 774 else 776 775 { 777 tolORMin = fRmin2 - 0.5*kRadTolerance*secRMin ;778 tolIRMin = fRmin2 + 0.5*kRadTolerance*secRMin ;779 tolIRMax = fRmax2 - 0.5*kRadTolerance*secRMin ;780 tolORMax2 = (fRmax2 + 0.5*kRadTolerance*secRMax)*781 (fRmax2 + 0.5*kRadTolerance*secRMax) ;776 tolORMin = fRmin2 - halfRadTolerance*secRMin ; 777 tolIRMin = fRmin2 + halfRadTolerance*secRMin ; 778 tolIRMax = fRmax2 - halfRadTolerance*secRMin ; 779 tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)* 780 (fRmax2 + halfRadTolerance*secRMax) ; 782 781 } 783 782 if ( tolORMin > 0 ) … … 791 790 tolIRMin2 = 0.0 ; 792 791 } 793 if ( tolIRMax > 0 ) tolIRMax2 = tolIRMax*tolIRMax ;794 else tolIRMax2 = 0.0 ;792 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; } 793 else { tolIRMax2 = 0.0; } 795 794 796 if (tolIRMin2 <= rhoi2 && rhoi2 <= tolIRMax2) 797 { 798 if ( seg && rhoi2 ) 799 { 800 // Psi = angle made with central (average) phi of shape 801 802 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 803 804 if (cosPsi >= cosHDPhiIT) return s ; 805 } 806 else return s ; 807 } 808 /* 809 else if (tolORMin2 <= rhoi2 && rhoi2 <= tolORMax2) 810 { 811 if ( seg && rhoi2 ) 812 { 813 // Psi = angle made with central (average) phi of shape 814 815 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 816 817 if (cosPsi >= cosHDPhiIT) return s ; 818 } 819 else return s ; 820 } 821 */ 795 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) ) 796 { 797 if ( !fPhiFullCone && rhoi2 ) 798 { 799 // Psi = angle made with central (average) phi of shape 800 801 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 802 803 if (cosPsi >= cosHDPhiIT) { return s; } 804 } 805 else 806 { 807 return s; 808 } 809 } 822 810 } 823 811 else // On/outside extent, and heading away -> cannot intersect … … 861 849 if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots 862 850 { 863 b = nt2/nt1;864 c = nt3/nt1;865 d = b*b-c;866 if ( nt3 > rout*kRadTolerance*secRMax || rout < 0)851 b = nt2/nt1; 852 c = nt3/nt1; 853 d = b*b-c ; 854 if ( (nt3 > rout*kRadTolerance*secRMax) || (rout < 0) ) 867 855 { 868 856 // If outside real cone (should be rho-rout>kRadTolerance*0.5 869 857 // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy 870 858 871 872 859 if (d >= 0) 873 860 { 874 861 875 if ( rout < 0 && nt3 <= 0)862 if ((rout < 0) && (nt3 <= 0)) 876 863 { 877 864 // Inside `shadow cone' with -ve radius … … 882 869 else 883 870 { 884 if ( b <= 0 && c >= 0) // both >=0, try smaller root871 if ((b <= 0) && (c >= 0)) // both >=0, try smaller root 885 872 { 886 873 s = -b - std::sqrt(d) ; … … 906 893 // Z ok. Check phi intersection if reqd 907 894 908 if ( ! seg ) return s ;895 if ( fPhiFullCone ) { return s; } 909 896 else 910 897 { … … 914 901 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 915 902 916 if ( cosPsi >= cosHDPhiIT ) return s ;903 if ( cosPsi >= cosHDPhiIT ) { return s; } 917 904 } 918 905 } … … 925 912 // check not inside, and heading through G4Cons (-> 0 to in) 926 913 927 if ( t3 > (rin + kRadTolerance*0.5*secRMin)* 928 (rin + kRadTolerance*0.5*secRMin) && 929 nt2 < 0 && 930 d >= 0 && 931 // nt2 < -kCarTolerance*secRMax/2/fDz && 932 // t2 < std::sqrt(t3)*v.z()*tanRMax && 933 // d > kCarTolerance*secRMax*(rout-b*tanRMax*v.z())/nt1 && 934 std::fabs(p.z()) <= tolIDz ) 914 if ( ( t3 > (rin + halfRadTolerance*secRMin)* 915 (rin + halfRadTolerance*secRMin) ) 916 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) ) 935 917 { 936 918 // Inside cones, delta r -ve, inside z extent 937 919 938 if ( seg)920 if ( !fPhiFullCone ) 939 921 { 940 922 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ; 941 923 942 if (cosPsi >= cosHDPhiIT) return 0.0 ;943 } 944 else return 0.0 ;924 if (cosPsi >= cosHDPhiIT) { return 0.0; } 925 } 926 else { return 0.0; } 945 927 } 946 928 } … … 952 934 s = -0.5*nt3/nt2 ; 953 935 954 if ( s < 0 ) return kInfinity ;// travel away936 if ( s < 0 ) { return kInfinity; } // travel away 955 937 else // s >= 0, If 'forwards'. Check z intersection 956 938 { 957 939 zi = p.z() + s*v.z() ; 958 940 959 if ( std::fabs(zi) <= tolODz && nt2 < 0)941 if ((std::fabs(zi) <= tolODz) && (nt2 < 0)) 960 942 { 961 943 // Z ok. Check phi intersection if reqd 962 944 963 if ( ! seg ) return s ;945 if ( fPhiFullCone ) { return s; } 964 946 else 965 947 { … … 969 951 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 970 952 971 if (cosPsi >= cosHDPhiIT) return s ;953 if (cosPsi >= cosHDPhiIT) { return s; } 972 954 } 973 955 } … … 1015 997 if ( std::fabs(zi) <= tolODz ) 1016 998 { 1017 if ( seg)999 if ( !fPhiFullCone ) 1018 1000 { 1019 1001 xi = p.x() + s*v.x() ; … … 1022 1004 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 1023 1005 1024 if (cosPsi >= cosHDPhiIT) snxt = s ;1006 if (cosPsi >= cosHDPhiIT) { snxt = s; } 1025 1007 } 1026 else return s ;1008 else { return s; } 1027 1009 } 1028 1010 } … … 1046 1028 ri = rMinAv + zi*tanRMin ; 1047 1029 1048 if ( ri > =0 )1030 if ( ri > 0 ) 1049 1031 { 1050 if ( s >= 0 && std::fabs(zi) <= tolODz) // s > 01032 if ( (s >= 0) && (std::fabs(zi) <= tolODz) ) // s > 0 1051 1033 { 1052 if ( seg)1034 if ( !fPhiFullCone ) 1053 1035 { 1054 1036 xi = p.x() + s*v.x() ; … … 1056 1038 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 1057 1039 1058 if (cosPsi >= cosHDPhiOT) snxt = s ;1040 if (cosPsi >= cosHDPhiOT) { snxt = s; } 1059 1041 } 1060 else return s ;1042 else { return s; } 1061 1043 } 1062 1044 } … … 1067 1049 ri = rMinAv + zi*tanRMin ; 1068 1050 1069 if ( s >= 0 && ri >= 0 && std::fabs(zi) <= tolODz) // s>01051 if ( (s >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // s>0 1070 1052 { 1071 if ( seg)1053 if ( !fPhiFullCone ) 1072 1054 { 1073 1055 xi = p.x() + s*v.x() ; … … 1075 1057 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 1076 1058 1077 if (cosPsi >= cosHDPhiIT) snxt = s ;1059 if (cosPsi >= cosHDPhiIT) { snxt = s; } 1078 1060 } 1079 else return s ;1061 else { return s; } 1080 1062 } 1081 1063 } … … 1095 1077 // Inside inner real cone, heading outwards, inside z range 1096 1078 1097 if ( seg)1079 if ( !fPhiFullCone ) 1098 1080 { 1099 1081 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ; 1100 1082 1101 if (cosPsi >= cosHDPhiIT) return 0.0 ;1083 if (cosPsi >= cosHDPhiIT) { return 0.0; } 1102 1084 } 1103 else return 0.0 ;1085 else { return 0.0; } 1104 1086 } 1105 1087 else … … 1123 1105 zi = p.z() + s*v.z() ; 1124 1106 1125 if ( s >= 0 && std::fabs(zi) <= tolODz) // s>01107 if ( (s >= 0) && (std::fabs(zi) <= tolODz) ) // s>0 1126 1108 { 1127 if ( seg)1109 if ( !fPhiFullCone ) 1128 1110 { 1129 1111 xi = p.x() + s*v.x() ; … … 1132 1114 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 1133 1115 1134 if ( cosPsi >= cosHDPhiIT ) snxt = s ;1116 if ( cosPsi >= cosHDPhiIT ) { snxt = s; } 1135 1117 } 1136 else return s ;1118 else { return s; } 1137 1119 } 1138 1120 } 1139 else return kInfinity ;1121 else { return kInfinity; } 1140 1122 } 1141 1123 } … … 1152 1134 zi = p.z() + s*v.z() ; 1153 1135 1154 if ( s >= 0 && std::fabs(zi) <= tolODz) // s>01136 if ( (s >= 0) && (std::fabs(zi) <= tolODz) ) // s>0 1155 1137 { 1156 if ( seg)1138 if ( !fPhiFullCone ) 1157 1139 { 1158 1140 xi = p.x() + s*v.x(); … … 1161 1143 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri; 1162 1144 1163 if (cosPsi >= cosHDPhiIT) snxt = s ;1145 if (cosPsi >= cosHDPhiIT) { snxt = s; } 1164 1146 } 1165 else return s ;1147 else { return s; } 1166 1148 } 1167 1149 } … … 1180 1162 // -> Should use some form of loop Construct 1181 1163 1182 if ( seg ) 1183 { 1184 // First phi surface (`S'tarting phi) 1185 1186 sinSPhi = std::sin(fSPhi) ; 1187 cosSPhi = std::cos(fSPhi) ; 1164 if ( !fPhiFullCone ) 1165 { 1166 // First phi surface (starting phi) 1167 1188 1168 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1189 1169 … … 1192 1172 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; 1193 1173 1194 if (Dist < kCarTolerance*0.5)1174 if (Dist < halfCarTolerance) 1195 1175 { 1196 1176 s = Dist/Comp ; … … 1198 1178 if ( s < snxt ) 1199 1179 { 1200 if ( s < 0 ) s = 0.0 ;1180 if ( s < 0 ) { s = 0.0; } 1201 1181 1202 1182 zi = p.z() + s*v.z() ; … … 1210 1190 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ; 1211 1191 1212 if ( rhoi2 >= tolORMin2 && rhoi2 <= tolORMax2)1192 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) ) 1213 1193 { 1214 1194 // z and r intersections good - check intersecting with 1215 1195 // correct half-plane 1216 1196 1217 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) snxt = s ;1218 } 1197 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = s; } 1198 } 1219 1199 } 1220 1200 } 1221 1201 } 1222 } 1223 // Second phi surface (`E'nding phi) 1224 1225 ePhi = fSPhi + fDPhi ; 1226 sinEPhi = std::sin(ePhi) ; 1227 cosEPhi = std::cos(ePhi) ; 1202 } 1203 1204 // Second phi surface (Ending phi) 1205 1228 1206 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 1229 1207 … … 1231 1209 { 1232 1210 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ; 1233 if (Dist < kCarTolerance*0.5)1211 if (Dist < halfCarTolerance) 1234 1212 { 1235 1213 s = Dist/Comp ; … … 1237 1215 if ( s < snxt ) 1238 1216 { 1239 if ( s < 0 ) s = 0.0 ;1217 if ( s < 0 ) { s = 0.0; } 1240 1218 1241 1219 zi = p.z() + s*v.z() ; … … 1249 1227 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ; 1250 1228 1251 if ( rhoi2 >= tolORMin2 && rhoi2 <= tolORMax2)1229 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) ) 1252 1230 { 1253 1231 // z and r intersections good - check intersecting with 1254 1232 // correct half-plane 1255 1233 1256 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) snxt = s ;1257 } 1234 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = s; } 1235 } 1258 1236 } 1259 1237 } … … 1261 1239 } 1262 1240 } 1263 if (snxt < kCarTolerance*0.5) snxt = 0.; 1264 1265 #ifdef consdebug 1266 G4cout.precision(24); 1267 G4cout<<"G4Cons::DistanceToIn(p,v) "<<G4endl; 1268 G4cout<<"position = "<<p<<G4endl; 1269 G4cout<<"direction = "<<v<<G4endl; 1270 G4cout<<"distance = "<<snxt<<G4endl; 1271 #endif 1241 if (snxt < halfCarTolerance) { snxt = 0.; } 1272 1242 1273 1243 return snxt ; … … 1283 1253 G4double G4Cons::DistanceToIn(const G4ThreeVector& p) const 1284 1254 { 1285 G4double safe=0.0, rho, safeR1, safeR2, safeZ ;1255 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ; 1286 1256 G4double tanRMin, secRMin, pRMin ; 1287 1257 G4double tanRMax, secRMax, pRMax ; 1288 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi ;1289 G4double cosPsi ;1290 1258 1291 1259 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; … … 1304 1272 safeR2 = (rho - pRMax)/secRMax ; 1305 1273 1306 if ( safeR1 > safeR2) safe = safeR1 ;1307 else safe = safeR2 ;1274 if ( safeR1 > safeR2) { safe = safeR1; } 1275 else { safe = safeR2; } 1308 1276 } 1309 1277 else … … 1314 1282 safe = (rho - pRMax)/secRMax ; 1315 1283 } 1316 if ( safeZ > safe ) safe = safeZ ; 1317 1318 if ( fDPhi < twopi && rho ) 1319 { 1320 phiC = fSPhi + fDPhi*0.5 ; 1321 cosPhiC = std::cos(phiC) ; 1322 sinPhiC = std::sin(phiC) ; 1323 1284 if ( safeZ > safe ) { safe = safeZ; } 1285 1286 if ( !fPhiFullCone && rho ) 1287 { 1324 1288 // Psi=angle from central phi to point 1325 1289 1326 cosPsi = (p.x()*cos PhiC + p.y()*sinPhiC)/rho ;1290 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ; 1327 1291 1328 1292 if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range 1329 1293 { 1330 if ( (p.y()*cos PhiC - p.x()*sinPhiC) <= 0.0 )1331 { 1332 safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi));1294 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 ) 1295 { 1296 safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi)); 1333 1297 } 1334 1298 else 1335 1299 { 1336 ePhi = fSPhi + fDPhi ; 1337 safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ; 1338 } 1339 if ( safePhi > safe ) safe = safePhi ; 1340 } 1341 } 1342 if ( safe < 0.0 ) safe = 0.0 ; 1300 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi); 1301 } 1302 if ( safePhi > safe ) { safe = safePhi; } 1303 } 1304 } 1305 if ( safe < 0.0 ) { safe = 0.0; } 1343 1306 1344 1307 return safe ; … … 1347 1310 /////////////////////////////////////////////////////////////// 1348 1311 // 1349 // Calculate distance to surface of shape from `inside', allowing for tolerance1312 // Calculate distance to surface of shape from 'inside', allowing for tolerance 1350 1313 // - Only Calc rmax intersection if no valid rmin intersection 1351 1314 1352 1315 G4double G4Cons::DistanceToOut( const G4ThreeVector& p, 1353 const G4ThreeVector& v,1354 const G4bool calcNorm,1355 G4bool *validNorm,1356 G4ThreeVector *n) const1316 const G4ThreeVector& v, 1317 const G4bool calcNorm, 1318 G4bool *validNorm, 1319 G4ThreeVector *n) const 1357 1320 { 1358 1321 ESide side = kNull, sider = kNull, sidephi = kNull; 1359 1322 1323 static const G4double halfCarTolerance=kCarTolerance*0.5; 1324 static const G4double halfRadTolerance=kRadTolerance*0.5; 1325 static const G4double halfAngTolerance=kAngTolerance*0.5; 1326 1360 1327 G4double snxt,sr,sphi,pdist ; 1361 1328 … … 1373 1340 // Vars for phi intersection: 1374 1341 1375 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;1376 G4double cPhi, sinCPhi, cosCPhi ;1377 1342 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ; 1378 1343 G4double zi, ri, deltaRoi2 ; … … 1384 1349 pdist = fDz - p.z() ; 1385 1350 1386 if (pdist > kCarTolerance*0.5)1351 if (pdist > halfCarTolerance) 1387 1352 { 1388 1353 snxt = pdist/v.z() ; … … 1396 1361 *validNorm = true ; 1397 1362 } 1398 return snxt = 0.0;1363 return snxt = 0.0; 1399 1364 } 1400 1365 } … … 1403 1368 pdist = fDz + p.z() ; 1404 1369 1405 if ( pdist > kCarTolerance*0.5)1370 if ( pdist > halfCarTolerance) 1406 1371 { 1407 1372 snxt = -pdist/v.z() ; … … 1465 1430 - fRmax1*(fRmax1 + kRadTolerance*secRMax); 1466 1431 } 1467 else deltaRoi2 = 1.0 ; 1468 1469 if ( nt1 && deltaRoi2 > 0.0 ) 1432 else 1433 { 1434 deltaRoi2 = 1.0; 1435 } 1436 1437 if ( nt1 && (deltaRoi2 > 0.0) ) 1470 1438 { 1471 1439 // Equation quadratic => 2 roots : second root must be leaving … … 1478 1446 { 1479 1447 // Check if on outer cone & heading outwards 1480 // NOTE: Should use rho-rout>-kRad tolerance*0.51448 // NOTE: Should use rho-rout>-kRadTolerance*0.5 1481 1449 1482 if (nt3 > - kRadTolerance*0.5&& nt2 >= 0 )1450 if (nt3 > -halfRadTolerance && nt2 >= 0 ) 1483 1451 { 1484 1452 if (calcNorm) … … 1486 1454 risec = std::sqrt(t3)*secRMax ; 1487 1455 *validNorm = true ; 1488 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) 1456 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 1489 1457 } 1490 1458 return snxt=0 ; … … 1497 1465 ri = tanRMax*zi + rMaxAv ; 1498 1466 1499 if ( (ri >= 0) && (-kRadTolerance*0.5 <= sr) && 1500 ( sr <= kRadTolerance*0.5) ) 1467 if ((ri >= 0) && (-halfRadTolerance <= sr) && (sr <= halfRadTolerance)) 1501 1468 { 1502 1469 // An intersection within the tolerance … … 1506 1473 sidetol = kRMax ; 1507 1474 } 1508 if ( (ri < 0) || (sr < kRadTolerance*0.5) )1475 if ( (ri < 0) || (sr < halfRadTolerance) ) 1509 1476 { 1510 1477 // Safety: if both roots -ve ensure that sr cannot `win' … … 1515 1482 ri = tanRMax*zi + rMaxAv ; 1516 1483 1517 if (ri >= 0 && sr2 > kRadTolerance*0.5) sr = sr2 ; 1484 if ((ri >= 0) && (sr2 > halfRadTolerance)) 1485 { 1486 sr = sr2; 1487 } 1518 1488 else 1519 1489 { 1520 1490 sr = kInfinity ; 1521 1491 1522 if( (-kRadTolerance*0.5 <= sr2) 1523 && ( sr2 <= kRadTolerance*0.5) ) 1492 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) ) 1524 1493 { 1525 1494 // An intersection within the tolerance. … … 1540 1509 if ( calcNorm ) 1541 1510 { 1542 risec = std::sqrt(t3)*secRMax 1543 *validNorm = true 1544 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) 1511 risec = std::sqrt(t3)*secRMax; 1512 *validNorm = true; 1513 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 1545 1514 } 1546 1515 return snxt = 0.0 ; 1547 1516 } 1548 1517 } 1549 else if ( nt2 && deltaRoi2 > 0.0)1518 else if ( nt2 && (deltaRoi2 > 0.0) ) 1550 1519 { 1551 1520 // Linear case (only one intersection) => point outside outer cone … … 1553 1522 if ( calcNorm ) 1554 1523 { 1555 risec = std::sqrt(t3)*secRMax 1556 *validNorm = true 1557 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax) 1524 risec = std::sqrt(t3)*secRMax; 1525 *validNorm = true; 1526 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 1558 1527 } 1559 1528 return snxt = 0.0 ; … … 1569 1538 // Check possible intersection within tolerance 1570 1539 1571 if ( slentol <= kCarTolerance*0.5)1540 if ( slentol <= halfCarTolerance ) 1572 1541 { 1573 1542 // An intersection within the tolerance was found. … … 1580 1549 // Calculate a normal vector, as below 1581 1550 1582 xi = p.x() + slentol*v.x() 1583 yi = p.y() + slentol*v.y() 1584 risec = std::sqrt(xi*xi + yi*yi)*secRMax 1585 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) 1551 xi = p.x() + slentol*v.x(); 1552 yi = p.y() + slentol*v.y(); 1553 risec = std::sqrt(xi*xi + yi*yi)*secRMax; 1554 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax); 1586 1555 1587 1556 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly … … 1595 1564 } 1596 1565 else // On the surface, but not heading out so we ignore this intersection 1597 { // (as it is within tolerance).1566 { // (as it is within tolerance). 1598 1567 slentol = kInfinity ; 1599 1568 } … … 1621 1590 d = b*b - c ; 1622 1591 1623 if ( d >= 0.0 )1592 if ( d >= 0.0 ) 1624 1593 { 1625 1594 // NOTE: should be rho-rin<kRadTolerance*0.5, … … 1630 1599 if ( nt2 < 0.0 ) 1631 1600 { 1632 if (calcNorm) *validNorm = false ;1633 return snxt = 0.0 1601 if (calcNorm) { *validNorm = false; } 1602 return snxt = 0.0; 1634 1603 } 1635 1604 } … … 1640 1609 ri = tanRMin*zi + rMinAv ; 1641 1610 1642 if( (ri >= 0.0) && (-kRadTolerance*0.5 <= sr2) && 1643 ( sr2 <= kRadTolerance*0.5) ) 1611 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) ) 1644 1612 { 1645 1613 // An intersection within the tolerance … … 1649 1617 sidetol = kRMax ; 1650 1618 } 1651 if( (ri<0) || (sr2 < kRadTolerance*0.5) )1619 if( (ri<0) || (sr2 < halfRadTolerance) ) 1652 1620 { 1653 1621 sr3 = -b + std::sqrt(d) ; … … 1656 1624 // distancetoout 1657 1625 1658 if ( sr3 > kCarTolerance*0.5)1626 if ( sr3 > halfRadTolerance ) 1659 1627 { 1660 1628 if( sr3 < sr ) … … 1670 1638 } 1671 1639 } 1672 else if ( sr3 > - kCarTolerance*0.5)1640 else if ( sr3 > -halfRadTolerance ) 1673 1641 { 1674 1642 // Intersection in tolerance. Store to check if it's good … … 1678 1646 } 1679 1647 } 1680 else if ( sr2 < sr && sr2 > kCarTolerance*0.5)1648 else if ( (sr2 < sr) && (sr2 > halfCarTolerance) ) 1681 1649 { 1682 1650 sr = sr2 ; 1683 1651 sider = kRMin ; 1684 1652 } 1685 else if (sr2 > - kCarTolerance*0.5)1653 else if (sr2 > -halfCarTolerance) 1686 1654 { 1687 1655 // Intersection in tolerance. Store to check if it's good … … 1690 1658 sidetol = kRMin ; 1691 1659 } 1692 if( slentol <= kCarTolerance*0.5)1660 if( slentol <= halfCarTolerance ) 1693 1661 { 1694 1662 // An intersection within the tolerance was found. … … 1706 1674 if( Normal.dot(v) > 0 ) 1707 1675 { 1708 // We will leave the Cone immediatelly 1676 // We will leave the cone immediately 1677 1709 1678 if( calcNorm ) 1710 1679 { … … 1731 1700 // Phi Intersection 1732 1701 1733 if ( fDPhi < twopi ) 1734 { 1735 sinSPhi = std::sin(fSPhi) ; 1736 cosSPhi = std::cos(fSPhi) ; 1737 ePhi = fSPhi + fDPhi ; 1738 sinEPhi = std::sin(ePhi) ; 1739 cosEPhi = std::cos(ePhi) ; 1740 cPhi = fSPhi + fDPhi*0.5 ; 1741 sinCPhi = std::sin(cPhi) ; 1742 cosCPhi = std::cos(cPhi) ; 1702 if ( !fPhiFullCone ) 1703 { 1743 1704 // add angle calculation with correction 1744 // of the difference in domain of atan2 and Sphi 1745 vphi = std::atan2(v.y(),v.x()) ; 1746 1747 if ( vphi < fSPhi - kAngTolerance*0.5 ) vphi += twopi ; 1748 else if ( vphi > fSPhi + fDPhi + kAngTolerance*0.5 ) vphi -= twopi; 1705 // of the difference in domain of atan2 and Sphi 1706 1707 vphi = std::atan2(v.y(),v.x()) ; 1708 1709 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; } 1710 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; } 1711 1749 1712 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later) 1750 1713 { … … 1761 1724 sidephi = kNull ; 1762 1725 1763 if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance) 1764 && (pDistE <= 0.5*kCarTolerance) ) ) 1765 || ( (fDPhi > pi) && !((pDistS > 0.5*kCarTolerance) 1766 && (pDistE > 0.5*kCarTolerance) ) ) ) 1767 { 1768 // Inside both phi *full* planes 1769 if ( compS < 0 ) 1726 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance) 1727 && (pDistE <= halfCarTolerance) ) ) 1728 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance) 1729 && (pDistE > halfCarTolerance) ) ) ) 1730 { 1731 // Inside both phi *full* planes 1732 if ( compS < 0 ) 1733 { 1734 sphi = pDistS/compS ; 1735 if (sphi >= -halfCarTolerance) 1770 1736 { 1771 sphi = pDistS/compS ; 1772 if (sphi >= -0.5*kCarTolerance) 1737 xi = p.x() + sphi*v.x() ; 1738 yi = p.y() + sphi*v.y() ; 1739 1740 // Check intersecting with correct half-plane 1741 // (if not -> no intersect) 1742 // 1743 if ( (std::abs(xi)<=kCarTolerance) 1744 && (std::abs(yi)<=kCarTolerance) ) 1773 1745 { 1774 xi = p.x() + sphi*v.x() ; 1775 yi = p.y() + sphi*v.y() ; 1776 1777 // Check intersecting with correct half-plane 1778 // (if not -> no intersect) 1779 // 1780 if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)){ 1781 sidephi= kSPhi; 1782 if(((fSPhi-0.5*kAngTolerance)<=vphi)&&((ePhi+0.5*kAngTolerance)>=vphi)) 1783 { sphi = kInfinity; } 1784 1746 sidephi= kSPhi; 1747 if ( ( fSPhi-halfAngTolerance <= vphi ) 1748 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) ) 1749 { 1750 sphi = kInfinity; 1785 1751 } 1786 else 1787 if ((yi*cosCPhi-xi*sinCPhi)>=0) 1788 { 1789 sphi = kInfinity ; 1790 } 1791 else 1792 { 1793 sidephi = kSPhi ; 1794 if ( pDistS > -kCarTolerance*0.5 ) 1795 { 1796 sphi = 0.0 ; // Leave by sphi immediately 1797 } 1798 } 1752 } 1753 else 1754 if ( (yi*cosCPhi-xi*sinCPhi)>=0 ) 1755 { 1756 sphi = kInfinity ; 1799 1757 } 1800 1758 else 1801 1759 { 1802 sphi = kInfinity ; 1803 } 1760 sidephi = kSPhi ; 1761 if ( pDistS > -halfCarTolerance ) 1762 { 1763 sphi = 0.0 ; // Leave by sphi immediately 1764 } 1765 } 1804 1766 } 1805 1767 else … … 1807 1769 sphi = kInfinity ; 1808 1770 } 1809 1810 if ( compE < 0 ) 1771 } 1772 else 1773 { 1774 sphi = kInfinity ; 1775 } 1776 1777 if ( compE < 0 ) 1778 { 1779 sphi2 = pDistE/compE ; 1780 1781 // Only check further if < starting phi intersection 1782 // 1783 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) ) 1811 1784 { 1812 sphi2 = pDistE/compE ; 1813 1814 // Only check further if < starting phi intersection 1815 // 1816 if ( (sphi2 > -0.5*kCarTolerance) && (sphi2 < sphi) ) 1785 xi = p.x() + sphi2*v.x() ; 1786 yi = p.y() + sphi2*v.y() ; 1787 1788 // Check intersecting with correct half-plane 1789 1790 if ( (std::abs(xi)<=kCarTolerance) 1791 && (std::abs(yi)<=kCarTolerance) ) 1817 1792 { 1818 xi = p.x() + sphi2*v.x() ; 1819 yi = p.y() + sphi2*v.y() ; 1820 1821 // Check intersecting with correct half-plane 1822 if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)){ 1823 // Leaving via ending phi 1824 if(!(((fSPhi-0.5*kAngTolerance)<=vphi)&&((ePhi+0.5*kAngTolerance)>=vphi))){ 1825 sidephi = kEPhi ; 1826 if ( pDistE <= -kCarTolerance*0.5 ) sphi = sphi2 ; 1827 else sphi = 0.0 ; 1828 } 1829 } 1830 else // Check intersecting with correct half-plane 1831 if ( (yi*cosCPhi-xi*sinCPhi) >= 0) 1793 // Leaving via ending phi 1794 1795 if(!( (fSPhi-halfAngTolerance <= vphi) 1796 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) ) 1832 1797 { 1833 // Leaving via ending phi1834 1835 1798 sidephi = kEPhi ; 1836 if ( pDistE <= - kCarTolerance*0.5 ) sphi = sphi2 ;1837 else sphi = 0.0 ;1799 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; } 1800 else { sphi = 0.0; } 1838 1801 } 1839 1802 } 1803 else // Check intersecting with correct half-plane 1804 if ( yi*cosCPhi-xi*sinCPhi >= 0 ) 1805 { 1806 // Leaving via ending phi 1807 1808 sidephi = kEPhi ; 1809 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; } 1810 else { sphi = 0.0; } 1811 } 1840 1812 } 1841 1813 } 1842 else 1843 { 1844 sphi = kInfinity ; 1845 } 1846 1847 1814 } 1815 else 1816 { 1817 sphi = kInfinity ; 1818 } 1848 1819 } 1849 1820 else … … 1852 1823 // within phi of shape, Step limited by rmax, else Step =0 1853 1824 1854 // vphi = std::atan2(v.y(),v.x()) ; 1855 1856 // if ( fSPhi < vphi && vphi < fSPhi + fDPhi ) sphi = kInfinity ; 1857 1858 if ( ((fSPhi-0.5*kAngTolerance) <= vphi) && (vphi <=( fSPhi + fDPhi)+0.5*kAngTolerance) ) 1859 { 1860 sphi = kInfinity ; 1861 } 1825 if ( (fSPhi-halfAngTolerance <= vphi) 1826 && (vphi <= fSPhi+fDPhi+halfAngTolerance) ) 1827 { 1828 sphi = kInfinity ; 1829 } 1862 1830 else 1863 1831 { … … 1868 1836 if ( sphi < snxt ) // Order intersecttions 1869 1837 { 1870 1871 1838 snxt=sphi ; 1839 side=sidephi ; 1872 1840 } 1873 1841 } … … 1880 1848 { 1881 1849 switch(side) 1882 { 1883 case kRMax: 1884 // Note: returned vector not normalised 1885 // (divide by frmax for unit vector) 1850 { // Note: returned vector not normalised 1851 case kRMax: // (divide by frmax for unit vector) 1886 1852 xi = p.x() + snxt*v.x() ; 1887 1853 yi = p.y() + snxt*v.y() ; … … 1891 1857 break ; 1892 1858 case kRMin: 1893 *validNorm =false ; // Rmin is inconvex1859 *validNorm = false ; // Rmin is inconvex 1894 1860 break ; 1895 1861 case kSPhi: 1896 1862 if ( fDPhi <= pi ) 1897 1863 { 1898 *n = G4ThreeVector(s td::sin(fSPhi),-std::cos(fSPhi),0);1864 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0); 1899 1865 *validNorm = true ; 1900 1866 } 1901 else *validNorm = false ; 1867 else 1868 { 1869 *validNorm = false ; 1870 } 1902 1871 break ; 1903 1872 case kEPhi: 1904 1873 if ( fDPhi <= pi ) 1905 1874 { 1906 *n = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);1875 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0); 1907 1876 *validNorm = true ; 1908 1877 } 1909 else *validNorm = false ; 1878 else 1879 { 1880 *validNorm = false ; 1881 } 1910 1882 break ; 1911 1883 case kPZ: … … 1925 1897 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 1926 1898 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 1927 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm << " mm"1928 << G4endl << G4endl ;1899 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm 1900 << " mm" << G4endl << G4endl ; 1929 1901 if( p.x() != 0. || p.x() != 0.) 1930 1902 { 1931 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree << " degree"1932 << G4endl << G4endl ;1903 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree 1904 << " degree" << G4endl << G4endl ; 1933 1905 } 1934 1906 G4cout << "Direction:" << G4endl << G4endl ; … … 1943 1915 } 1944 1916 } 1945 if (snxt < kCarTolerance*0.5) snxt = 0.; 1946 #ifdef consdebug 1947 G4cout.precision(24); 1948 G4cout<<"G4Cons::DistanceToOut(p,v,...) "<<G4endl; 1949 G4cout<<"position = "<<p<<G4endl; 1950 G4cout<<"direction = "<<v<<G4endl; 1951 G4cout<<"distance = "<<snxt<<G4endl; 1952 #endif 1917 if (snxt < halfCarTolerance) { snxt = 0.; } 1918 1953 1919 return snxt ; 1954 1920 } … … 1960 1926 G4double G4Cons::DistanceToOut(const G4ThreeVector& p) const 1961 1927 { 1962 G4double safe=0.0,rho,safeR1,safeR2,safeZ ; 1963 G4double tanRMin,secRMin,pRMin ; 1964 G4double tanRMax,secRMax,pRMax ; 1965 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi ; 1928 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi; 1929 G4double tanRMin, secRMin, pRMin; 1930 G4double tanRMax, secRMax, pRMax; 1966 1931 1967 1932 #ifdef G4CSGDEBUG … … 1975 1940 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 1976 1941 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 1977 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm << " mm"1978 << G4endl << G4endl ;1979 if( p.x() != 0. || p.x() != 0.)1980 { 1981 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree << " degree"1982 << G4endl << G4endl ;1983 } 1984 G4Exception("G4Cons::DistanceToOut(p)", "Notification", JustWarning,1985 "Point p is outside !?" );1942 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm 1943 << " mm" << G4endl << G4endl ; 1944 if( (p.x() != 0.) || (p.x() != 0.) ) 1945 { 1946 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree 1947 << " degree" << G4endl << G4endl ; 1948 } 1949 G4Exception("G4Cons::DistanceToOut(p)", "Notification", 1950 JustWarning, "Point p is outside !?" ); 1986 1951 } 1987 1952 #endif … … 1997 1962 safeR1 = (rho - pRMin)/secRMin ; 1998 1963 } 1999 else safeR1 = kInfinity ; 1964 else 1965 { 1966 safeR1 = kInfinity ; 1967 } 2000 1968 2001 1969 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; … … 2004 1972 safeR2 = (pRMax - rho)/secRMax ; 2005 1973 2006 if (safeR1 < safeR2) safe = safeR1 ;2007 else safe = safeR2 ;2008 if (safeZ < safe) safe = safeZ ;1974 if (safeR1 < safeR2) { safe = safeR1; } 1975 else { safe = safeR2; } 1976 if (safeZ < safe) { safe = safeZ ; } 2009 1977 2010 1978 // Check if phi divided, Calc distances closest phi plane 2011 1979 2012 if ( fDPhi < twopi)1980 if (!fPhiFullCone) 2013 1981 { 2014 1982 // Above/below central phi of G4Cons? 2015 1983 2016 phiC = fSPhi + fDPhi*0.5 ; 2017 cosPhiC = std::cos(phiC) ; 2018 sinPhiC = std::sin(phiC) ; 2019 2020 if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 ) 2021 { 2022 safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ; 1984 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 ) 1985 { 1986 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ; 2023 1987 } 2024 1988 else 2025 1989 { 2026 ePhi = fSPhi + fDPhi;2027 safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;2028 }2029 if (safePhi < safe) safe = safePhi ;2030 }2031 if ( safe < 0 ) safe = 0 ; 2032 return safe ; 1990 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ; 1991 } 1992 if (safePhi < safe) { safe = safePhi; } 1993 } 1994 if ( safe < 0 ) { safe = 0; } 1995 1996 return safe ; 2033 1997 } 2034 1998 … … 2068 2032 meshAngle = fDPhi/(noCrossSections - 1) ; 2069 2033 2070 // G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ;2071 2072 2034 meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ; 2073 2035 meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ; … … 2076 2038 // on the x axis. Will give better extent calculations when not rotated. 2077 2039 2078 if ( fDPhi == twopi && fSPhi == 0.0)2040 if ( fPhiFullCone && (fSPhi == 0.0) ) 2079 2041 { 2080 2042 sAngle = -meshAngle*0.5 ; … … 2102 2064 rMaxY2 = meshRMax2*sinCrossAngle ; 2103 2065 2104 // G4double RMin = (fRmin2 <= fRmin1) ? fRmin2 : fRmin1 ;2105 2106 2066 rMinX1 = fRmin1*cosCrossAngle ; 2107 2067 rMinY1 = fRmin1*sinCrossAngle ; … … 2127 2087 "Error in allocation of vertices. Out of memory !"); 2128 2088 } 2089 2129 2090 return vertices ; 2130 2091 } … … 2192 2153 rRand2 = RandFlat::shoot(fRmin2,fRmax2); 2193 2154 2194 if (fSPhi == 0. && fDPhi == twopi){ Afive = 0.; }2155 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; } 2195 2156 chose = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive); 2196 2157 … … 2225 2186 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) 2226 2187 { 2227 return G4ThreeVector (rRand1*cosu, rRand1*sinu,-1*fDz);2188 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz); 2228 2189 } 2229 2190 else if( (chose >= Aone + Atwo + Athree) -
trunk/source/geometry/solids/CSG/src/G4Tubs.cc
r850 r921 25 25 // 26 26 // 27 // $Id: G4Tubs.cc,v 1. 68 2008/06/23 13:37:39gcosmo Exp $28 // GEANT4 tag $Name: HEAD$27 // $Id: G4Tubs.cc,v 1.74 2008/11/06 15:26:53 gcosmo Exp $ 28 // GEANT4 tag $Name: geant4-09-02-cand-01 $ 29 29 // 30 30 // … … 108 108 "Invalid Z half-length"); 109 109 } 110 if ( pRMin < pRMax && pRMin >= 0) // Check radii110 if ( (pRMin < pRMax) && (pRMin >= 0) ) // Check radii 111 111 { 112 112 fRMin = pRMin ; … … 121 121 "Invalid radii."); 122 122 } 123 if ( pDPhi >= twopi ) // Check angles 123 124 fPhiFullTube = true; 125 if ( pDPhi >= twopi-kAngTolerance*0.5 ) // Check angles 124 126 { 125 127 fDPhi=twopi; 128 fSPhi=0; 126 129 } 127 130 else 128 131 { 132 fPhiFullTube = false; 129 133 if ( pDPhi > 0 ) 130 134 { … … 136 140 << " Negative delta-Phi ! - " 137 141 << pDPhi << G4endl; 138 G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException, 139 "Invalid dphi."); 140 } 141 } 142 143 // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0 144 145 fSPhi = pSPhi; 146 147 if ( fSPhi < 0 ) 148 { 149 fSPhi = twopi - std::fmod(std::fabs(fSPhi),twopi) ; 150 } 151 else 152 { 153 fSPhi = std::fmod(fSPhi,twopi) ; 154 } 155 if (fSPhi + fDPhi > twopi ) 156 { 157 fSPhi -= twopi ; 158 } 142 G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", 143 FatalException, "Invalid dphi."); 144 } 145 146 // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0 147 148 if ( pSPhi < 0 ) 149 { 150 fSPhi = twopi - std::fmod(std::fabs(pSPhi),twopi); 151 } 152 else 153 { 154 fSPhi = std::fmod(pSPhi,twopi) ; 155 } 156 if ( fSPhi+fDPhi > twopi ) 157 { 158 fSPhi -= twopi ; 159 } 160 } 161 InitializeTrigonometry(); 159 162 } 160 163 … … 290 293 yoff2 = yMax - yoffset ; 291 294 292 if ( yoff1 >= 0 && yoff2 >= 0 ) // Y limits cross max/min x => no change293 { 295 if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x 296 { // => no change 294 297 pMin = xMin ; 295 298 pMax = xMax ; … … 317 320 xoff2 = xMax - xoffset ; 318 321 319 if ( xoff1 >= 0 && xoff2 >= 0 ) // X limits cross max/min y => no change320 { 322 if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y 323 { // => no change 321 324 pMin = yMin ; 322 325 pMax = yMax ; … … 363 366 noEntries = vertices->size() ; 364 367 noBetweenSections4 = noEntries - 4 ; 365 /* 366 G4cout << "vertices = " << noEntries << "\t" 367 << "v-4 = " << noBetweenSections4 << G4endl; 368 G4cout << G4endl; 369 for(i = 0 ; i < noEntries ; i++ ) 370 { 371 G4cout << i << "\t" << "v.x = " << ((*vertices)[i]).x() << "\t" 372 << "v.y = " << ((*vertices)[i]).y() << "\t" 373 << "v.z = " << ((*vertices)[i]).z() << "\t" << G4endl; 374 } 375 G4cout << G4endl; 376 G4cout << "ClipCrossSection" << G4endl; 377 */ 368 378 369 for (i = 0 ; i < noEntries ; i += 4 ) 379 370 { 380 // G4cout << "section = " << i << G4endl; 381 ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ; 382 } 383 // G4cout << "ClipBetweenSections" << G4endl; 371 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ; 372 } 384 373 for (i = 0 ; i < noBetweenSections4 ; i += 4 ) 385 374 { 386 // G4cout << "between sections = " << i << G4endl; 387 ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ; 388 } 389 if (pMin != kInfinity || pMax != -kInfinity ) 375 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ; 376 } 377 if ((pMin != kInfinity) || (pMax != -kInfinity) ) 390 378 { 391 379 existsAfterClip = true ; … … 426 414 G4double r2,pPhi,tolRMin,tolRMax; 427 415 EInside in = kOutside ; 428 429 if (std::fabs(p.z()) <= fDz - kCarTolerance*0.5) 416 static const G4double halfCarTolerance=kCarTolerance*0.5; 417 static const G4double halfRadTolerance=kRadTolerance*0.5; 418 static const G4double halfAngTolerance=kAngTolerance*0.5; 419 420 if (std::fabs(p.z()) <= fDz - halfCarTolerance) 430 421 { 431 422 r2 = p.x()*p.x() + p.y()*p.y() ; 432 423 433 if (fRMin) { tolRMin = fRMin + kRadTolerance*0.5; }424 if (fRMin) { tolRMin = fRMin + halfRadTolerance ; } 434 425 else { tolRMin = 0 ; } 435 426 436 tolRMax = fRMax - kRadTolerance*0.5;427 tolRMax = fRMax - halfRadTolerance ; 437 428 438 if ( r2 >= tolRMin*tolRMin && r2 <= tolRMax*tolRMax)439 { 440 if ( f DPhi == twopi)429 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax)) 430 { 431 if ( fPhiFullTube ) 441 432 { 442 433 in = kInside ; … … 447 438 // if not inside, try outer tolerant phi boundaries 448 439 449 pPhi = std::atan2(p.y(),p.x()) ; 450 if ((tolRMin==0)&&(p.x()==0)&&(p.y()==0)) 440 if ((tolRMin==0)&&(p.x()<=halfCarTolerance)&&(p.y()<=halfCarTolerance)) 451 441 { 452 442 in=kSurface; … … 454 444 else 455 445 { 456 if ( pPhi < -kAngTolerance*0.5 ) { pPhi += twopi; } // 0<=pPhi<2pi 446 pPhi = std::atan2(p.y(),p.x()) ; 447 if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi 457 448 458 449 if ( fSPhi >= 0 ) 459 450 { 460 if ( (std::abs(pPhi) < kAngTolerance*0.5)461 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )451 if ( (std::abs(pPhi) < halfAngTolerance) 452 && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) ) 462 453 { 463 454 pPhi += twopi ; // 0 <= pPhi < 2pi 464 455 } 465 if ( (pPhi >= fSPhi + kAngTolerance*0.5)466 && (pPhi <= fSPhi + fDPhi - kAngTolerance*0.5) )456 if ( (pPhi >= fSPhi + halfAngTolerance) 457 && (pPhi <= fSPhi + fDPhi - halfAngTolerance) ) 467 458 { 468 459 in = kInside ; 469 460 } 470 else if ( (pPhi >= fSPhi - kAngTolerance*0.5)471 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )461 else if ( (pPhi >= fSPhi - halfAngTolerance) 462 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) ) 472 463 { 473 464 in = kSurface ; … … 476 467 else // fSPhi < 0 477 468 { 478 if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)479 && (pPhi >= fSPhi + fDPhi + kAngTolerance*0.5) ) {;}480 else if ( (pPhi <= fSPhi + twopi + kAngTolerance*0.5)481 && (pPhi >= fSPhi + fDPhi - kAngTolerance*0.5) )469 if ( (pPhi <= fSPhi + twopi - halfAngTolerance) 470 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} //kOutside 471 else if ( (pPhi <= fSPhi + twopi + halfAngTolerance) 472 && (pPhi >= fSPhi + fDPhi - halfAngTolerance) ) 482 473 { 483 474 in = kSurface ; … … 493 484 else // Try generous boundaries 494 485 { 495 tolRMin = fRMin - kRadTolerance*0.5;496 tolRMax = fRMax + kRadTolerance*0.5;486 tolRMin = fRMin - halfRadTolerance ; 487 tolRMax = fRMax + halfRadTolerance ; 497 488 498 489 if ( tolRMin < 0 ) { tolRMin = 0; } … … 500 491 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) ) 501 492 { 502 if ( fDPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis503 { 493 if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) ) 494 { // Continuous in phi or on z-axis 504 495 in = kSurface ; 505 496 } … … 508 499 pPhi = std::atan2(p.y(),p.x()) ; 509 500 510 if ( pPhi < - kAngTolerance*0.5) { pPhi += twopi; } // 0<=pPhi<2pi501 if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi 511 502 if ( fSPhi >= 0 ) 512 503 { 513 if ( (std::abs(pPhi) < kAngTolerance*0.5)514 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )504 if ( (std::abs(pPhi) < halfAngTolerance) 505 && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) ) 515 506 { 516 507 pPhi += twopi ; // 0 <= pPhi < 2pi 517 508 } 518 if ( (pPhi >= fSPhi - kAngTolerance*0.5)519 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )509 if ( (pPhi >= fSPhi - halfAngTolerance) 510 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) ) 520 511 { 521 512 in = kSurface ; … … 524 515 else // fSPhi < 0 525 516 { 526 if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)527 && (pPhi >= fSPhi + fDPhi + kAngTolerance*0.5) ) {;}517 if ( (pPhi <= fSPhi + twopi - halfAngTolerance) 518 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} // kOutside 528 519 else 529 520 { … … 535 526 } 536 527 } 537 else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5)528 else if (std::fabs(p.z()) <= fDz + halfCarTolerance) 538 529 { // Check within tolerant r limits 539 530 r2 = p.x()*p.x() + p.y()*p.y() ; 540 tolRMin = fRMin - kRadTolerance*0.5;541 tolRMax = fRMax + kRadTolerance*0.5;531 tolRMin = fRMin - halfRadTolerance ; 532 tolRMax = fRMax + halfRadTolerance ; 542 533 543 534 if ( tolRMin < 0 ) { tolRMin = 0; } … … 545 536 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) ) 546 537 { 547 if (f DPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis548 { 538 if (fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance)) 539 { // Continuous in phi or on z-axis 549 540 in = kSurface ; 550 541 } … … 553 544 pPhi = std::atan2(p.y(),p.x()) ; 554 545 555 if ( pPhi < - kAngTolerance*0.5) { pPhi += twopi; } // 0<=pPhi<2pi546 if ( pPhi < -halfAngTolerance ) { pPhi += twopi; } // 0<=pPhi<2pi 556 547 if ( fSPhi >= 0 ) 557 548 { 558 if ( (std::abs(pPhi) < kAngTolerance*0.5)559 && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )549 if ( (std::abs(pPhi) < halfAngTolerance) 550 && (std::abs(fSPhi + fDPhi - twopi) < halfAngTolerance) ) 560 551 { 561 552 pPhi += twopi ; // 0 <= pPhi < 2pi 562 553 } 563 if ( (pPhi >= fSPhi - kAngTolerance*0.5)564 && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )554 if ( (pPhi >= fSPhi - halfAngTolerance) 555 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) ) 565 556 { 566 557 in = kSurface; … … 569 560 else // fSPhi < 0 570 561 { 571 if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)572 && (pPhi >= fSPhi + fDPhi + kAngTolerance*0.5) ) {;}562 if ( (pPhi <= fSPhi + twopi - halfAngTolerance) 563 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;} 573 564 else 574 565 { … … 589 580 590 581 G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p ) const 591 { G4int noSurfaces = 0; 582 { 583 G4int noSurfaces = 0; 592 584 G4double rho, pPhi; 593 G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;594 585 G4double distZ, distRMin, distRMax; 595 586 G4double distSPhi = kInfinity, distEPhi = kInfinity; 587 588 static const G4double halfCarTolerance = 0.5*kCarTolerance; 589 static const G4double halfAngTolerance = 0.5*kAngTolerance; 590 596 591 G4ThreeVector norm, sumnorm(0.,0.,0.); 597 592 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0); … … 604 599 distZ = std::fabs(std::fabs(p.z()) - fDz); 605 600 606 if ( fDPhi < twopi) // && rho )// Protected against (0,0,z)607 { 608 if ( rho )601 if (!fPhiFullTube) // Protected against (0,0,z) 602 { 603 if ( rho > halfCarTolerance ) 609 604 { 610 605 pPhi = std::atan2(p.y(),p.x()); 611 606 612 if(pPhi < fSPhi- delta) { pPhi += twopi; }613 else if(pPhi > fSPhi+fDPhi+ delta) { pPhi -= twopi; }614 615 distSPhi = std::fabs( pPhi - fSPhi);607 if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; } 608 else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; } 609 610 distSPhi = std::fabs(pPhi - fSPhi); 616 611 distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 617 612 } … … 624 619 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 625 620 } 626 if ( rho > delta ){ nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }627 628 if( distRMax <= delta)621 if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); } 622 623 if( distRMax <= halfCarTolerance ) 629 624 { 630 625 noSurfaces ++; 631 626 sumnorm += nR; 632 627 } 633 if( fRMin && distRMin <= delta)628 if( fRMin && (distRMin <= halfCarTolerance) ) 634 629 { 635 630 noSurfaces ++; … … 638 633 if( fDPhi < twopi ) 639 634 { 640 if (distSPhi <= dAngle) // delta)635 if (distSPhi <= halfAngTolerance) 641 636 { 642 637 noSurfaces ++; 643 638 sumnorm += nPs; 644 639 } 645 if (distEPhi <= dAngle) // delta)640 if (distEPhi <= halfAngTolerance) 646 641 { 647 642 noSurfaces ++; … … 649 644 } 650 645 } 651 if (distZ <= delta)646 if (distZ <= halfCarTolerance) 652 647 { 653 648 noSurfaces ++; … … 668 663 else if ( noSurfaces == 1 ) { norm = sumnorm; } 669 664 else { norm = sumnorm.unit(); } 665 670 666 return norm; 671 667 } … … 714 710 } 715 711 } 716 if ( fDPhi < twopi&& rho ) // Protected against (0,0,z)712 if (!fPhiFullTube && rho ) // Protected against (0,0,z) 717 713 { 718 714 phi = std::atan2(p.y(),p.x()) ; … … 749 745 case kNRMin : // Inner radius 750 746 { 751 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho,0) ;747 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ; 752 748 break ; 753 749 } 754 750 case kNRMax : // Outer radius 755 751 { 756 norm = G4ThreeVector(p.x()/rho, p.y()/rho,0) ;752 norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ; 757 753 break ; 758 754 } 759 755 case kNZ : // + or - dz 760 756 { 761 if ( p.z() > 0 ) norm = G4ThreeVector(0,0,1) ;762 else norm = G4ThreeVector(0,0,-1) ;757 if ( p.z() > 0 ) { norm = G4ThreeVector(0,0,1) ; } 758 else { norm = G4ThreeVector(0,0,-1); } 763 759 break ; 764 760 } 765 761 case kNSPhi: 766 762 { 767 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi),0) ;763 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ; 768 764 break ; 769 765 } 770 766 case kNEPhi: 771 767 { 772 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi),0) ;768 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ; 773 769 break; 774 770 } … … 804 800 // 805 801 // NOTE: 806 // - Precalculations for phi trigonometry are Done `just in time' 807 // - `if valid' implies tolerant checking of intersection points 802 // - 'if valid' implies tolerant checking of intersection points 808 803 809 804 G4double G4Tubs::DistanceToIn( const G4ThreeVector& p, 810 805 const G4ThreeVector& v ) const 811 806 { 812 G4double snxt = kInfinity ; // snxt = default return value 813 814 // Precalculated trig for phi intersections - used by r,z intersections to 815 // check validity 816 817 G4bool seg ; // true if segmented 818 819 G4double hDPhi, hDPhiOT, hDPhiIT, cosHDPhiOT=0., cosHDPhiIT=0. ; 820 // half dphi + outer tolerance 821 822 G4double cPhi, sinCPhi=0., cosCPhi=0. ; // central phi 823 824 G4double tolORMin2, tolIRMax2 ; // `generous' radii squared 825 807 G4double snxt = kInfinity ; // snxt = default return value 808 G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared 826 809 G4double tolORMax2, tolIRMin2, tolODz, tolIDz ; 810 811 static const G4double halfCarTolerance = 0.5*kCarTolerance; 812 static const G4double halfRadTolerance = 0.5*kRadTolerance; 827 813 828 814 // Intersection point variables 829 815 // 830 G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi ; 831 832 G4double t1, t2, t3, b, c, d ; // Quadratic solver variables 833 834 G4double Comp ; 835 G4double cosSPhi, sinSPhi ; // Trig for phi start intersect 836 837 G4double ePhi, cosEPhi, sinEPhi ; // for phi end intersect 838 839 // Set phi divided flag and precalcs 840 841 if ( fDPhi < twopi ) 842 { 843 seg = true ; 844 hDPhi = 0.5*fDPhi ; // half delta phi 845 cPhi = fSPhi + hDPhi ; 846 hDPhiOT = hDPhi + 0.5*kAngTolerance ; // outers tol' half delta phi 847 hDPhiIT = hDPhi - 0.5*kAngTolerance ; 848 sinCPhi = std::sin(cPhi) ; 849 cosCPhi = std::cos(cPhi) ; 850 cosHDPhiOT = std::cos(hDPhiOT) ; 851 cosHDPhiIT = std::cos(hDPhiIT) ; 852 } 853 else 854 { 855 seg = false ; 856 } 857 816 G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ; 817 G4double t1, t2, t3, b, c, d ; // Quadratic solver variables 818 858 819 // Calculate tolerant rmin and rmax 859 820 860 821 if (fRMin > kRadTolerance) 861 822 { 862 tolORMin2 = (fRMin - 0.5*kRadTolerance)*(fRMin - 0.5*kRadTolerance) ;863 tolIRMin2 = (fRMin + 0.5*kRadTolerance)*(fRMin + 0.5*kRadTolerance) ;823 tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ; 824 tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ; 864 825 } 865 826 else … … 868 829 tolIRMin2 = 0.0 ; 869 830 } 870 tolORMax2 = (fRMax + 0.5*kRadTolerance)*(fRMax + 0.5*kRadTolerance) ;871 tolIRMax2 = (fRMax - 0.5*kRadTolerance)*(fRMax - 0.5*kRadTolerance) ;831 tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ; 832 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ; 872 833 873 834 // Intersection with Z surfaces 874 835 875 tolIDz = fDz - kCarTolerance*0.5;876 tolODz = fDz + kCarTolerance*0.5;836 tolIDz = fDz - halfCarTolerance ; 837 tolODz = fDz + halfCarTolerance ; 877 838 878 839 if (std::fabs(p.z()) >= tolIDz) … … 880 841 if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa 881 842 { 882 s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; 843 s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance 883 844 884 845 if(s < 0.0) { s = 0.0; } … … 890 851 // Check validity of intersection 891 852 892 if ( tolIRMin2 <= rho2 && rho2 <= tolIRMax2)893 { 894 if ( seg&& rho2)853 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2)) 854 { 855 if (!fPhiFullTube && rho2) 895 856 { 896 857 // Psi = angle made with central (average) phi of shape … … 899 860 iden = std::sqrt(rho2) ; 900 861 cosPsi = inum/iden ; 901 if (cosPsi >= cosHDPhiIT) return s ;862 if (cosPsi >= cosHDPhiIT) { return s ; } 902 863 } 903 864 else … … 909 870 else 910 871 { 911 if ( snxt< kCarTolerance*0.5) { snxt=0; }872 if ( snxt<halfCarTolerance ) { snxt=0; } 912 873 return snxt ; // On/outside extent, and heading away 913 874 // -> cannot intersect … … 934 895 b = t2/t1 ; 935 896 c = t3 - fRMax*fRMax ; 936 if ( t3 >= tolORMax2 && t2<0) // This also handles the tangent case897 if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case 937 898 { 938 899 // Try outer cylinder intersection … … 954 915 // Z ok. Check phi intersection if reqd 955 916 // 956 if ( !seg)917 if (fPhiFullTube) 957 918 { 958 919 return s ; … … 963 924 yi = p.y() + s*v.y() ; 964 925 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ; 965 if (cosPsi >= cosHDPhiIT) return s ;926 if (cosPsi >= cosHDPhiIT) { return s ; } 966 927 } 967 928 } // end if std::fabs(zi) … … 974 935 // check not inside, and heading through tubs (-> 0 to in) 975 936 976 if ( t3 > tolIRMin2 && t2 < 0 && std::fabs(p.z()) <= tolIDz)937 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz)) 977 938 { 978 939 // Inside both radii, delta r -ve, inside z extent 979 940 980 if ( seg)941 if (!fPhiFullTube) 981 942 { 982 943 inum = p.x()*cosCPhi + p.y()*sinCPhi ; … … 1004 965 snxt = c/(-b+std::sqrt(d)); // using safe solution 1005 966 // for quadratic equation 1006 if ( snxt <kCarTolerance*0.5) { snxt=0; }967 if ( snxt < halfCarTolerance ) { snxt=0; } 1007 968 return snxt ; 1008 969 } … … 1035 996 snxt= c/(-b+std::sqrt(d)); // using safe solution 1036 997 // for quadratic equation 1037 if ( snxt <kCarTolerance*0.5) { snxt=0; }998 if ( snxt < halfCarTolerance ) { snxt=0; } 1038 999 return snxt ; 1039 1000 } … … 1043 1004 } 1044 1005 } 1045 } // end if ( seg)1006 } // end if (!fPhiFullTube) 1046 1007 } // end if (t3>tolIRMin2) 1047 1008 } // end if (Inside Outer Radius) … … 1056 1017 1057 1018 s = -b + std::sqrt(d) ; 1058 if (s >= - 0.5*kCarTolerance) // check forwards1019 if (s >= -halfCarTolerance) // check forwards 1059 1020 { 1060 1021 // Check z intersection … … 1066 1027 // Z ok. Check phi 1067 1028 // 1068 if ( !seg)1029 if ( fPhiFullTube ) 1069 1030 { 1070 1031 return s ; … … 1098 1059 // -> use some form of loop Construct ? 1099 1060 // 1100 if ( seg)1061 if ( !fPhiFullTube ) 1101 1062 { 1102 1063 // First phi surface (Starting phi) 1103 1104 sinSPhi = std::sin(fSPhi) ; 1105 cosSPhi = std::cos(fSPhi) ; 1064 // 1106 1065 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 1107 1066 … … 1110 1069 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; 1111 1070 1112 if ( Dist < kCarTolerance*0.5)1071 if ( Dist < halfCarTolerance ) 1113 1072 { 1114 1073 s = Dist/Comp ; … … 1135 1094 // - check intersecting with correct half-plane 1136 1095 // 1137 if ((yi*cosCPhi-xi*sinCPhi) <= 0) snxt = s ;1138 } 1096 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = s; } 1097 } 1139 1098 } 1140 1099 } … … 1142 1101 } 1143 1102 1144 // Second phi surface (`E'nding phi) 1145 1146 ePhi = fSPhi + fDPhi ; 1147 sinEPhi = std::sin(ePhi) ; 1148 cosEPhi = std::cos(ePhi) ; 1103 // Second phi surface (Ending phi) 1104 1149 1105 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 1150 1106 … … 1153 1109 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ; 1154 1110 1155 if ( Dist < kCarTolerance*0.5)1111 if ( Dist < halfCarTolerance ) 1156 1112 { 1157 1113 s = Dist/Comp ; … … 1177 1133 // - check intersecting with correct half-plane 1178 1134 // 1179 if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) 1180 } 1135 if ( (yi*cosCPhi-xi*sinCPhi) >= 0 ) { snxt = s; } 1136 } //?? >=-halfCarTolerance 1181 1137 } 1182 1138 } 1183 1139 } 1184 1140 } // Comp < 0 1185 } // seg != 01186 if ( snxt< kCarTolerance*0.5) { snxt=0; }1141 } // !fPhiFullTube 1142 if ( snxt<halfCarTolerance ) { snxt=0; } 1187 1143 return snxt ; 1188 1144 } … … 1217 1173 { 1218 1174 G4double safe=0.0, rho, safe1, safe2, safe3 ; 1219 G4double phiC, cosPhiC, sinPhiC, safePhi,ePhi, cosPsi ;1175 G4double safePhi, cosPsi ; 1220 1176 1221 1177 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; … … 1228 1184 if ( safe3 > safe ) { safe = safe3; } 1229 1185 1230 if (fDPhi < twopi && rho) 1231 { 1232 phiC = fSPhi + fDPhi*0.5 ; 1233 cosPhiC = std::cos(phiC) ; 1234 sinPhiC = std::sin(phiC) ; 1235 1186 if ( (!fPhiFullTube) && (rho) ) 1187 { 1236 1188 // Psi=angle from central phi to point 1237 1189 // 1238 cosPsi = (p.x()*cos PhiC + p.y()*sinPhiC)/rho ;1239 1190 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ; 1191 1240 1192 if ( cosPsi < std::cos(fDPhi*0.5) ) 1241 1193 { 1242 1194 // Point lies outside phi range 1243 1195 1244 if ( (p.y()*cos PhiC - p.x()*sinPhiC) <= 0 )1245 { 1246 safePhi = std::fabs(p.x()*s td::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;1196 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 ) 1197 { 1198 safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ; 1247 1199 } 1248 1200 else 1249 1201 { 1250 ePhi = fSPhi + fDPhi ; 1251 safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ; 1202 safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ; 1252 1203 } 1253 1204 if ( safePhi > safe ) { safe = safePhi; } … … 1269 1220 G4ThreeVector *n ) const 1270 1221 { 1271 ESide side = kNull , sider = kNull, sidephi =kNull ;1272 G4double snxt, sr = kInfinity, sphi =kInfinity, pdist ;1222 ESide side=kNull , sider=kNull, sidephi=kNull ; 1223 G4double snxt, sr=kInfinity, sphi=kInfinity, pdist ; 1273 1224 G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ; 1274 1225 1226 static const G4double halfCarTolerance = kCarTolerance*0.5; 1227 static const G4double halfAngTolerance = kAngTolerance*0.5; 1228 1275 1229 // Vars for phi intersection: 1276 1230 1277 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;1278 G4double cPhi, sinCPhi, cosCPhi ;1279 1231 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ; 1280 1232 … … 1284 1236 { 1285 1237 pdist = fDz - p.z() ; 1286 if ( pdist > kCarTolerance*0.5)1238 if ( pdist > halfCarTolerance ) 1287 1239 { 1288 1240 snxt = pdist/v.z() ; … … 1303 1255 pdist = fDz + p.z() ; 1304 1256 1305 if ( pdist > kCarTolerance*0.5)1257 if ( pdist > halfCarTolerance ) 1306 1258 { 1307 1259 snxt = -pdist/v.z() ; … … 1326 1278 // Radial Intersections 1327 1279 // 1328 // Find inters ction with cylinders at rmax/rmin1280 // Find intersection with cylinders at rmax/rmin 1329 1281 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 1330 1282 // … … 1359 1311 b = t2/t1 ; 1360 1312 c = deltaR/t1 ; 1361 d2 = b*b-c;1362 if( d2>=0.){sr = -b + std::sqrt(d2);}1363 else {sr=0.;};1313 d2 = b*b-c; 1314 if( d2 >= 0 ) { sr = -b + std::sqrt(d2); } 1315 else { sr = 0.; } 1364 1316 sider = kRMax ; 1365 1317 } … … 1371 1323 if ( calcNorm ) 1372 1324 { 1373 // if ( p.x() || p.y() )1374 // {1375 // *n=G4ThreeVector(p.x(),p.y(),0);1376 // }1377 // else1378 // {1379 // *n=v;1380 // }1381 1325 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ; 1382 1326 *validNorm = true ; … … 1417 1361 c = deltaR/t1 ; 1418 1362 d2 = b*b-c; 1419 if( d2>=0.)1363 if( d2 >=0. ) 1420 1364 { 1421 1365 sr = -b + std::sqrt(d2) ; … … 1423 1367 } 1424 1368 else // Case: On the border+t2<kRadTolerance 1425 // (v is perpendicula ir to the surface)1369 // (v is perpendicular to the surface) 1426 1370 { 1427 1371 if (calcNorm) … … 1441 1385 c = deltaR/t1; 1442 1386 d2 = b*b-c; 1443 if( d2>=0.)1387 if( d2 >= 0 ) 1444 1388 { 1445 1389 sr = -b + std::sqrt(d2) ; … … 1447 1391 } 1448 1392 else // Case: On the border+t2<kRadTolerance 1449 // (v is perpendicula ir to the surface)1393 // (v is perpendicular to the surface) 1450 1394 { 1451 1395 if (calcNorm) … … 1461 1405 // Phi Intersection 1462 1406 1463 if ( fDPhi < twopi ) 1464 { 1465 sinSPhi = std::sin(fSPhi) ; 1466 cosSPhi = std::cos(fSPhi) ; 1467 ePhi = fSPhi + fDPhi ; 1468 sinEPhi = std::sin(ePhi) ; 1469 cosEPhi = std::cos(ePhi) ; 1470 cPhi = fSPhi + fDPhi*0.5 ; 1471 sinCPhi = std::sin(cPhi) ; 1472 cosCPhi = std::cos(cPhi) ; 1473 1407 if ( !fPhiFullTube ) 1408 { 1474 1409 // add angle calculation with correction 1475 1410 // of the difference in domain of atan2 and Sphi 1476 1411 // 1477 1412 vphi = std::atan2(v.y(),v.x()) ; 1478 1479 if ( vphi < fSPhi - kAngTolerance*0.5) { vphi += twopi; }1480 else if ( vphi > fSPhi + fDPhi + kAngTolerance*0.5) { vphi -= twopi; }1413 1414 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; } 1415 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; } 1481 1416 1482 1417 … … 1492 1427 compS = -sinSPhi*v.x() + cosSPhi*v.y() ; 1493 1428 compE = sinEPhi*v.x() - cosEPhi*v.y() ; 1429 1494 1430 sidephi = kNull; 1495 1431 1496 if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance)1497 && (pDistE <= 0.5*kCarTolerance) ) )1498 || ( (fDPhi > pi) && !((pDistS > 0.5*kCarTolerance)1499 && (pDistE > 0.5*kCarTolerance) ) ) )1432 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance) 1433 && (pDistE <= halfCarTolerance) ) ) 1434 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance) 1435 && (pDistE > halfCarTolerance) ) ) ) 1500 1436 { 1501 1437 // Inside both phi *full* planes … … 1505 1441 sphi = pDistS/compS ; 1506 1442 1507 if (sphi >= - 0.5*kCarTolerance)1443 if (sphi >= -halfCarTolerance) 1508 1444 { 1509 1445 xi = p.x() + sphi*v.x() ; … … 1513 1449 // (if not -> no intersect) 1514 1450 // 1515 if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance)) 1516 { sidephi = kSPhi; 1517 if (((fSPhi-0.5*kAngTolerance)<=vphi) 1518 &&((ePhi+0.5*kAngTolerance)>=vphi)) 1451 if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) ) 1452 { 1453 sidephi = kSPhi; 1454 if (((fSPhi-halfAngTolerance)<=vphi) 1455 &&((fSPhi+fDPhi+halfAngTolerance)>=vphi)) 1519 1456 { 1520 1457 sphi = kInfinity; 1521 1458 } 1522 1459 } 1523 else if ( (yi*cosCPhi-xi*sinCPhi)>=0)1460 else if ( yi*cosCPhi-xi*sinCPhi >=0 ) 1524 1461 { 1525 1462 sphi = kInfinity ; … … 1528 1465 { 1529 1466 sidephi = kSPhi ; 1530 if ( pDistS > - kCarTolerance*0.5)1467 if ( pDistS > -halfCarTolerance ) 1531 1468 { 1532 1469 sphi = 0.0 ; // Leave by sphi immediately … … 1550 1487 // Only check further if < starting phi intersection 1551 1488 // 1552 if ( (sphi2 > - 0.5*kCarTolerance) && (sphi2 < sphi) )1489 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) ) 1553 1490 { 1554 1491 xi = p.x() + sphi2*v.x() ; … … 1559 1496 // Leaving via ending phi 1560 1497 // 1561 if( !(((fSPhi-0.5*kAngTolerance)<=vphi)1562 &&((ePhi+0.5*kAngTolerance)>=vphi)))1498 if( !((fSPhi-halfAngTolerance <= vphi) 1499 &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) ) 1563 1500 { 1564 1501 sidephi = kEPhi ; 1565 if ( pDistE <= - kCarTolerance*0.5) { sphi = sphi2 ; }1566 else { sphi = 0.0 ;}1502 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; } 1503 else { sphi = 0.0 ; } 1567 1504 } 1568 1505 } … … 1574 1511 // 1575 1512 sidephi = kEPhi ; 1576 if ( pDistE <= - kCarTolerance*0.5) { sphi = sphi2 ; }1577 else { sphi = 0.0 ;}1513 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; } 1514 else { sphi = 0.0 ; } 1578 1515 } 1579 1516 } … … 1589 1526 // On z axis + travel not || to z axis -> if phi of vector direction 1590 1527 // within phi of shape, Step limited by rmax, else Step =0 1591 1592 // vphi = std::atan2(v.y(),v.x()) ;//defined previosly1593 // G4cout<<"In axis vphi="<<vphi1594 // <<" Sphi="<<fSPhi<<" Ephi="<<ePhi<<G4endl;1595 // old if ( (fSPhi < vphi) && (vphi < fSPhi + fDPhi) )1596 // new : correction for if statement, must be '<='1597 1528 1598 if ( ( (fSPhi-0.5*kAngTolerance)<= vphi)1599 && (vphi <= ( ePhi+0.5*kAngTolerance) ))1529 if ( (fSPhi - halfAngTolerance <= vphi) 1530 && (vphi <= fSPhi + fDPhi + halfAngTolerance ) ) 1600 1531 { 1601 1532 sphi = kInfinity ; … … 1640 1571 if ( fDPhi <= pi ) 1641 1572 { 1642 *n = G4ThreeVector(s td::sin(fSPhi),-std::cos(fSPhi),0) ;1573 *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ; 1643 1574 *validNorm = true ; 1644 1575 } … … 1652 1583 if (fDPhi <= pi) 1653 1584 { 1654 *n = G4ThreeVector(-s td::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;1585 *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ; 1655 1586 *validNorm = true ; 1656 1587 } … … 1662 1593 1663 1594 case kPZ: 1664 *n =G4ThreeVector(0,0,1) ;1665 *validNorm =true ;1595 *n = G4ThreeVector(0,0,1) ; 1596 *validNorm = true ; 1666 1597 break ; 1667 1598 … … 1690 1621 } 1691 1622 } 1692 if ( snxt< kCarTolerance*0.5) { snxt=0 ; }1623 if ( snxt<halfCarTolerance ) { snxt=0 ; } 1693 1624 1694 1625 return snxt ; … … 1701 1632 G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const 1702 1633 { 1703 G4double safe=0.0, rho, safeR1, safeR2, safeZ ; 1704 G4double safePhi, phiC, cosPhiC, sinPhiC, ePhi ; 1634 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ; 1705 1635 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 1706 1636 … … 1738 1668 // Check if phi divided, Calc distances closest phi plane 1739 1669 // 1740 if ( fDPhi < twopi ) 1741 { 1742 // Above/below central phi of Tubs? 1743 1744 phiC = fSPhi + fDPhi*0.5 ; 1745 cosPhiC = std::cos(phiC) ; 1746 sinPhiC = std::sin(phiC) ; 1747 1748 if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 ) 1749 { 1750 safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ; 1670 if ( !fPhiFullTube ) 1671 { 1672 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 ) 1673 { 1674 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ; 1751 1675 } 1752 1676 else 1753 1677 { 1754 ePhi = fSPhi + fDPhi ; 1755 safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ; 1678 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ; 1756 1679 } 1757 1680 if (safePhi < safe) { safe = safePhi ; } … … 1806 1729 // on the x axis. Will give better extent calculations when not rotated. 1807 1730 1808 if (f DPhi == pi*2.0 && fSPhi == 0) { sAngle = -meshAngle*0.5 ; }1809 else 1731 if (fPhiFullTube && (fSPhi == 0) ) { sAngle = -meshAngle*0.5 ; } 1732 else { sAngle = fSPhi ; } 1810 1733 1811 1734 vertices = new G4ThreeVectorList(); … … 1975 1898 if (fRMin != 0) 1976 1899 { 1977 if (f DPhi >= twopi)1900 if (fPhiFullTube) 1978 1901 { 1979 1902 pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ; … … 1986 1909 else 1987 1910 { 1988 if (f DPhi >= twopi)1911 if (fPhiFullTube) 1989 1912 { 1990 1913 pNURBS = new G4NURBScylinder (fRMax,fDz) ;
Note: See TracChangeset
for help on using the changeset viewer.