50 using namespace CLHEP;
61 :
G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
73 std::ostringstream message;
74 message <<
"Negative Z half-length (" << pDz <<
") in solid: " <<
GetName();
77 if ( (pRMin >= pRMax) || (pRMin < 0) )
79 std::ostringstream message;
80 message <<
"Invalid values for radii in solid: " <<
GetName()
82 <<
" pRMin = " << pRMin <<
", pRMax = " << pRMax;
97 :
G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
98 fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
99 sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
100 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
101 fPhiFullTube(false), halfCarTolerance(0.), halfRadTolerance(0.),
121 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
122 fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
123 fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
124 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
125 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
126 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
127 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube),
128 halfCarTolerance(rhs.halfCarTolerance),
129 halfRadTolerance(rhs.halfRadTolerance),
130 halfAngTolerance(rhs.halfAngTolerance)
142 if (
this == &rhs) {
return *
this; }
187 G4double diff1, diff2, maxDiff, newMin, newMax;
188 G4double xoff1, xoff2, yoff1, yoff2, delta;
191 xMin = xoffset -
fRMax;
192 xMax = xoffset +
fRMax;
214 yMin = yoffset -
fRMax;
215 yMax = yoffset +
fRMax;
237 zMin = zoffset -
fDz;
238 zMax = zoffset +
fDz;
263 yoff1 = yoffset - yMin;
264 yoff2 = yMax - yoffset;
266 if ( (yoff1 >= 0) && (yoff2 >= 0) )
276 delta = fRMax*fRMax - yoff1*yoff1;
277 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
278 delta = fRMax*fRMax - yoff2*yoff2;
279 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
280 maxDiff = (diff1 > diff2) ? diff1:diff2;
281 newMin = xoffset - maxDiff;
282 newMax = xoffset + maxDiff;
283 pMin = (newMin < xMin) ? xMin : newMin;
284 pMax = (newMax > xMax) ? xMax : newMax;
290 xoff1 = xoffset - xMin;
291 xoff2 = xMax - xoffset;
293 if ( (xoff1 >= 0) && (xoff2 >= 0) )
303 delta = fRMax*fRMax - xoff1*xoff1;
304 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
305 delta = fRMax*fRMax - xoff2*xoff2;
306 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
307 maxDiff = (diff1 > diff2) ? diff1 : diff2;
308 newMin = yoffset - maxDiff;
309 newMax = yoffset + maxDiff;
310 pMin = (newMin < yMin) ? yMin : newMin;
311 pMax = (newMax > yMax) ? yMax : newMax;
330 G4int i, noEntries, noBetweenSections4;
331 G4bool existsAfterClip =
false;
337 noEntries = vertices->size();
338 noBetweenSections4 = noEntries - 4;
340 for ( i = 0 ; i < noEntries ; i += 4 )
344 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
350 existsAfterClip =
true;
368 existsAfterClip =
true;
374 return existsAfterClip;
390 r2 = p.x()*p.x() + p.y()*p.y() ;
393 else { tolRMin = 0 ; }
397 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
415 pPhi = std::atan2(p.y(),p.x()) ;
458 if ( tolRMin < 0 ) { tolRMin = 0; }
460 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
462 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
468 pPhi = std::atan2(p.y(),p.x()) ;
499 r2 = p.x()*p.x() + p.y()*p.y() ;
503 if ( tolRMin < 0 ) { tolRMin = 0; }
505 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
507 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
513 pPhi = std::atan2(p.y(),p.x()) ;
552 G4int noSurfaces = 0;
561 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
563 distRMin = std::fabs(rho -
fRMin);
564 distRMax = std::fabs(rho -
fRMax);
565 distZ = std::fabs(std::fabs(p.z()) -
fDz);
571 pPhi = std::atan2(p.y(),p.x());
576 distSPhi = std::fabs(pPhi -
fSPhi);
615 if ( p.z() >= 0.) { sumnorm += nZ; }
616 else { sumnorm -= nZ; }
618 if ( noSurfaces == 0 )
621 G4Exception(
"G4Tubs::SurfaceNormal(p)",
"GeomSolids1002",
624 G4cout<<
"G4Tubs::SN ( "<<p.x()<<
", "<<p.y()<<
", "<<p.z()<<
" ); "
626 G4cout.precision(oldprc) ;
630 else if ( noSurfaces == 1 ) { norm = sumnorm; }
631 else { norm = sumnorm.unit(); }
646 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
648 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
650 distRMin = std::fabs(rho -
fRMin) ;
651 distRMax = std::fabs(rho -
fRMax) ;
652 distZ = std::fabs(std::fabs(p.z()) -
fDz) ;
654 if (distRMin < distRMax)
656 if ( distZ < distRMin )
669 if ( distZ < distRMax )
682 phi = std::atan2(p.y(),p.x()) ;
684 if ( phi < 0 ) { phi +=
twopi; }
688 distSPhi = std::fabs(phi - (
fSPhi +
twopi))*rho ;
692 distSPhi = std::fabs(phi -
fSPhi)*rho ;
694 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
696 if (distSPhi < distEPhi)
698 if ( distSPhi < distMin )
705 if ( distEPhi < distMin )
744 "Undefined side for valid surface normal to solid.");
778 G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
783 G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
806 if (std::fabs(p.z()) >= tolIDz)
808 if ( p.z()*v.z() < 0 )
810 sd = (std::fabs(p.z()) -
fDz)/std::fabs(v.z()) ;
812 if(sd < 0.0) { sd = 0.0; }
814 xi = p.x() + sd*v.x() ;
815 yi = p.y() + sd*v.y() ;
816 rho2 = xi*xi + yi*yi ;
820 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
827 iden = std::sqrt(rho2) ;
839 if ( snxt<halfCarTolerance ) { snxt=0; }
856 t1 = 1.0 - v.z()*v.z() ;
857 t2 = p.x()*v.x() + p.y()*v.y() ;
858 t3 = p.x()*p.x() + p.y()*p.y() ;
864 if ((t3 >= tolORMax2) && (t2<0))
874 sd = c/(-b+std::sqrt(d));
879 G4double fTerm = sd-std::fmod(sd,dRmax);
884 zi = p.z() + sd*v.z() ;
885 if (std::fabs(zi)<=tolODz)
895 xi = p.x() + sd*v.x() ;
896 yi = p.y() + sd*v.y() ;
909 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.z()) <= tolIDz))
916 iden = std::sqrt(t3) ;
937 snxt = c/(-b+std::sqrt(d));
939 if ( snxt < halfCarTolerance ) { snxt=0; }
957 c = t3 - fRMax*
fRMax;
968 snxt= c/(-b+std::sqrt(d));
970 if ( snxt < halfCarTolerance ) { snxt=0; }
990 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
991 if (sd >= -halfCarTolerance)
995 if(sd < 0.0) { sd = 0.0; }
998 G4double fTerm = sd-std::fmod(sd,dRmax);
1001 zi = p.z() + sd*v.z() ;
1002 if (std::fabs(zi) <= tolODz)
1012 xi = p.x() + sd*v.x() ;
1013 yi = p.y() + sd*v.y() ;
1048 if ( Dist < halfCarTolerance )
1054 if ( sd < 0 ) { sd = 0.0; }
1055 zi = p.z() + sd*v.z() ;
1056 if ( std::fabs(zi) <= tolODz )
1058 xi = p.x() + sd*v.x() ;
1059 yi = p.y() + sd*v.y() ;
1060 rho2 = xi*xi + yi*yi ;
1062 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1063 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1066 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1088 if ( Dist < halfCarTolerance )
1094 if ( sd < 0 ) { sd = 0; }
1095 zi = p.z() + sd*v.z() ;
1096 if ( std::fabs(zi) <= tolODz )
1098 xi = p.x() + sd*v.x() ;
1099 yi = p.y() + sd*v.y() ;
1100 rho2 = xi*xi + yi*yi ;
1101 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1102 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1105 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1119 if ( snxt<halfCarTolerance ) { snxt=0; }
1151 G4double safe=0.0, rho, safe1, safe2, safe3 ;
1154 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1155 safe1 =
fRMin - rho ;
1156 safe2 = rho -
fRMax ;
1157 safe3 = std::fabs(p.z()) -
fDz ;
1159 if ( safe1 > safe2 ) { safe = safe1; }
1160 else { safe = safe2; }
1161 if ( safe3 > safe ) { safe = safe3; }
1169 if ( cosPsi < std::cos(
fDPhi*0.5) )
1181 if ( safePhi > safe ) { safe = safePhi; }
1184 if ( safe < 0 ) { safe = 0; }
1201 G4double deltaR, t1, t2, t3, b, c,
d2, roMin2 ;
1205 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1211 pdist =
fDz - p.z() ;
1214 snxt = pdist/v.z() ;
1227 else if ( v.z() < 0 )
1229 pdist =
fDz + p.z() ;
1233 snxt = -pdist/v.z() ;
1263 t1 = 1.0 - v.z()*v.z() ;
1264 t2 = p.x()*v.x() + p.y()*v.y() ;
1265 t3 = p.x()*p.x() + p.y()*p.y() ;
1268 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }
1288 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1307 roMin2 = t3 - t2*t2/t1 ;
1323 srd = c/(-b+std::sqrt(d2));
1328 if ( calcNorm ) { *validNorm =
false; }
1339 srd = -b + std::sqrt(d2) ;
1363 srd = -b + std::sqrt(d2) ;
1386 vphi = std::atan2(v.y(),v.x()) ;
1392 if ( p.x() || p.y() )
1415 sphi = pDistS/compS ;
1419 xi = p.x() + sphi*v.x() ;
1420 yi = p.y() + sphi*v.y() ;
1459 sphi2 = pDistE/compE ;
1465 xi = p.x() + sphi2*v.x() ;
1466 yi = p.y() + sphi2*v.y() ;
1477 else { sphi = 0.0 ; }
1488 else { sphi = 0.0 ; }
1534 xi = p.x() + snxt*v.x() ;
1535 yi = p.y() + snxt*v.y() ;
1541 *validNorm = false ;
1552 *validNorm = false ;
1564 *validNorm = false ;
1581 std::ostringstream message;
1582 G4int oldprc = message.precision(16);
1583 message <<
"Undefined side for valid surface normal to solid."
1585 <<
"Position:" << G4endl << G4endl
1586 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1587 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1588 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1589 <<
"Direction:" << G4endl << G4endl
1590 <<
"v.x() = " << v.x() << G4endl
1591 <<
"v.y() = " << v.y() << G4endl
1592 <<
"v.z() = " << v.z() << G4endl << G4endl
1593 <<
"Proposed distance :" << G4endl << G4endl
1594 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
1595 message.precision(oldprc) ;
1596 G4Exception(
"G4Tubs::DistanceToOut(p,v,..)",
"GeomSolids1002",
1612 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1613 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1625 G4cout.precision(oldprc) ;
1626 G4Exception(
"G4Tubs::DistanceToOut(p)",
"GeomSolids1002",
1633 safeR1 = rho -
fRMin ;
1634 safeR2 =
fRMax - rho ;
1636 if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1637 else { safe = safeR2 ; }
1641 safe =
fRMax - rho ;
1643 safeZ =
fDz - std::fabs(p.z()) ;
1645 if ( safeZ < safe ) { safe = safeZ ; }
1659 if (safePhi < safe) { safe = safePhi ; }
1661 if ( safe < 0 ) { safe = 0 ; }
1682 G4double meshAngle, meshRMax, crossAngle,
1683 cosCrossAngle, sinCrossAngle, sAngle;
1684 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1685 G4int crossSection, noCrossSections;
1701 meshAngle =
fDPhi/(noCrossSections - 1) ;
1711 else { sAngle =
fSPhi ; }
1717 vertices->reserve(noCrossSections*4);
1718 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1722 crossAngle = sAngle + crossSection*meshAngle ;
1723 cosCrossAngle = std::cos(crossAngle) ;
1724 sinCrossAngle = std::sin(crossAngle) ;
1726 rMaxX = meshRMax*cosCrossAngle ;
1727 rMaxY = meshRMax*sinCrossAngle ;
1736 rMinX = meshRMin*cosCrossAngle ;
1737 rMinY = meshRMin*sinCrossAngle ;
1755 "Error in allocation of vertices. Out of memory !");
1784 G4int oldprc = os.precision(16);
1785 os <<
"-----------------------------------------------------------\n"
1786 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1787 <<
" ===================================================\n"
1788 <<
" Solid type: G4Tubs\n"
1789 <<
" Parameters: \n"
1790 <<
" inner radius : " <<
fRMin/
mm <<
" mm \n"
1791 <<
" outer radius : " <<
fRMax/
mm <<
" mm \n"
1792 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1793 <<
" starting phi : " <<
fSPhi/
degree <<
" degrees \n"
1794 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1795 <<
"-----------------------------------------------------------\n";
1796 os.precision(oldprc);
1807 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1808 aOne, aTwo, aThr, aFou;
1817 cosphi = std::cos(phi);
1818 sinphi = std::sin(phi);
1826 if( (chose >=0) && (chose < aOne) )
1828 xRand = fRMax*cosphi;
1829 yRand = fRMax*sinphi;
1833 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1835 xRand = fRMin*cosphi;
1836 yRand = fRMin*sinphi;
1840 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1842 xRand = rRand*cosphi;
1843 yRand = rRand*sinphi;
1847 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1849 xRand = rRand*cosphi;
1850 yRand = rRand*sinphi;
1854 else if( (chose >= aOne + aTwo + 2.*aThr)
1855 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1857 xRand = rRand*std::cos(
fSPhi);
1858 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)
static const double twopi
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