66 #if !defined(G4GEOM_USE_UTUBS)
82 using namespace CLHEP;
93 :
G4CSGSolid(pName), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
106 message <<
"Negative Z half-length (" << pDz <<
") in solid: " <<
GetName();
109 if ( (pRMin >= pRMax) || (pRMin < 0) )
112 message <<
"Invalid values for radii in solid: " <<
GetName()
114 <<
" pRMin = " << pRMin <<
", pRMax = " << pRMax;
129 :
G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
130 fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
131 sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
132 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
133 fPhiFullTube(false), halfCarTolerance(0.), halfRadTolerance(0.),
152 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
153 fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
154 fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
155 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
156 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
157 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
158 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi), fPhiFullTube(rhs.fPhiFullTube),
159 halfCarTolerance(rhs.halfCarTolerance),
160 halfRadTolerance(rhs.halfRadTolerance),
161 halfAngTolerance(rhs.halfAngTolerance)
173 if (
this == &rhs) {
return *
this; }
227 pMin.
set(vmin.
x(),vmin.
y(),-dz);
228 pMax.
set(vmax.
x(),vmax.
y(), dz);
232 pMin.
set(-rmax,-rmax,-dz);
233 pMax.
set( rmax, rmax, dz);
238 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
241 message <<
"Bad bounding box (min >= max) for solid: "
243 <<
"\npMin = " << pMin
244 <<
"\npMax = " << pMax;
245 G4Exception(
"G4Tubs::BoundingLimits()",
"GeomMgt0001",
270 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
274 return exist = (pMin < pMax) ?
true :
false;
285 const G4int NSTEPS = 24;
287 G4int ksteps = (dphi <= astep) ? 1 : (
G4int)((dphi-
deg)/astep) + 1;
290 G4double sinHalf = std::sin(0.5*ang);
291 G4double cosHalf = std::cos(0.5*ang);
292 G4double sinStep = 2.*sinHalf*cosHalf;
293 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
298 if (rmin == 0 && dphi ==
twopi)
304 for (
G4int k=0; k<NSTEPS; ++k)
306 baseA[k].set(rext*cosCur,rext*sinCur,-dz);
307 baseB[k].set(rext*cosCur,rext*sinCur, dz);
310 sinCur = sinCur*cosStep + cosCur*sinStep;
311 cosCur = cosCur*cosStep - sinTmp*sinStep;
313 std::vector<const G4ThreeVectorList *> polygons(2);
314 polygons[0] = &baseA;
315 polygons[1] = &baseB;
325 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
326 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
330 for (
G4int k=0; k<ksteps+2; ++k) pols[k].resize(4);
331 pols[0][0].set(rmin*cosStart,rmin*sinStart, dz);
332 pols[0][1].set(rmin*cosStart,rmin*sinStart,-dz);
333 pols[0][2].set(rmax*cosStart,rmax*sinStart,-dz);
334 pols[0][3].set(rmax*cosStart,rmax*sinStart, dz);
335 for (
G4int k=1; k<ksteps+1; ++k)
337 pols[k][0].set(rmin*cosCur,rmin*sinCur, dz);
338 pols[k][1].set(rmin*cosCur,rmin*sinCur,-dz);
339 pols[k][2].set(rext*cosCur,rext*sinCur,-dz);
340 pols[k][3].set(rext*cosCur,rext*sinCur, dz);
343 sinCur = sinCur*cosStep + cosCur*sinStep;
344 cosCur = cosCur*cosStep - sinTmp*sinStep;
346 pols[ksteps+1][0].set(rmin*cosEnd,rmin*sinEnd, dz);
347 pols[ksteps+1][1].set(rmin*cosEnd,rmin*sinEnd,-dz);
348 pols[ksteps+1][2].set(rmax*cosEnd,rmax*sinEnd,-dz);
349 pols[ksteps+1][3].set(rmax*cosEnd,rmax*sinEnd, dz);
352 std::vector<const G4ThreeVectorList *> polygons;
353 polygons.resize(ksteps+2);
354 for (
G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
372 r2 = p.
x()*p.
x() + p.
y()*p.
y() ;
375 else { tolRMin = 0 ; }
379 if ((r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax))
397 pPhi = std::atan2(p.
y(),p.
x()) ;
440 if ( tolRMin < 0 ) { tolRMin = 0; }
442 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
444 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance) )
450 pPhi = std::atan2(p.
y(),p.
x()) ;
481 r2 = p.
x()*p.
x() + p.
y()*p.
y() ;
485 if ( tolRMin < 0 ) { tolRMin = 0; }
487 if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
489 if (
fPhiFullTube || (r2 <=halfRadTolerance*halfRadTolerance))
495 pPhi = std::atan2(p.
y(),p.
x()) ;
534 G4int noSurfaces = 0;
543 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y());
545 distRMin = std::fabs(rho -
fRMin);
546 distRMax = std::fabs(rho -
fRMax);
547 distZ = std::fabs(std::fabs(p.
z()) -
fDz);
553 pPhi = std::atan2(p.
y(),p.
x());
558 distSPhi = std::fabs(pPhi -
fSPhi);
597 if ( p.
z() >= 0.) { sumnorm += nZ; }
598 else { sumnorm -= nZ; }
600 if ( noSurfaces == 0 )
603 G4Exception(
"G4Tubs::SurfaceNormal(p)",
"GeomSolids1002",
606 G4cout<<
"G4Tubs::SN ( "<<p.
x()<<
", "<<p.
y()<<
", "<<p.
z()<<
" ); "
608 G4cout.precision(oldprc) ;
612 else if ( noSurfaces == 1 ) { norm = sumnorm; }
613 else { norm = sumnorm.
unit(); }
628 G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
630 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
632 distRMin = std::fabs(rho -
fRMin) ;
633 distRMax = std::fabs(rho -
fRMax) ;
634 distZ = std::fabs(std::fabs(p.
z()) -
fDz) ;
636 if (distRMin < distRMax)
638 if ( distZ < distRMin )
651 if ( distZ < distRMax )
664 phi = std::atan2(p.
y(),p.
x()) ;
666 if ( phi < 0 ) { phi +=
twopi; }
670 distSPhi = std::fabs(phi - (
fSPhi +
twopi))*rho ;
674 distSPhi = std::fabs(phi -
fSPhi)*rho ;
676 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
678 if (distSPhi < distEPhi)
680 if ( distSPhi < distMin )
687 if ( distEPhi < distMin )
726 "Undefined side for valid surface normal to solid.");
760 G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
765 G4double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp ;
788 if (std::fabs(p.
z()) >= tolIDz)
790 if ( p.
z()*v.
z() < 0 )
792 sd = (std::fabs(p.
z()) -
fDz)/std::fabs(v.
z()) ;
794 if(sd < 0.0) { sd = 0.0; }
796 xi = p.
x() + sd*v.
x() ;
797 yi = p.
y() + sd*v.
y() ;
798 rho2 = xi*xi + yi*yi ;
802 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
809 iden = std::sqrt(rho2) ;
821 if ( snxt<halfCarTolerance ) { snxt=0; }
838 t1 = 1.0 - v.
z()*v.
z() ;
839 t2 = p.
x()*v.
x() + p.
y()*v.
y() ;
840 t3 = p.
x()*p.
x() + p.
y()*p.
y() ;
846 if ((t3 >= tolORMax2) && (t2<0))
856 sd = c/(-b+std::sqrt(d));
861 G4double fTerm = sd-std::fmod(sd,dRmax);
866 zi = p.
z() + sd*v.
z() ;
867 if (std::fabs(zi)<=tolODz)
877 xi = p.
x() + sd*v.
x() ;
878 yi = p.
y() + sd*v.
y() ;
891 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.
z()) <= tolIDz))
898 iden = std::sqrt(t3) ;
919 snxt = c/(-b+std::sqrt(d));
921 if ( snxt < halfCarTolerance ) { snxt=0; }
939 c = t3 - fRMax*
fRMax;
950 snxt= c/(-b+std::sqrt(d));
952 if ( snxt < halfCarTolerance ) { snxt=0; }
972 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
973 if (sd >= -halfCarTolerance)
977 if(sd < 0.0) { sd = 0.0; }
980 G4double fTerm = sd-std::fmod(sd,dRmax);
983 zi = p.
z() + sd*v.
z() ;
984 if (std::fabs(zi) <= tolODz)
994 xi = p.
x() + sd*v.
x() ;
995 yi = p.
y() + sd*v.
y() ;
1030 if ( Dist < halfCarTolerance )
1036 if ( sd < 0 ) { sd = 0.0; }
1037 zi = p.
z() + sd*v.
z() ;
1038 if ( std::fabs(zi) <= tolODz )
1040 xi = p.
x() + sd*v.
x() ;
1041 yi = p.
y() + sd*v.
y() ;
1042 rho2 = xi*xi + yi*yi ;
1044 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1045 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1048 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1070 if ( Dist < halfCarTolerance )
1076 if ( sd < 0 ) { sd = 0; }
1077 zi = p.
z() + sd*v.
z() ;
1078 if ( std::fabs(zi) <= tolODz )
1080 xi = p.
x() + sd*v.
x() ;
1081 yi = p.
y() + sd*v.
y() ;
1082 rho2 = xi*xi + yi*yi ;
1083 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1084 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1087 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1101 if ( snxt<halfCarTolerance ) { snxt=0; }
1133 G4double safe=0.0, rho, safe1, safe2, safe3 ;
1136 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
1137 safe1 =
fRMin - rho ;
1138 safe2 = rho -
fRMax ;
1139 safe3 = std::fabs(p.
z()) -
fDz ;
1141 if ( safe1 > safe2 ) { safe = safe1; }
1142 else { safe = safe2; }
1143 if ( safe3 > safe ) { safe = safe3; }
1151 if ( cosPsi < std::cos(
fDPhi*0.5) )
1163 if ( safePhi > safe ) { safe = safePhi; }
1166 if ( safe < 0 ) { safe = 0; }
1187 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1193 pdist =
fDz - p.
z() ;
1196 snxt = pdist/v.
z() ;
1209 else if ( v.
z() < 0 )
1211 pdist =
fDz + p.
z() ;
1215 snxt = -pdist/v.
z() ;
1245 t1 = 1.0 - v.
z()*v.
z() ;
1246 t2 = p.
x()*v.
x() + p.
y()*v.
y() ;
1247 t3 = p.
x()*p.
x() + p.
y()*p.
y() ;
1250 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }
1270 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1289 roMin2 = t3 - t2*t2/
t1 ;
1305 srd = c/(-b+std::sqrt(d2));
1310 if ( calcNorm ) { *validNorm =
false; }
1321 srd = -b + std::sqrt(d2) ;
1345 srd = -b + std::sqrt(d2) ;
1368 vphi = std::atan2(v.
y(),v.
x()) ;
1374 if ( p.
x() || p.
y() )
1397 sphi = pDistS/compS ;
1401 xi = p.
x() + sphi*v.
x() ;
1402 yi = p.
y() + sphi*v.
y() ;
1441 sphi2 = pDistE/compE ;
1447 xi = p.
x() + sphi2*v.
x() ;
1448 yi = p.
y() + sphi2*v.
y() ;
1459 else { sphi = 0.0 ; }
1470 else { sphi = 0.0 ; }
1516 xi = p.
x() + snxt*v.
x() ;
1517 yi = p.
y() + snxt*v.
y() ;
1523 *validNorm = false ;
1534 *validNorm = false ;
1546 *validNorm = false ;
1564 G4int oldprc = message.precision(16);
1565 message <<
"Undefined side for valid surface normal to solid."
1567 <<
"Position:" << G4endl << G4endl
1568 <<
"p.x() = " << p.
x()/
mm <<
" mm" << G4endl
1569 <<
"p.y() = " << p.
y()/
mm <<
" mm" << G4endl
1570 <<
"p.z() = " << p.
z()/
mm <<
" mm" << G4endl << G4endl
1571 <<
"Direction:" << G4endl << G4endl
1572 <<
"v.x() = " << v.
x() << G4endl
1573 <<
"v.y() = " << v.
y() << G4endl
1574 <<
"v.z() = " << v.
z() << G4endl << G4endl
1575 <<
"Proposed distance :" << G4endl << G4endl
1576 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
1577 message.precision(oldprc) ;
1578 G4Exception(
"G4Tubs::DistanceToOut(p,v,..)",
"GeomSolids1002",
1594 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi ;
1595 rho = std::sqrt(p.
x()*p.
x() + p.
y()*p.
y()) ;
1607 G4cout.precision(oldprc) ;
1608 G4Exception(
"G4Tubs::DistanceToOut(p)",
"GeomSolids1002",
1615 safeR1 = rho -
fRMin ;
1616 safeR2 =
fRMax - rho ;
1618 if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1619 else { safe = safeR2 ; }
1623 safe =
fRMax - rho ;
1625 safeZ =
fDz - std::fabs(p.
z()) ;
1627 if ( safeZ < safe ) { safe = safeZ ; }
1641 if (safePhi < safe) { safe = safePhi ; }
1643 if ( safe < 0 ) { safe = 0 ; }
1663 return new G4Tubs(*
this);
1672 G4int oldprc = os.precision(16);
1673 os <<
"-----------------------------------------------------------\n"
1674 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1675 <<
" ===================================================\n"
1676 <<
" Solid type: G4Tubs\n"
1677 <<
" Parameters: \n"
1678 <<
" inner radius : " <<
fRMin/
mm <<
" mm \n"
1679 <<
" outer radius : " <<
fRMax/
mm <<
" mm \n"
1680 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1681 <<
" starting phi : " <<
fSPhi/
degree <<
" degrees \n"
1682 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1683 <<
"-----------------------------------------------------------\n";
1684 os.precision(oldprc);
1695 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1696 aOne, aTwo, aThr, aFou;
1705 cosphi = std::cos(phi);
1706 sinphi = std::sin(phi);
1714 if( (chose >=0) && (chose < aOne) )
1716 xRand = fRMax*cosphi;
1717 yRand = fRMax*sinphi;
1721 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1723 xRand = fRMin*cosphi;
1724 yRand = fRMin*sinphi;
1728 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1730 xRand = rRand*cosphi;
1731 yRand = rRand*sinphi;
1735 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1737 xRand = rRand*cosphi;
1738 yRand = rRand*sinphi;
1742 else if( (chose >= aOne + aTwo + 2.*aThr)
1743 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1745 xRand = rRand*std::cos(
fSPhi);
1746 yRand = rRand*std::sin(
fSPhi);
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void set(double x, double y, double z)
EInside Inside(const G4ThreeVector &p) const
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
virtual G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
G4ThreeVector GetPointOnSurface() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
static constexpr double mm
G4double GetSinStartPhi() const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4Tubs & operator=(const G4Tubs &rhs)
void message(RunManager *runmanager)
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
G4double GetAngularTolerance() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double GetRadialTolerance() const
G4double GetSinEndPhi() const
G4double GetDeltaPhiAngle() const
void CheckPhiAngles(G4double sPhi, G4double dPhi)
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
std::vector< G4ThreeVector > G4ThreeVectorList
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4double halfCarTolerance
static constexpr double deg
G4GeometryType GetEntityType() const
G4double GetInnerRadius() const
std::ostream & StreamInfo(std::ostream &os) const
G4double halfRadTolerance
static constexpr double twopi
G4double GetZHalfLength() const
ThreeVector shoot(const G4int Ap, const G4int Af)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double halfAngTolerance
G4Polyhedron * CreatePolyhedron() const
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
static constexpr double degree
G4double GetOuterRadius() const
G4GLOB_DLL std::ostream G4cout
static G4GeometryTolerance * GetInstance()
G4double GetCosEndPhi() const
static constexpr double pi
G4double GetCosStartPhi() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4Tubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi)