73 using namespace CLHEP;
120 std::ostringstream message;
121 message <<
"Invalid swept radius for Solid: " <<
GetName() <<
G4endl
122 <<
" pRtor = " << pRtor <<
", pRmax = " << pRmax;
129 if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
131 if (pRmin >= 1.
e2*kCarTolerance) {
fRmin = pRmin ; }
132 else {
fRmin = 0.0 ; }
137 std::ostringstream message;
138 message <<
"Invalid values of radii for Solid: " <<
GetName() <<
G4endl
139 <<
" pRmin = " << pRmin <<
", pRmax = " << pRmax;
152 if ( pDPhi >= twopi ) {
fDPhi = twopi ; }
155 if (pDPhi > 0) {
fDPhi = pDPhi ; }
158 std::ostringstream message;
159 message <<
"Invalid Z delta-Phi for Solid: " <<
GetName() <<
G4endl
160 <<
" pDPhi = " << pDPhi;
182 :
G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
183 fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
184 kRadTolerance(0.), kAngTolerance(0.),
185 halfCarTolerance(0.), halfAngTolerance(0.)
201 :
G4CSGSolid(rhs), fRmin(rhs.fRmin),fRmax(rhs.fRmax),
202 fRtor(rhs.fRtor),fSPhi(rhs.fSPhi),fDPhi(rhs.fDPhi),
203 fRminTolerance(rhs.fRminTolerance), fRmaxTolerance(rhs.fRmaxTolerance),
204 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
205 halfCarTolerance(rhs.halfCarTolerance),
206 halfAngTolerance(rhs.halfAngTolerance)
219 if (
this == &rhs) {
return *
this; }
260 std::vector<G4double>& roots )
const
268 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
269 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
273 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - r2 + 2*Rtor2*v.z()*v.z()) ;
274 c[3] = 4*(pDotV*(pRad2 - Rtor2 - r2) + 2*Rtor2*p.z()*v.z()) ;
275 c[4] = pRad2*pRad2 - 2*pRad2*(Rtor2+r2)
276 + 4*Rtor2*p.z()*p.z() + (Rtor2-r2)*(Rtor2-r2) ;
280 num = torusEq.
FindRoots( c, 4, srd, si );
282 for ( i = 0; i < num; i++ )
284 if( si[i] == 0. ) { roots.push_back(srd[i]) ; }
287 std::sort(roots.begin() , roots.end() ) ;
300 G4bool IsDistanceToIn )
const
309 std::vector<G4double> roots ;
310 std::vector<G4double> rootsrefined ;
317 for (
size_t k = 0 ; k<roots.size() ; k++ )
327 if ( rootsrefined.size()==roots.size() )
329 t = t + rootsrefined[k] ;
335 G4double theta = std::atan2(ptmp.y(),ptmp.x());
357 if ( IsDistanceToIn ==
true )
366 p.y()*(1-
fRtor/std::sqrt(p.x()*p.x()
372 if ( r ==
GetRmin() ) { scal = -scal ; }
373 if ( scal < 0 ) {
return 0.0 ; }
380 if ( IsDistanceToIn ==
false )
388 p.y()*(1-
fRtor/std::sqrt(p.x()*p.x()
394 if ( r ==
GetRmin() ) { scal = -scal ; }
395 if ( scal > 0 ) {
return 0.0 ; }
432 G4double RTorus,delta,diff1,diff2,maxDiff,newMin,newMax;
480 zMin = zoffset -
fRmax ;
481 zMax = zoffset +
fRmax ;
510 if ( yoff1 >= 0 && yoff2 >= 0 )
524 delta = RTorus*RTorus - yoff1*yoff1;
525 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
526 delta = RTorus*RTorus - yoff2*yoff2;
527 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
528 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
529 newMin = xoffset - maxDiff ;
530 newMax = xoffset + maxDiff ;
531 pMin = (newMin < xMin) ? xMin : newMin ;
532 pMax = (newMax > xMax) ? xMax : newMax ;
537 xoff1 = xoffset - xMin ;
538 xoff2 = xMax - xoffset ;
539 if (xoff1 >= 0 && xoff2 >= 0 )
552 delta = RTorus*RTorus - xoff1*xoff1;
553 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
554 delta = RTorus*RTorus - xoff2*xoff2;
555 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
556 maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
557 newMin = yoffset - maxDiff ;
558 newMax = yoffset + maxDiff ;
559 pMin = (newMin < yMin) ? yMin : newMin ;
560 pMax = (newMax > yMax) ? yMax : newMax ;
579 G4int i, noEntries, noBetweenSections4 ;
580 G4bool existsAfterClip = false ;
585 G4int noPolygonVertices ;
591 noEntries = vertices->size() ;
592 noBetweenSections4 = noEntries - noPolygonVertices ;
594 for (i=0;i<noEntries;i+=noPolygonVertices)
598 for (i=0;i<noBetweenSections4;i+=noPolygonVertices)
604 existsAfterClip = true ;
622 existsAfterClip = true ;
628 return existsAfterClip;
638 G4double r2, pt2, pPhi, tolRMin, tolRMax ;
644 r2 = p.x()*p.x() + p.y()*p.y() ;
652 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
654 if (
fDPhi == twopi || pt2 == 0 )
663 pPhi = std::atan2(p.y(),p.x()) ;
700 if (tolRMin < 0 ) { tolRMin = 0 ; }
702 if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
704 if ( (
fDPhi == twopi) || (pt2 == 0) )
710 pPhi = std::atan2(p.y(),p.x()) ;
749 G4int noSurfaces = 0;
763 rho2 = p.x()*p.x() + p.y()*p.y();
764 rho = std::sqrt(rho2);
765 pt2 = rho2+p.z()*p.z() +
fRtor * (
fRtor-2*rho);
767 pt = std::sqrt(pt2) ;
772 if( rho > delta && pt != 0.0 )
785 pPhi = std::atan2(p.y(),p.x());
787 if(pPhi <
fSPhi-delta) { pPhi += twopi; }
788 else if(pPhi >
fSPhi+
fDPhi+delta) { pPhi -= twopi; }
790 distSPhi = std::fabs( pPhi -
fSPhi );
796 if( distRMax <= delta )
801 else if(
fRmin && (distRMin <= delta) )
812 if (distSPhi <= dAngle)
817 if (distEPhi <= dAngle)
823 if ( noSurfaces == 0 )
833 ed <<
" ERROR> Surface Normal was called for Torus,"
834 <<
" with point not on surface." <<
G4endl;
838 ed <<
" ERROR> Surface Normal has not found a surface, "
839 <<
" despite the point being on the surface. " <<
G4endl;
850 ed <<
" Coordinates of point : " << p <<
G4endl;
851 ed <<
" Parameters of solid : " << G4endl << *
this <<
G4endl;
855 G4Exception(
"G4Torus::SurfaceNormal(p)",
"GeomSolids1002",
857 "Failing to find normal, even though point is on surface!" );
861 static const char* NameInside[3]= {
"Inside",
"Surface",
"Outside" };
862 ed <<
" The point is " << NameInside[inIt] <<
" the solid. "<<
G4endl;
863 G4Exception(
"G4Torus::SurfaceNormal(p)",
"GeomSolids1002",
869 else if ( noSurfaces == 1 ) { norm = sumnorm; }
870 else { norm = sumnorm.unit(); }
887 G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
889 rho2 = p.x()*p.x() + p.y()*p.y();
890 rho = std::sqrt(rho2) ;
892 pt = std::sqrt(pt2) ;
895 G4cout <<
" G4Torus::ApproximateSurfaceNormal called for point " << p
899 distRMax = std::fabs(pt -
fRmax) ;
903 distRMin = std::fabs(pt -
fRmin) ;
905 if (distRMin < distRMax)
921 if ( (
fDPhi < twopi) && rho )
923 phi = std::atan2(p.y(),p.x()) ;
925 if (phi < 0) { phi += twopi ; }
927 if (
fSPhi < 0 ) { distSPhi = std::fabs(phi-(
fSPhi+twopi))*rho ; }
928 else { distSPhi = std::fabs(phi-
fSPhi)*rho ; }
930 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
932 if (distSPhi < distEPhi)
934 if (distSPhi<distMin) side =
kNSPhi ;
938 if (distEPhi < distMin) { side =
kNEPhi ; }
945 -p.y()*(1-
fRtor/rho)/pt,
950 p.y()*(1-
fRtor/rho)/pt,
963 "Undefined side for valid surface normal to solid.");
1004 G4double cPhi,sinCPhi=0.,cosCPhi=0.;
1017 if (
fDPhi < twopi )
1021 cPhi =
fSPhi + hDPhi ;
1022 sinCPhi = std::sin(cPhi) ;
1023 cosCPhi = std::cos(cPhi) ;
1049 if ( sd[0] < snxt ) { snxt = sd[0] ; }
1064 sinSPhi = std::sin(
fSPhi) ;
1065 cosSPhi = std::cos(
fSPhi) ;
1066 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
1070 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1077 if ( sphi < 0 ) { sphi = 0 ; }
1079 xi = p.x() + sphi*v.x() ;
1080 yi = p.y() + sphi*v.y() ;
1081 zi = p.z() + sphi*v.z() ;
1082 rhoi2 = xi*xi + yi*yi ;
1083 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
1085 if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
1090 if ((yi*cosCPhi-xi*sinCPhi)<=0) { snxt=sphi; }
1096 sinEPhi=std::sin(ePhi);
1097 cosEPhi=std::cos(ePhi);
1098 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
1102 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1110 if (sphi < 0 ) { sphi = 0 ; }
1112 xi = p.x() + sphi*v.x() ;
1113 yi = p.y() + sphi*v.y() ;
1114 zi = p.z() + sphi*v.z() ;
1115 rhoi2 = xi*xi + yi*yi ;
1116 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
1118 if (it2 >= tolORMin2 && it2 <= tolORMax2)
1123 if ((yi*cosCPhi-xi*sinCPhi)>=0) { snxt=sphi; }
1144 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1147 rho2 = p.x()*p.x() + p.y()*p.y() ;
1148 rho = std::sqrt(rho2) ;
1150 pt = std::sqrt(pt2) ;
1152 safe1 =
fRmin - pt ;
1153 safe2 = pt -
fRmax ;
1155 if (safe1 > safe2) { safe = safe1; }
1156 else { safe = safe2; }
1158 if (
fDPhi < twopi && rho )
1161 cosPhiC = std::cos(phiC) ;
1162 sinPhiC = std::sin(phiC) ;
1163 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1165 if (cosPsi < std::cos(
fDPhi*0.5) )
1167 if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1169 safePhi = std::fabs(p.x()*std::sin(
fSPhi) - p.y()*std::cos(
fSPhi)) ;
1174 safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1176 if (safePhi > safe) { safe = safePhi ; }
1179 if (safe < 0 ) { safe = 0 ; }
1200 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
1202 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
1214 G4double rho2 = p.x()*p.x()+p.y()*p.y();
1222 pt2= std::fabs( pt2 );
1227 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
1231 G4double vDotNmax = pDotV -
fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
1234 if( (pt2 > tolRMax*tolRMax) && (vDotNmax >= 0) )
1240 if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) )
1243 p.y()*(1 -
fRtor/rho)/pt,
1260 if ( (pt2 < tolRMin*tolRMin) && (vDotNmax < 0) )
1262 if (calcNorm) { *validNorm = false ; }
1292 if ( calcNorm && (snxt == 0.0) )
1294 *validNorm = false ;
1302 sinSPhi = std::sin(
fSPhi) ;
1303 cosSPhi = std::cos(
fSPhi) ;
1305 sinEPhi = std::sin(ePhi) ;
1306 cosEPhi = std::cos(ePhi) ;
1307 cPhi =
fSPhi + fDPhi*0.5 ;
1308 sinCPhi = std::sin(cPhi) ;
1309 cosCPhi = std::cos(cPhi) ;
1314 vphi = std::atan2(v.y(),v.x()) ;
1319 if ( p.x() || p.y() )
1321 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1322 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1326 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
1327 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
1339 sphi = pDistS/compS ;
1343 xi = p.x() + sphi*v.x() ;
1344 yi = p.y() + sphi*v.y() ;
1359 else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1380 sphi2 = pDistE/compE ;
1386 xi = p.x() + sphi2*v.x() ;
1387 yi = p.y() + sphi2*v.y() ;
1403 if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1425 vphi = std::atan2(v.y(),v.x());
1448 G4double rhoi2,rhoi,it2,it,iDotxyNmax ;
1457 xi = p.x() + snxt*v.x() ;
1458 yi =p.y() + snxt*v.y() ;
1459 zi = p.z() + snxt*v.z() ;
1460 rhoi2 = xi*xi + yi*yi ;
1461 rhoi = std::sqrt(rhoi2) ;
1463 it = std::sqrt(it2) ;
1464 iDotxyNmax = (1-
fRtor/rhoi) ;
1465 if(iDotxyNmax >= -2.*fRmaxTolerance)
1468 yi*(1-
fRtor/rhoi)/it,
1474 *validNorm = false ;
1479 *validNorm = false ;
1490 *validNorm = false ;
1502 *validNorm = false ;
1512 std::ostringstream message;
1513 G4int oldprc = message.precision(16);
1514 message <<
"Undefined side for valid surface normal to solid."
1516 <<
"Position:" << G4endl << G4endl
1517 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1518 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1519 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1520 <<
"Direction:" << G4endl << G4endl
1521 <<
"v.x() = " << v.x() << G4endl
1522 <<
"v.y() = " << v.y() << G4endl
1523 <<
"v.z() = " << v.z() << G4endl << G4endl
1524 <<
"Proposed distance :" << G4endl << G4endl
1525 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl;
1526 message.precision(oldprc);
1545 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1546 rho2 = p.x()*p.x() + p.y()*p.y() ;
1547 rho = std::sqrt(rho2) ;
1549 pt = std::sqrt(pt2) ;
1561 G4cout.precision(oldprc);
1562 G4Exception(
"G4Torus::DistanceToOut(p)",
"GeomSolids1002",
1569 safeR1 = pt -
fRmin ;
1570 safeR2 =
fRmax - pt ;
1572 if (safeR1 < safeR2) { safe = safeR1 ; }
1573 else { safe = safeR2 ; }
1585 cosPhiC = std::cos(phiC) ;
1586 sinPhiC = std::sin(phiC) ;
1588 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
1590 safePhi = -(p.x()*std::sin(
fSPhi) - p.y()*std::cos(
fSPhi)) ;
1595 safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1597 if (safePhi < safe) { safe = safePhi ; }
1599 if (safe < 0) { safe = 0 ; }
1616 G4int& noPolygonVertices )
const
1620 G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle;
1622 G4int crossSection,noCrossSections;
1636 meshAngle =
fDPhi/(noCrossSections - 1) ;
1637 meshRMax = (
fRtor +
fRmax)/std::cos(meshAngle*0.5) ;
1644 sAngle = -meshAngle*0.5 ;
1654 vertices->reserve(noCrossSections*4) ;
1655 for (crossSection=0;crossSection<noCrossSections;crossSection++)
1659 crossAngle=sAngle+crossSection*meshAngle;
1660 cosCrossAngle=std::cos(crossAngle);
1661 sinCrossAngle=std::sin(crossAngle);
1663 rMaxX=meshRMax*cosCrossAngle;
1664 rMaxY=meshRMax*sinCrossAngle;
1677 noPolygonVertices = 4 ;
1684 "Error in allocation of vertices. Out of memory !");
1713 G4int oldprc = os.precision(16);
1714 os <<
"-----------------------------------------------------------\n"
1715 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1716 <<
" ===================================================\n"
1717 <<
" Solid type: G4Torus\n"
1718 <<
" Parameters: \n"
1719 <<
" inner radius: " <<
fRmin/
mm <<
" mm \n"
1720 <<
" outer radius: " <<
fRmax/
mm <<
" mm \n"
1721 <<
" swept radius: " <<
fRtor/
mm <<
" mm \n"
1722 <<
" starting phi: " <<
fSPhi/
degree <<
" degrees \n"
1723 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1724 <<
"-----------------------------------------------------------\n";
1725 os.precision(oldprc);
1736 G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
1741 cosu = std::cos(phi); sinu = std::sin(phi);
1742 cosv = std::cos(theta); sinv = std::sin(theta);
1750 if ((
fSPhi == 0) && (
fDPhi == twopi)){ aSide = 0; }
1758 else if( (chose >= aOut) && (chose < aOut + aIn) )
1763 else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
1767 (
fRtor+rRand*cosv)*std::sin(
fSPhi), rRand*sinv);
EInside Inside(const G4ThreeVector &p) const
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)
G4Polyhedron * CreatePolyhedron() const
static const G4double kInfinity
G4double GetMinYExtent() const
std::ostringstream G4ExceptionDescription
CLHEP::Hep3Vector G4ThreeVector
G4Torus & operator=(const G4Torus &rhs)
std::ostream & StreamInfo(std::ostream &os) const
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double SolveNumericJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, G4bool IsDistanceToIn) const
G4bool IsXLimited() const
const G4int kMinMeshSections
virtual void AddSolid(const G4Box &)=0
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4double GetMaxXExtent() const
G4double GetMinZExtent() const
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform, G4int &noPolygonVertices) const
G4GLOB_DLL std::ostream G4cout
G4Polyhedron * fpPolyhedron
G4double GetRadialTolerance() const
virtual G4Polyhedron * GetPolyhedron() const
std::vector< G4ThreeVector > G4ThreeVectorList
G4Torus(const G4String &pName, G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4ThreeVector GetPointOnSurface() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4double GetMinXExtent() const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double halfCarTolerance
G4double GetMaxZExtent() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double halfAngTolerance
static const double degree
G4GeometryType GetEntityType() const
void TorusRootsJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, std::vector< G4double > &roots) const
G4double GetMaxYExtent() const
void SetAllParameters(G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
G4double GetMaxExtent(const EAxis pAxis) const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
G4bool IsZLimited() const
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
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()
const G4int kMaxMeshSections
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const