64 #if !defined(G4GEOM_USE_UTUBS)
78 using namespace CLHEP;
89 :
G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
101 std::ostringstream message;
102 message <<
"Negative Z half-length (" << pDz <<
") in solid: " <<
GetName();
105 if ( (pRMin >= pRMax) || (pRMin < 0) )
107 std::ostringstream message;
108 message <<
"Invalid values for radii in solid: " <<
GetName()
110 <<
" pRMin = " << pRMin <<
", pRMax = " << pRMax;
125 :
G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
126 fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
127 sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
128 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
129 fPhiFullTube(false), halfCarTolerance(0.), halfRadTolerance(0.),
149 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
150 fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
151 fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
152 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.sinCPhi),
153 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiOT),
154 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
155 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube),
156 halfCarTolerance(rhs.halfCarTolerance),
157 halfRadTolerance(rhs.halfRadTolerance),
158 halfAngTolerance(rhs.halfAngTolerance)
171 if (
this == &rhs) {
return *
this; }
229 G4double diff1, diff2, maxDiff, newMin, newMax;
230 G4double xoff1, xoff2, yoff1, yoff2, delta;
233 xMin = xoffset -
fRMax;
234 xMax = xoffset +
fRMax;
256 yMin = yoffset -
fRMax;
257 yMax = yoffset +
fRMax;
279 zMin = zoffset -
fDz;
280 zMax = zoffset +
fDz;
305 yoff1 = yoffset - yMin;
306 yoff2 = yMax - yoffset;
308 if ( (yoff1 >= 0) && (yoff2 >= 0) )
318 delta = fRMax*fRMax - yoff1*yoff1;
319 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
320 delta = fRMax*fRMax - yoff2*yoff2;
321 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
322 maxDiff = (diff1 > diff2) ? diff1:diff2;
323 newMin = xoffset - maxDiff;
324 newMax = xoffset + maxDiff;
325 pMin = (newMin < xMin) ? xMin : newMin;
326 pMax = (newMax > xMax) ? xMax : newMax;
332 xoff1 = xoffset - xMin;
333 xoff2 = xMax - xoffset;
335 if ( (xoff1 >= 0) && (xoff2 >= 0) )
345 delta = fRMax*fRMax - xoff1*xoff1;
346 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
347 delta = fRMax*fRMax - xoff2*xoff2;
348 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
349 maxDiff = (diff1 > diff2) ? diff1 : diff2;
350 newMin = yoffset - maxDiff;
351 newMax = yoffset + maxDiff;
352 pMin = (newMin < yMin) ? yMin : newMin;
353 pMax = (newMax > yMax) ? yMax : newMax;
372 G4int i, noEntries, noBetweenSections4;
373 G4bool existsAfterClip =
false;
379 noEntries = vertices->size();
380 noBetweenSections4 = noEntries - 4;
382 for ( i = 0 ; i < noEntries ; i += 4 )
386 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
392 existsAfterClip =
true;
410 existsAfterClip =
true;
416 return existsAfterClip;
432 r2 = p.x()*p.x() + p.y()*p.y() ;
435 else { tolRMin = 0 ; }
439 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
457 pPhi = std::atan2(p.y(),p.x()) ;
500 if ( tolRMin < 0 ) { tolRMin = 0; }
502 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
504 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
510 pPhi = std::atan2(p.y(),p.x()) ;
541 r2 = p.x()*p.x() + p.y()*p.y() ;
545 if ( tolRMin < 0 ) { tolRMin = 0; }
547 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
549 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
555 pPhi = std::atan2(p.y(),p.x()) ;
594 G4int noSurfaces = 0;
603 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
605 distRMin = std::fabs(rho -
fRMin);
606 distRMax = std::fabs(rho -
fRMax);
607 distZ = std::fabs(std::fabs(p.z()) -
fDz);
613 pPhi = std::atan2(p.y(),p.x());
618 distSPhi = std::fabs(pPhi -
fSPhi);
657 if ( p.z() >= 0.) { sumnorm += nZ; }
658 else { sumnorm -= nZ; }
660 if ( noSurfaces == 0 )
663 G4Exception(
"G4Tubs::SurfaceNormal(p)",
"GeomSolids1002",
666 G4cout<<
"G4Tubs::SN ( "<<p.x()<<
", "<<p.y()<<
", "<<p.z()<<
" ); "
668 G4cout.precision(oldprc) ;
672 else if ( noSurfaces == 1 ) { norm = sumnorm; }
673 else { norm = sumnorm.unit(); }
688 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
690 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
692 distRMin = std::fabs(rho -
fRMin) ;
693 distRMax = std::fabs(rho -
fRMax) ;
694 distZ = std::fabs(std::fabs(p.z()) -
fDz) ;
696 if (distRMin < distRMax)
698 if ( distZ < distRMin )
711 if ( distZ < distRMax )
724 phi = std::atan2(p.y(),p.x()) ;
726 if ( phi < 0 ) { phi += twopi; }
730 distSPhi = std::fabs(phi - (
fSPhi + twopi))*rho ;
734 distSPhi = std::fabs(phi -
fSPhi)*rho ;
736 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
738 if (distSPhi < distEPhi)
740 if ( distSPhi < distMin )
747 if ( distEPhi < distMin )
786 "Undefined side for valid surface normal to solid.");
820 G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
825 G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
848 if (std::fabs(p.z()) >= tolIDz)
850 if ( p.z()*v.z() < 0 )
852 sd = (std::fabs(p.z()) -
fDz)/std::fabs(v.z()) ;
854 if(sd < 0.0) { sd = 0.0; }
856 xi = p.x() + sd*v.x() ;
857 yi = p.y() + sd*v.y() ;
858 rho2 = xi*xi + yi*yi ;
862 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
869 iden = std::sqrt(rho2) ;
881 if ( snxt<halfCarTolerance ) { snxt=0; }
898 t1 = 1.0 - v.z()*v.z() ;
899 t2 = p.x()*v.x() + p.y()*v.y() ;
900 t3 = p.x()*p.x() + p.y()*p.y() ;
906 if ((t3 >= tolORMax2) && (t2<0))
916 sd = c/(-b+std::sqrt(d));
921 G4double fTerm = sd-std::fmod(sd,dRmax);
926 zi = p.z() + sd*v.z() ;
927 if (std::fabs(zi)<=tolODz)
937 xi = p.x() + sd*v.x() ;
938 yi = p.y() + sd*v.y() ;
951 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
958 iden = std::sqrt(t3) ;
979 snxt = c/(-b+std::sqrt(d));
981 if ( snxt < halfCarTolerance ) { snxt=0; }
999 c = t3 - fRMax*
fRMax;
1010 snxt= c/(-b+std::sqrt(d));
1012 if ( snxt < halfCarTolerance ) { snxt=0; }
1032 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1033 if (sd >= -halfCarTolerance)
1037 if(sd < 0.0) { sd = 0.0; }
1040 G4double fTerm = sd-std::fmod(sd,dRmax);
1043 zi = p.z() + sd*v.z() ;
1044 if (std::fabs(zi) <= tolODz)
1054 xi = p.x() + sd*v.x() ;
1055 yi = p.y() + sd*v.y() ;
1090 if ( Dist < halfCarTolerance )
1096 if ( sd < 0 ) { sd = 0.0; }
1097 zi = p.z() + sd*v.z() ;
1098 if ( std::fabs(zi) <= tolODz )
1100 xi = p.x() + sd*v.x() ;
1101 yi = p.y() + sd*v.y() ;
1102 rho2 = xi*xi + yi*yi ;
1104 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1105 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1108 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1130 if ( Dist < halfCarTolerance )
1136 if ( sd < 0 ) { sd = 0; }
1137 zi = p.z() + sd*v.z() ;
1138 if ( std::fabs(zi) <= tolODz )
1140 xi = p.x() + sd*v.x() ;
1141 yi = p.y() + sd*v.y() ;
1142 rho2 = xi*xi + yi*yi ;
1143 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1144 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1147 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1161 if ( snxt<halfCarTolerance ) { snxt=0; }
1193 G4double safe=0.0, rho, safe1, safe2, safe3 ;
1196 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1197 safe1 =
fRMin - rho ;
1198 safe2 = rho -
fRMax ;
1199 safe3 = std::fabs(p.z()) -
fDz ;
1201 if ( safe1 > safe2 ) { safe = safe1; }
1202 else { safe = safe2; }
1203 if ( safe3 > safe ) { safe = safe3; }
1211 if ( cosPsi < std::cos(
fDPhi*0.5) )
1223 if ( safePhi > safe ) { safe = safePhi; }
1226 if ( safe < 0 ) { safe = 0; }
1243 G4double deltaR, t1, t2, t3, b, c,
d2, roMin2 ;
1247 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1253 pdist =
fDz - p.z() ;
1256 snxt = pdist/v.z() ;
1269 else if ( v.z() < 0 )
1271 pdist =
fDz + p.z() ;
1275 snxt = -pdist/v.z() ;
1305 t1 = 1.0 - v.z()*v.z() ;
1306 t2 = p.x()*v.x() + p.y()*v.y() ;
1307 t3 = p.x()*p.x() + p.y()*p.y() ;
1310 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }
1330 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1349 roMin2 = t3 - t2*t2/t1 ;
1365 srd = c/(-b+std::sqrt(d2));
1370 if ( calcNorm ) { *validNorm =
false; }
1381 srd = -b + std::sqrt(d2) ;
1405 srd = -b + std::sqrt(d2) ;
1428 vphi = std::atan2(v.y(),v.x()) ;
1434 if ( p.x() || p.y() )
1457 sphi = pDistS/compS ;
1461 xi = p.x() + sphi*v.x() ;
1462 yi = p.y() + sphi*v.y() ;
1501 sphi2 = pDistE/compE ;
1507 xi = p.x() + sphi2*v.x() ;
1508 yi = p.y() + sphi2*v.y() ;
1519 else { sphi = 0.0 ; }
1530 else { sphi = 0.0 ; }
1576 xi = p.x() + snxt*v.x() ;
1577 yi = p.y() + snxt*v.y() ;
1583 *validNorm = false ;
1594 *validNorm = false ;
1606 *validNorm = false ;
1623 std::ostringstream message;
1624 G4int oldprc = message.precision(16);
1625 message <<
"Undefined side for valid surface normal to solid."
1627 <<
"Position:" << G4endl << G4endl
1628 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1629 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1630 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1631 <<
"Direction:" << G4endl << G4endl
1632 <<
"v.x() = " << v.x() << G4endl
1633 <<
"v.y() = " << v.y() << G4endl
1634 <<
"v.z() = " << v.z() << G4endl << G4endl
1635 <<
"Proposed distance :" << G4endl << G4endl
1636 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
1637 message.precision(oldprc) ;
1638 G4Exception(
"G4Tubs::DistanceToOut(p,v,..)",
"GeomSolids1002",
1654 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1655 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1667 G4cout.precision(oldprc) ;
1668 G4Exception(
"G4Tubs::DistanceToOut(p)",
"GeomSolids1002",
1675 safeR1 = rho -
fRMin ;
1676 safeR2 =
fRMax - rho ;
1678 if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1679 else { safe = safeR2 ; }
1683 safe =
fRMax - rho ;
1685 safeZ =
fDz - std::fabs(p.z()) ;
1687 if ( safeZ < safe ) { safe = safeZ ; }
1701 if (safePhi < safe) { safe = safePhi ; }
1703 if ( safe < 0 ) { safe = 0 ; }
1724 G4double meshAngle, meshRMax, crossAngle,
1725 cosCrossAngle, sinCrossAngle, sAngle;
1726 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1727 G4int crossSection, noCrossSections;
1743 meshAngle =
fDPhi/(noCrossSections - 1) ;
1753 else { sAngle =
fSPhi ; }
1759 vertices->reserve(noCrossSections*4);
1760 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1764 crossAngle = sAngle + crossSection*meshAngle ;
1765 cosCrossAngle = std::cos(crossAngle) ;
1766 sinCrossAngle = std::sin(crossAngle) ;
1768 rMaxX = meshRMax*cosCrossAngle ;
1769 rMaxY = meshRMax*sinCrossAngle ;
1778 rMinX = meshRMin*cosCrossAngle ;
1779 rMinY = meshRMin*sinCrossAngle ;
1797 "Error in allocation of vertices. Out of memory !");
1817 return new G4Tubs(*
this);
1826 G4int oldprc = os.precision(16);
1827 os <<
"-----------------------------------------------------------\n"
1828 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1829 <<
" ===================================================\n"
1830 <<
" Solid type: G4Tubs\n"
1831 <<
" Parameters: \n"
1832 <<
" inner radius : " <<
fRMin/
mm <<
" mm \n"
1833 <<
" outer radius : " <<
fRMax/
mm <<
" mm \n"
1834 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1835 <<
" starting phi : " <<
fSPhi/
degree <<
" degrees \n"
1836 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1837 <<
"-----------------------------------------------------------\n";
1838 os.precision(oldprc);
1849 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1850 aOne, aTwo, aThr, aFou;
1859 cosphi = std::cos(phi);
1860 sinphi = std::sin(phi);
1864 if( (
fSPhi == 0) && (
fDPhi == twopi) ) { aFou = 0; }
1868 if( (chose >=0) && (chose < aOne) )
1870 xRand = fRMax*cosphi;
1871 yRand = fRMax*sinphi;
1875 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1877 xRand = fRMin*cosphi;
1878 yRand = fRMin*sinphi;
1882 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1884 xRand = rRand*cosphi;
1885 yRand = rRand*sinphi;
1889 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1891 xRand = rRand*cosphi;
1892 yRand = rRand*sinphi;
1896 else if( (chose >= aOne + aTwo + 2.*aThr)
1897 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1899 xRand = rRand*std::cos(
fSPhi);
1900 yRand = rRand*std::sin(
fSPhi);
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
ThreeVector shoot(const G4int Ap, const G4int Af)
static const G4double kInfinity
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4double halfCarTolerance
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
G4Tubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool IsXLimited() const
const G4int kMinMeshSections
G4ThreeVector GetPointOnSurface() const
virtual void AddSolid(const G4Box &)=0
virtual G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
std::ostream & StreamInfo(std::ostream &os) const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
G4GeometryType GetEntityType() const
G4GLOB_DLL std::ostream G4cout
G4Polyhedron * fpPolyhedron
G4Polyhedron * CreatePolyhedron() const
G4double GetRadialTolerance() const
virtual G4Polyhedron * GetPolyhedron() const
std::vector< G4ThreeVector > G4ThreeVectorList
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4double halfAngTolerance
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double GetMinXExtent() const
G4double GetMaxZExtent() const
EInside Inside(const G4ThreeVector &p) const
G4double halfRadTolerance
static const double degree
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4double GetMaxYExtent() const
G4double GetMaxExtent(const EAxis pAxis) const
void CheckPhiAngles(G4double sPhi, G4double dPhi)
G4CSGSolid & operator=(const G4CSGSolid &rhs)
G4bool IsZLimited() const
G4double GetAngularTolerance() const
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
G4double GetMinExtent(const EAxis pAxis) const
const G4double kMeshAngleDefault
static G4GeometryTolerance * GetInstance()
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4Tubs & operator=(const G4Tubs &rhs)
const G4int kMaxMeshSections
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const