38 #if !defined(G4GEOM_USE_UTUBS)
53 using namespace CLHEP;
64 :
G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
76 std::ostringstream message;
77 message <<
"Negative Z half-length (" << pDz <<
") in solid: " <<
GetName();
80 if ( (pRMin >= pRMax) || (pRMin < 0) )
82 std::ostringstream message;
83 message <<
"Invalid values for radii in solid: " <<
GetName()
85 <<
" pRMin = " << pRMin <<
", pRMax = " << pRMax;
100 :
G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
101 fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
102 sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
103 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
104 fPhiFullTube(false), halfCarTolerance(0.), halfRadTolerance(0.),
124 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
125 fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
126 fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
127 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.sinCPhi),
128 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiOT),
129 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
130 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube),
131 halfCarTolerance(rhs.halfCarTolerance),
132 halfRadTolerance(rhs.halfRadTolerance),
133 halfAngTolerance(rhs.halfAngTolerance)
145 if (
this == &rhs) {
return *
this; }
190 G4double diff1, diff2, maxDiff, newMin, newMax;
191 G4double xoff1, xoff2, yoff1, yoff2, delta;
194 xMin = xoffset -
fRMax;
195 xMax = xoffset +
fRMax;
217 yMin = yoffset -
fRMax;
218 yMax = yoffset +
fRMax;
240 zMin = zoffset -
fDz;
241 zMax = zoffset +
fDz;
266 yoff1 = yoffset - yMin;
267 yoff2 = yMax - yoffset;
269 if ( (yoff1 >= 0) && (yoff2 >= 0) )
279 delta = fRMax*fRMax - yoff1*yoff1;
280 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
281 delta = fRMax*fRMax - yoff2*yoff2;
282 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
283 maxDiff = (diff1 > diff2) ? diff1:diff2;
284 newMin = xoffset - maxDiff;
285 newMax = xoffset + maxDiff;
286 pMin = (newMin < xMin) ? xMin : newMin;
287 pMax = (newMax > xMax) ? xMax : newMax;
293 xoff1 = xoffset - xMin;
294 xoff2 = xMax - xoffset;
296 if ( (xoff1 >= 0) && (xoff2 >= 0) )
306 delta = fRMax*fRMax - xoff1*xoff1;
307 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
308 delta = fRMax*fRMax - xoff2*xoff2;
309 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
310 maxDiff = (diff1 > diff2) ? diff1 : diff2;
311 newMin = yoffset - maxDiff;
312 newMax = yoffset + maxDiff;
313 pMin = (newMin < yMin) ? yMin : newMin;
314 pMax = (newMax > yMax) ? yMax : newMax;
333 G4int i, noEntries, noBetweenSections4;
334 G4bool existsAfterClip =
false;
340 noEntries = vertices->size();
341 noBetweenSections4 = noEntries - 4;
343 for ( i = 0 ; i < noEntries ; i += 4 )
347 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
353 existsAfterClip =
true;
371 existsAfterClip =
true;
377 return existsAfterClip;
393 r2 = p.x()*p.x() + p.y()*p.y() ;
396 else { tolRMin = 0 ; }
400 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
418 pPhi = std::atan2(p.y(),p.x()) ;
461 if ( tolRMin < 0 ) { tolRMin = 0; }
463 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
465 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
471 pPhi = std::atan2(p.y(),p.x()) ;
502 r2 = p.x()*p.x() + p.y()*p.y() ;
506 if ( tolRMin < 0 ) { tolRMin = 0; }
508 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
510 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
516 pPhi = std::atan2(p.y(),p.x()) ;
555 G4int noSurfaces = 0;
564 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
566 distRMin = std::fabs(rho -
fRMin);
567 distRMax = std::fabs(rho -
fRMax);
568 distZ = std::fabs(std::fabs(p.z()) -
fDz);
574 pPhi = std::atan2(p.y(),p.x());
579 distSPhi = std::fabs(pPhi -
fSPhi);
618 if ( p.z() >= 0.) { sumnorm += nZ; }
619 else { sumnorm -= nZ; }
621 if ( noSurfaces == 0 )
624 G4Exception(
"G4Tubs::SurfaceNormal(p)",
"GeomSolids1002",
627 G4cout<<
"G4Tubs::SN ( "<<p.x()<<
", "<<p.y()<<
", "<<p.z()<<
" ); "
629 G4cout.precision(oldprc) ;
633 else if ( noSurfaces == 1 ) { norm = sumnorm; }
634 else { norm = sumnorm.unit(); }
649 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
651 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
653 distRMin = std::fabs(rho -
fRMin) ;
654 distRMax = std::fabs(rho -
fRMax) ;
655 distZ = std::fabs(std::fabs(p.z()) -
fDz) ;
657 if (distRMin < distRMax)
659 if ( distZ < distRMin )
672 if ( distZ < distRMax )
685 phi = std::atan2(p.y(),p.x()) ;
687 if ( phi < 0 ) { phi += twopi; }
691 distSPhi = std::fabs(phi - (
fSPhi + twopi))*rho ;
695 distSPhi = std::fabs(phi -
fSPhi)*rho ;
697 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
699 if (distSPhi < distEPhi)
701 if ( distSPhi < distMin )
708 if ( distEPhi < distMin )
747 "Undefined side for valid surface normal to solid.");
781 G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
786 G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
809 if (std::fabs(p.z()) >= tolIDz)
811 if ( p.z()*v.z() < 0 )
813 sd = (std::fabs(p.z()) -
fDz)/std::fabs(v.z()) ;
815 if(sd < 0.0) { sd = 0.0; }
817 xi = p.x() + sd*v.x() ;
818 yi = p.y() + sd*v.y() ;
819 rho2 = xi*xi + yi*yi ;
823 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
830 iden = std::sqrt(rho2) ;
842 if ( snxt<halfCarTolerance ) { snxt=0; }
859 t1 = 1.0 - v.z()*v.z() ;
860 t2 = p.x()*v.x() + p.y()*v.y() ;
861 t3 = p.x()*p.x() + p.y()*p.y() ;
867 if ((t3 >= tolORMax2) && (t2<0))
877 sd = c/(-b+std::sqrt(d));
882 G4double fTerm = sd-std::fmod(sd,dRmax);
887 zi = p.z() + sd*v.z() ;
888 if (std::fabs(zi)<=tolODz)
898 xi = p.x() + sd*v.x() ;
899 yi = p.y() + sd*v.y() ;
912 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
919 iden = std::sqrt(t3) ;
940 snxt = c/(-b+std::sqrt(d));
942 if ( snxt < halfCarTolerance ) { snxt=0; }
960 c = t3 - fRMax*
fRMax;
971 snxt= c/(-b+std::sqrt(d));
973 if ( snxt < halfCarTolerance ) { snxt=0; }
993 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
994 if (sd >= -halfCarTolerance)
998 if(sd < 0.0) { sd = 0.0; }
1001 G4double fTerm = sd-std::fmod(sd,dRmax);
1004 zi = p.z() + sd*v.z() ;
1005 if (std::fabs(zi) <= tolODz)
1015 xi = p.x() + sd*v.x() ;
1016 yi = p.y() + sd*v.y() ;
1051 if ( Dist < halfCarTolerance )
1057 if ( sd < 0 ) { sd = 0.0; }
1058 zi = p.z() + sd*v.z() ;
1059 if ( std::fabs(zi) <= tolODz )
1061 xi = p.x() + sd*v.x() ;
1062 yi = p.y() + sd*v.y() ;
1063 rho2 = xi*xi + yi*yi ;
1065 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1066 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1069 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1091 if ( Dist < halfCarTolerance )
1097 if ( sd < 0 ) { sd = 0; }
1098 zi = p.z() + sd*v.z() ;
1099 if ( std::fabs(zi) <= tolODz )
1101 xi = p.x() + sd*v.x() ;
1102 yi = p.y() + sd*v.y() ;
1103 rho2 = xi*xi + yi*yi ;
1104 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1105 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1108 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1122 if ( snxt<halfCarTolerance ) { snxt=0; }
1154 G4double safe=0.0, rho, safe1, safe2, safe3 ;
1157 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1158 safe1 =
fRMin - rho ;
1159 safe2 = rho -
fRMax ;
1160 safe3 = std::fabs(p.z()) -
fDz ;
1162 if ( safe1 > safe2 ) { safe = safe1; }
1163 else { safe = safe2; }
1164 if ( safe3 > safe ) { safe = safe3; }
1172 if ( cosPsi < std::cos(
fDPhi*0.5) )
1184 if ( safePhi > safe ) { safe = safePhi; }
1187 if ( safe < 0 ) { safe = 0; }
1204 G4double deltaR, t1, t2, t3, b, c,
d2, roMin2 ;
1208 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1214 pdist =
fDz - p.z() ;
1217 snxt = pdist/v.z() ;
1230 else if ( v.z() < 0 )
1232 pdist =
fDz + p.z() ;
1236 snxt = -pdist/v.z() ;
1266 t1 = 1.0 - v.z()*v.z() ;
1267 t2 = p.x()*v.x() + p.y()*v.y() ;
1268 t3 = p.x()*p.x() + p.y()*p.y() ;
1271 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }
1291 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1310 roMin2 = t3 - t2*t2/t1 ;
1326 srd = c/(-b+std::sqrt(d2));
1331 if ( calcNorm ) { *validNorm =
false; }
1342 srd = -b + std::sqrt(d2) ;
1366 srd = -b + std::sqrt(d2) ;
1389 vphi = std::atan2(v.y(),v.x()) ;
1395 if ( p.x() || p.y() )
1418 sphi = pDistS/compS ;
1422 xi = p.x() + sphi*v.x() ;
1423 yi = p.y() + sphi*v.y() ;
1462 sphi2 = pDistE/compE ;
1468 xi = p.x() + sphi2*v.x() ;
1469 yi = p.y() + sphi2*v.y() ;
1480 else { sphi = 0.0 ; }
1491 else { sphi = 0.0 ; }
1537 xi = p.x() + snxt*v.x() ;
1538 yi = p.y() + snxt*v.y() ;
1544 *validNorm = false ;
1555 *validNorm = false ;
1567 *validNorm = false ;
1584 std::ostringstream message;
1585 G4int oldprc = message.precision(16);
1586 message <<
"Undefined side for valid surface normal to solid."
1588 <<
"Position:" << G4endl << G4endl
1589 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1590 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1591 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1592 <<
"Direction:" << G4endl << G4endl
1593 <<
"v.x() = " << v.x() << G4endl
1594 <<
"v.y() = " << v.y() << G4endl
1595 <<
"v.z() = " << v.z() << G4endl << G4endl
1596 <<
"Proposed distance :" << G4endl << G4endl
1597 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
1598 message.precision(oldprc) ;
1599 G4Exception(
"G4Tubs::DistanceToOut(p,v,..)",
"GeomSolids1002",
1615 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1616 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1628 G4cout.precision(oldprc) ;
1629 G4Exception(
"G4Tubs::DistanceToOut(p)",
"GeomSolids1002",
1636 safeR1 = rho -
fRMin ;
1637 safeR2 =
fRMax - rho ;
1639 if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1640 else { safe = safeR2 ; }
1644 safe =
fRMax - rho ;
1646 safeZ =
fDz - std::fabs(p.z()) ;
1648 if ( safeZ < safe ) { safe = safeZ ; }
1662 if (safePhi < safe) { safe = safePhi ; }
1664 if ( safe < 0 ) { safe = 0 ; }
1685 G4double meshAngle, meshRMax, crossAngle,
1686 cosCrossAngle, sinCrossAngle, sAngle;
1687 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1688 G4int crossSection, noCrossSections;
1704 meshAngle =
fDPhi/(noCrossSections - 1) ;
1714 else { sAngle =
fSPhi ; }
1720 vertices->reserve(noCrossSections*4);
1721 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1725 crossAngle = sAngle + crossSection*meshAngle ;
1726 cosCrossAngle = std::cos(crossAngle) ;
1727 sinCrossAngle = std::sin(crossAngle) ;
1729 rMaxX = meshRMax*cosCrossAngle ;
1730 rMaxY = meshRMax*sinCrossAngle ;
1739 rMinX = meshRMin*cosCrossAngle ;
1740 rMinY = meshRMin*sinCrossAngle ;
1758 "Error in allocation of vertices. Out of memory !");
1787 G4int oldprc = os.precision(16);
1788 os <<
"-----------------------------------------------------------\n"
1789 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1790 <<
" ===================================================\n"
1791 <<
" Solid type: G4Tubs\n"
1792 <<
" Parameters: \n"
1793 <<
" inner radius : " <<
fRMin/
mm <<
" mm \n"
1794 <<
" outer radius : " <<
fRMax/
mm <<
" mm \n"
1795 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1796 <<
" starting phi : " <<
fSPhi/
degree <<
" degrees \n"
1797 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1798 <<
"-----------------------------------------------------------\n";
1799 os.precision(oldprc);
1810 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1811 aOne, aTwo, aThr, aFou;
1820 cosphi = std::cos(phi);
1821 sinphi = std::sin(phi);
1825 if( (
fSPhi == 0) && (
fDPhi == twopi) ) { aFou = 0; }
1829 if( (chose >=0) && (chose < aOne) )
1831 xRand = fRMax*cosphi;
1832 yRand = fRMax*sinphi;
1836 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1838 xRand = fRMin*cosphi;
1839 yRand = fRMin*sinphi;
1843 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1845 xRand = rRand*cosphi;
1846 yRand = rRand*sinphi;
1850 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1852 xRand = rRand*cosphi;
1853 yRand = rRand*sinphi;
1857 else if( (chose >= aOne + aTwo + 2.*aThr)
1858 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1860 xRand = rRand*std::cos(
fSPhi);
1861 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)
G4GeometryType GetEntityType() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
static const G4double kInfinity
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
G4ThreeVector GetPointOnSurface() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
void CheckPhiAngles(G4double sPhi, G4double dPhi)
G4bool IsXLimited() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
const G4int kMinMeshSections
virtual void AddSolid(const G4Box &)=0
G4double GetMaxXExtent() const
G4OTubs & operator=(const G4OTubs &rhs)
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
G4double halfRadTolerance
G4double GetRadialTolerance() const
G4OTubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)
std::vector< G4ThreeVector > G4ThreeVectorList
G4Polyhedron * CreatePolyhedron() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double GetMinXExtent() const
G4double GetMaxZExtent() const
G4double halfAngTolerance
G4double halfCarTolerance
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
static const double degree
G4double GetMaxYExtent() const
EInside Inside(const G4ThreeVector &p) const
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
virtual G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
G4double GetMaxExtent(const EAxis pAxis) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
G4bool IsZLimited() const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) 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()
std::ostream & StreamInfo(std::ostream &os) const
const G4int kMaxMeshSections