54 using namespace CLHEP;
66 :
G4OTubs(pName, pRMin, pRMax, pDz, pSPhi, pDPhi),
81 if ( ( !pLowNorm.x()) && ( !pLowNorm.y())
82 && ( !pHighNorm.x()) && (!pHighNorm.y()) )
84 std::ostringstream message;
85 message <<
"Inexisting Low/High Normal to Z plane or Parallel to Z."
87 <<
"Normals to Z plane are (" << pLowNorm <<
" and "
88 << pHighNorm <<
") in solid: " <<
GetName();
89 G4Exception(
"G4CutTubs::G4CutTubs()",
"GeomSolids1001",
95 if (pLowNorm.mag2() == 0.) { pLowNorm.setZ(-1.); }
96 if (pHighNorm.mag2()== 0.) { pHighNorm.setZ(1.); }
101 if (pLowNorm.mag2() != 1.) { pLowNorm = pLowNorm.unit(); }
102 if (pHighNorm.mag2()!= 1.) { pHighNorm = pHighNorm.unit(); }
106 if( (pLowNorm.mag2() != 0.) && (pHighNorm.mag2()!= 0. ) )
108 if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z() <= 0.))
110 std::ostringstream message;
111 message <<
"Invalid Low or High Normal to Z plane; "
112 "has to point outside Solid." <<
G4endl
113 <<
"Invalid Norm to Z plane (" << pLowNorm <<
" or "
114 << pHighNorm <<
") in solid: " <<
GetName();
115 G4Exception(
"G4CutTubs::G4CutTubs()",
"GeomSolids0002",
126 std::ostringstream message;
127 message <<
"Invalid Low or High Normal to Z plane; "
128 <<
"Crossing Cutted Planes." <<
G4endl
129 <<
"Invalid Norm to Z plane (" << pLowNorm <<
" and "
130 << pHighNorm <<
") in solid: " <<
GetName();
131 G4Exception(
"G4CutTubs::G4CutTubs()",
"GeomSolids0002",
144 halfCarTolerance(0.), halfRadTolerance(0.), halfAngTolerance(0.)
161 :
G4OTubs(rhs), fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm),
162 fPhiFullCutTube(rhs.fPhiFullCutTube),
163 halfCarTolerance(rhs.halfCarTolerance),
164 halfRadTolerance(rhs.halfRadTolerance),
165 halfAngTolerance(rhs.halfAngTolerance)
178 if (
this == &rhs) {
return *
this; }
218 G4double diff1, diff2, maxDiff, newMin, newMax;
219 G4double xoff1, xoff2, yoff1, yoff2, delta;
222 xMin = xoffset -
fRMax;
223 xMax = xoffset +
fRMax;
245 yMin = yoffset -
fRMax;
246 yMax = yoffset +
fRMax;
295 yoff1 = yoffset - yMin;
296 yoff2 = yMax - yoffset;
298 if ( (yoff1 >= 0) && (yoff2 >= 0) )
308 delta = fRMax*fRMax - yoff1*yoff1;
309 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
310 delta = fRMax*fRMax - yoff2*yoff2;
311 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
312 maxDiff = (diff1 > diff2) ? diff1:diff2;
313 newMin = xoffset - maxDiff;
314 newMax = xoffset + maxDiff;
315 pMin = (newMin < xMin) ? xMin : newMin;
316 pMax = (newMax > xMax) ? xMax : newMax;
322 xoff1 = xoffset - xMin;
323 xoff2 = xMax - xoffset;
325 if ( (xoff1 >= 0) && (xoff2 >= 0) )
335 delta = fRMax*fRMax - xoff1*xoff1;
336 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
337 delta = fRMax*fRMax - xoff2*xoff2;
338 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
339 maxDiff = (diff1 > diff2) ? diff1 : diff2;
340 newMin = yoffset - maxDiff;
341 newMax = yoffset + maxDiff;
342 pMin = (newMin < yMin) ? yMin : newMin;
343 pMax = (newMax > yMax) ? yMax : newMax;
362 G4int i, noEntries, noBetweenSections4;
363 G4bool existsAfterClip =
false;
369 noEntries = vertices->size();
370 noBetweenSections4 = noEntries - 4;
372 for ( i = 0 ; i < noEntries ; i += 4 )
376 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
382 existsAfterClip =
true;
400 existsAfterClip =
true;
406 return existsAfterClip;
435 r2 = p.x()*p.x() + p.y()*p.y() ;
441 if ( tolRMin < 0 ) { tolRMin = 0; }
443 if ( ((r2 < tolRMin*tolRMin) || (r2 > tolRMax*tolRMax))
444 && (r2 >=halfRadTolerance*halfRadTolerance) ) {
return kOutside; }
458 pPhi = std::atan2(p.y(),p.x()) ;
504 else { tolRMin = 0 ; }
506 if ( ((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax))&&
523 G4int noSurfaces = 0;
525 G4double distZLow,distZHigh, distRMin, distRMax;
533 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
535 distRMin = std::fabs(rho -
fRMin);
536 distRMax = std::fabs(rho -
fRMax);
540 distZLow =std::fabs((p+vZ).dot(
fLowNorm));
544 distZHigh = std::fabs((p-vZ).dot(
fHighNorm));
550 pPhi = std::atan2(p.y(),p.x());
555 distSPhi = std::fabs(pPhi -
fSPhi);
601 if ( noSurfaces == 0 )
604 G4Exception(
"G4CutTubs::SurfaceNormal(p)",
"GeomSolids1002",
607 G4cout<<
"G4CutTubs::SN ( "<<p.x()<<
", "<<p.y()<<
", "<<p.z()<<
" ); "
609 G4cout.precision(oldprc) ;
613 else if ( noSurfaces == 1 ) { norm = sumnorm; }
614 else { norm = sumnorm.unit(); }
630 G4double distRMin, distRMax, distSPhi, distEPhi, distMin ;
633 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
635 distRMin = std::fabs(rho -
fRMin) ;
636 distRMax = std::fabs(rho -
fRMax) ;
640 distZLow =std::fabs((p+vZ).dot(
fLowNorm));
644 distZHigh = std::fabs((p-vZ).dot(
fHighNorm));
647 if (distRMin < distRMax)
649 if ( distZ < distRMin )
662 if ( distZ < distRMax )
675 phi = std::atan2(p.y(),p.x()) ;
677 if ( phi < 0 ) { phi += twopi; }
681 distSPhi = std::fabs(phi - (
fSPhi + twopi))*rho ;
685 distSPhi = std::fabs(phi -
fSPhi)*rho ;
687 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
689 if (distSPhi < distEPhi)
691 if ( distSPhi < distMin )
698 if ( distEPhi < distMin )
718 if ( distZHigh > distZLow ) { norm =
fHighNorm ; }
737 "Undefined side for valid surface normal to solid.");
777 G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ;
811 if(sd < 0.0) { sd = 0.0; }
813 xi = p.x() + sd*v.x() ;
814 yi = p.y() + sd*v.y() ;
815 rho2 = xi*xi + yi*yi ;
819 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
826 iden = std::sqrt(rho2) ;
851 sd = -distZHigh/calf;
853 if(sd < 0.0) { sd = 0.0; }
855 xi = p.x() + sd*v.x() ;
856 yi = p.y() + sd*v.y() ;
857 rho2 = xi*xi + yi*yi ;
861 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
868 iden = std::sqrt(rho2) ;
899 t1 = 1.0 - v.z()*v.z() ;
900 t2 = p.x()*v.x() + p.y()*v.y() ;
901 t3 = p.x()*p.x() + p.y()*p.y() ;
907 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 xi = p.x() + sd*v.x() ;
928 yi = p.y() + sd*v.y() ;
943 xi = p.x() + sd*v.x() ;
944 yi = p.y() + sd*v.y() ;
957 if ((t3 > tolIRMin2) && (t2 < 0)
965 iden = std::sqrt(t3) ;
986 snxt = c/(-b+std::sqrt(d));
1006 c = t3 - fRMax*
fRMax;
1017 snxt= c/(-b+std::sqrt(d));
1040 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1045 if (sd < 0.0) { sd = 0.0; }
1048 G4double fTerm = sd-std::fmod(sd,dRmax);
1051 zi = p.z() + sd*v.z() ;
1052 xi = p.x() + sd*v.x() ;
1053 yi = p.y() + sd*v.y() ;
1109 if ( sd < 0 ) { sd = 0.0; }
1110 zi = p.z() + sd*v.z() ;
1111 xi = p.x() + sd*v.x() ;
1112 yi = p.y() + sd*v.y() ;
1119 rho2 = xi*xi + yi*yi ;
1120 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1121 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1124 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1153 if ( sd < 0 ) { sd = 0; }
1154 zi = p.z() + sd*v.z() ;
1155 xi = p.x() + sd*v.x() ;
1156 yi = p.y() + sd*v.y() ;
1163 xi = p.x() + sd*v.x() ;
1164 yi = p.y() + sd*v.y() ;
1165 rho2 = xi*xi + yi*yi ;
1166 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1167 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1170 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1221 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1226 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1228 safRMin =
fRMin- rho ;
1229 safRMax = rho -
fRMax ;
1243 if ( safRMin > safe ) { safe = safRMin; }
1244 if ( safRMax> safe ) { safe = safRMax; }
1254 if ( cosPsi < std::cos(
fDPhi*0.5) )
1266 if ( safePhi > safe ) { safe = safePhi; }
1269 if ( safe < 0 ) { safe = 0; }
1287 G4double deltaR, t1, t2, t3, b, c,
d2, roMin2 ;
1288 G4double distZLow,distZHigh,calfH,calfL;
1293 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1313 snxt = -distZHigh/calfH ;
1331 sz = -distZLow/calfL ;
1348 if((calfH<=0)&&(calfL<=0))
1364 t1 = 1.0 - v.z()*v.z() ;
1365 t2 = p.x()*v.x() + p.y()*v.y() ;
1366 t3 = p.x()*p.x() + p.y()*p.y() ;
1369 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }
1389 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1408 roMin2 = t3 - t2*t2/t1 ;
1424 srd = c/(-b+std::sqrt(d2));
1429 if ( calcNorm ) { *validNorm =
false; }
1440 srd = -b + std::sqrt(d2) ;
1464 srd = -b + std::sqrt(d2) ;
1486 vphi = std::atan2(v.y(),v.x()) ;
1492 if ( p.x() || p.y() )
1515 sphi = pDistS/compS ;
1519 xi = p.x() + sphi*v.x() ;
1520 yi = p.y() + sphi*v.y() ;
1560 sphi2 = pDistE/compE ;
1566 xi = p.x() + sphi2*v.x() ;
1567 yi = p.y() + sphi2*v.y() ;
1578 else { sphi = 0.0 ; }
1589 else { sphi = 0.0 ; }
1635 xi = p.x() + snxt*v.x() ;
1636 yi = p.y() + snxt*v.y() ;
1642 *validNorm = false ;
1653 *validNorm = false ;
1665 *validNorm = false ;
1682 std::ostringstream message;
1683 G4int oldprc = message.precision(16);
1684 message <<
"Undefined side for valid surface normal to solid."
1686 <<
"Position:" << G4endl << G4endl
1687 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1688 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1689 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1690 <<
"Direction:" << G4endl << G4endl
1691 <<
"v.x() = " << v.x() << G4endl
1692 <<
"v.y() = " << v.y() << G4endl
1693 <<
"v.z() = " << v.z() << G4endl << G4endl
1694 <<
"Proposed distance :" << G4endl << G4endl
1695 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
1696 message.precision(oldprc) ;
1697 G4Exception(
"G4CutTubs::DistanceToOut(p,v,..)",
"GeomSolids1002",
1712 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1715 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1717 safRMin = rho -
fRMin ;
1718 safRMax =
fRMax - rho ;
1724 safZLow = std::fabs((p+vZ).dot(
fLowNorm));
1728 safZHigh = std::fabs((p-vZ).dot(
fHighNorm));
1731 if ( safRMin < safe ) { safe = safRMin; }
1732 if ( safRMax< safe ) { safe = safRMax; }
1746 if (safePhi < safe) { safe = safePhi ; }
1748 if ( safe < 0 ) { safe = 0; }
1767 G4double meshAngle, meshRMax, crossAngle,
1768 cosCrossAngle, sinCrossAngle, sAngle;
1769 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1770 G4int crossSection, noCrossSections;
1786 meshAngle =
fDPhi/(noCrossSections - 1) ;
1796 else { sAngle =
fSPhi ; }
1802 vertices->reserve(noCrossSections*4);
1803 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1807 crossAngle = sAngle + crossSection*meshAngle ;
1808 cosCrossAngle = std::cos(crossAngle) ;
1809 sinCrossAngle = std::sin(crossAngle) ;
1811 rMaxX = meshRMax*cosCrossAngle ;
1812 rMaxY = meshRMax*sinCrossAngle ;
1821 rMinX = meshRMin*cosCrossAngle ;
1822 rMinY = meshRMin*sinCrossAngle ;
1840 "Error in allocation of vertices. Out of memory !");
1869 G4int oldprc = os.precision(16);
1870 os <<
"-----------------------------------------------------------\n"
1871 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1872 <<
" ===================================================\n"
1873 <<
" Solid type: G4CutTubs\n"
1874 <<
" Parameters: \n"
1875 <<
" inner radius : " <<
fRMin/
mm <<
" mm \n"
1876 <<
" outer radius : " <<
fRMax/
mm <<
" mm \n"
1877 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1878 <<
" starting phi : " <<
fSPhi/
degree <<
" degrees \n"
1879 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1880 <<
" low Norm : " <<
fLowNorm <<
" \n"
1882 <<
"-----------------------------------------------------------\n";
1883 os.precision(oldprc);
1894 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1895 aOne, aTwo, aThr, aFou;
1904 cosphi = std::cos(phi);
1905 sinphi = std::sin(phi);
1909 if( (
fSPhi == 0) && (
fDPhi == twopi) ) { aFou = 0; }
1913 if( (chose >=0) && (chose < aOne) )
1915 xRand = fRMax*cosphi;
1916 yRand = fRMax*sinphi;
1921 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1923 xRand = fRMin*cosphi;
1924 yRand = fRMin*sinphi;
1929 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1931 xRand = rRand*cosphi;
1932 yRand = rRand*sinphi;
1936 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1938 xRand = rRand*cosphi;
1939 yRand = rRand*sinphi;
1943 else if( (chose >= aOne + aTwo + 2.*aThr)
1944 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1946 xRand = rRand*std::cos(
fSPhi);
1947 yRand = rRand*std::sin(
fSPhi);
1974 typedef G4int G4int4[4];
1978 G4int nn=ph1->GetNoVertices();
1979 G4int nf=ph1->GetNoFacets();
1980 G4double3* xyz =
new G4double3[
nn];
1981 G4int4* faces =
new G4int4[nf] ;
1985 xyz[i][0]=ph1->GetVertex(i+1).x();
1986 xyz[i][1]=ph1->GetVertex(i+1).y();
1987 G4double tmpZ=ph1->GetVertex(i+1).z();
2004 for(
G4int i=0;i<nf;i++)
2006 ph1->GetFacet(i+1,n,iNodes,iEdge);
2009 faces[i][k]=iNodes[k];
2011 for(
G4int k=n;k<4;k++)
2016 ph->createPolyhedron(nn,nf,xyz,faces);
2033 G4double zXLow1,zXLow2,zYLow1,zYLow2;
2034 G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
2044 if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
2045 || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) {
return true; }
2086 G4bool in_range_low =
false;
2087 G4bool in_range_hi =
false;
2092 if (phiLow<0) { phiLow+=twopi; }
2094 if (ddp<0) { ddp += twopi; }
2097 xc =
fRMin*std::cos(phiLow);
2098 yc =
fRMin*std::sin(phiLow);
2100 xc =
fRMax*std::cos(phiLow);
2101 yc =
fRMax*std::sin(phiLow);
2103 if (in_range_low) { zmin =
std::min(zmin, z1); }
2105 in_range_low =
true;
2108 if (phiLow>twopi) { phiLow-=twopi; }
2112 if (phiHigh<0) { phiHigh+=twopi; }
2114 if (ddp<0) { ddp += twopi; }
2117 xc =
fRMin*std::cos(phiHigh);
2118 yc =
fRMin*std::sin(phiHigh);
2120 xc =
fRMax*std::cos(phiHigh);
2121 yc =
fRMax*std::sin(phiHigh);
2123 if (in_range_hi) { zmax =
std::min(zmax, z1); }
2128 if (phiLow>twopi) { phiHigh-=twopi; }
2154 for (i = 1; i < 4; i++)
2156 if(z[i] < z[i-1])z1=z[i];
2168 for (i = 1; i < 4; i++)
2170 if(z[4+i] > z[4+i-1]) { z1=z[4+i]; }
2173 if (in_range_hi) { zmax =
std::max(zmax, z1); }
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
ThreeVector shoot(const G4int Ap, const G4int Af)
static const G4double kInfinity
G4double GetMinYExtent() const
CLHEP::Hep3Vector G4ThreeVector
G4bool IsYLimited() const
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
void GetMaxMinZ(G4double &zmin, G4double &zmax) const
G4bool IsXLimited() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
const G4int kMinMeshSections
virtual void AddSolid(const G4Box &)=0
EInside Inside(const G4ThreeVector &p) const
G4double GetMaxXExtent() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4CutTubs & operator=(const G4CutTubs &rhs)
G4OTubs & operator=(const G4OTubs &rhs)
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
G4Polyhedron * fpPolyhedron
G4double GetRadialTolerance() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
virtual G4Polyhedron * GetPolyhedron() const
std::ostream & StreamInfo(std::ostream &os) const
std::vector< G4ThreeVector > G4ThreeVectorList
G4double halfRadTolerance
G4Polyhedron * CreatePolyhedron() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double GetMinXExtent() const
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
G4bool IsCrossingCutPlanes() const
G4double GetMaxZExtent() const
G4GeometryType GetEntityType() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4double halfCarTolerance
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const double degree
G4double GetMaxYExtent() const
G4double halfAngTolerance
G4Polyhedron * CreatePolyhedron() const
G4ThreeVector GetPointOnSurface() const
G4double GetMaxExtent(const EAxis pAxis) const
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
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
const G4double kMeshAngleDefault
static G4GeometryTolerance * GetInstance()
G4CutTubs(const G4String &pName, G4double pRMin, G4double pRMax, G4double pDz, G4double pSPhi, G4double pDPhi, G4ThreeVector pLowNorm, G4ThreeVector pHighNorm)
G4double GetCutZ(const G4ThreeVector &p) const
const G4int kMaxMeshSections