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)
177 if (
this == &rhs) {
return *
this; }
215 G4double diff1, diff2, maxDiff, newMin, newMax;
216 G4double xoff1, xoff2, yoff1, yoff2, delta;
219 xMin = xoffset -
fRMax;
220 xMax = xoffset +
fRMax;
242 yMin = yoffset -
fRMax;
243 yMax = yoffset +
fRMax;
292 yoff1 = yoffset - yMin;
293 yoff2 = yMax - yoffset;
295 if ( (yoff1 >= 0) && (yoff2 >= 0) )
305 delta = fRMax*fRMax - yoff1*yoff1;
306 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
307 delta = fRMax*fRMax - yoff2*yoff2;
308 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
309 maxDiff = (diff1 > diff2) ? diff1:diff2;
310 newMin = xoffset - maxDiff;
311 newMax = xoffset + maxDiff;
312 pMin = (newMin < xMin) ? xMin : newMin;
313 pMax = (newMax > xMax) ? xMax : newMax;
319 xoff1 = xoffset - xMin;
320 xoff2 = xMax - xoffset;
322 if ( (xoff1 >= 0) && (xoff2 >= 0) )
332 delta = fRMax*fRMax - xoff1*xoff1;
333 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
334 delta = fRMax*fRMax - xoff2*xoff2;
335 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
336 maxDiff = (diff1 > diff2) ? diff1 : diff2;
337 newMin = yoffset - maxDiff;
338 newMax = yoffset + maxDiff;
339 pMin = (newMin < yMin) ? yMin : newMin;
340 pMax = (newMax > yMax) ? yMax : newMax;
359 G4int i, noEntries, noBetweenSections4;
360 G4bool existsAfterClip =
false;
366 noEntries = vertices->size();
367 noBetweenSections4 = noEntries - 4;
369 for ( i = 0 ; i < noEntries ; i += 4 )
373 for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
379 existsAfterClip =
true;
397 existsAfterClip =
true;
403 return existsAfterClip;
432 r2 = p.x()*p.x() + p.y()*p.y() ;
438 if ( tolRMin < 0 ) { tolRMin = 0; }
440 if ( ((r2 < tolRMin*tolRMin) || (r2 > tolRMax*tolRMax))
441 && (r2 >=halfRadTolerance*halfRadTolerance) ) {
return kOutside; }
455 pPhi = std::atan2(p.y(),p.x()) ;
501 else { tolRMin = 0 ; }
503 if ( ((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax))&&
520 G4int noSurfaces = 0;
522 G4double distZLow,distZHigh, distRMin, distRMax;
530 rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
532 distRMin = std::fabs(rho -
fRMin);
533 distRMax = std::fabs(rho -
fRMax);
537 distZLow =std::fabs((p+vZ).dot(
fLowNorm));
541 distZHigh = std::fabs((p-vZ).dot(
fHighNorm));
547 pPhi = std::atan2(p.y(),p.x());
552 distSPhi = std::fabs(pPhi -
fSPhi);
598 if ( noSurfaces == 0 )
601 G4Exception(
"G4CutTubs::SurfaceNormal(p)",
"GeomSolids1002",
604 G4cout<<
"G4CutTubs::SN ( "<<p.x()<<
", "<<p.y()<<
", "<<p.z()<<
" ); "
606 G4cout.precision(oldprc) ;
610 else if ( noSurfaces == 1 ) { norm = sumnorm; }
611 else { norm = sumnorm.unit(); }
627 G4double 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) ;
637 distZLow =std::fabs((p+vZ).dot(
fLowNorm));
641 distZHigh = std::fabs((p-vZ).dot(
fHighNorm));
644 if (distRMin < distRMax)
646 if ( distZ < distRMin )
659 if ( distZ < distRMax )
672 phi = std::atan2(p.y(),p.x()) ;
674 if ( phi < 0 ) { phi +=
twopi; }
678 distSPhi = std::fabs(phi - (
fSPhi +
twopi))*rho ;
682 distSPhi = std::fabs(phi -
fSPhi)*rho ;
684 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
686 if (distSPhi < distEPhi)
688 if ( distSPhi < distMin )
695 if ( distEPhi < distMin )
715 if ( distZHigh > distZLow ) { norm =
fHighNorm ; }
734 "Undefined side for valid surface normal to solid.");
774 G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ;
808 if(sd < 0.0) { sd = 0.0; }
810 xi = p.x() + sd*v.x() ;
811 yi = p.y() + sd*v.y() ;
812 rho2 = xi*xi + yi*yi ;
816 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
823 iden = std::sqrt(rho2) ;
848 sd = -distZHigh/calf;
850 if(sd < 0.0) { sd = 0.0; }
852 xi = p.x() + sd*v.x() ;
853 yi = p.y() + sd*v.y() ;
854 rho2 = xi*xi + yi*yi ;
858 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
865 iden = std::sqrt(rho2) ;
896 t1 = 1.0 - v.z()*v.z() ;
897 t2 = p.x()*v.x() + p.y()*v.y() ;
898 t3 = p.x()*p.x() + p.y()*p.y() ;
904 if ((t3 >= tolORMax2) && (t2<0))
913 sd = c/(-b+std::sqrt(d));
918 G4double fTerm = sd-std::fmod(sd,dRmax);
923 zi = p.z() + sd*v.z() ;
924 xi = p.x() + sd*v.x() ;
925 yi = p.y() + sd*v.y() ;
940 xi = p.x() + sd*v.x() ;
941 yi = p.y() + sd*v.y() ;
954 if ((t3 > tolIRMin2) && (t2 < 0)
962 iden = std::sqrt(t3) ;
983 snxt = c/(-b+std::sqrt(d));
1003 c = t3 - fRMax*
fRMax;
1014 snxt= c/(-b+std::sqrt(d));
1037 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d));
1042 if (sd < 0.0) { sd = 0.0; }
1045 G4double fTerm = sd-std::fmod(sd,dRmax);
1048 zi = p.z() + sd*v.z() ;
1049 xi = p.x() + sd*v.x() ;
1050 yi = p.y() + sd*v.y() ;
1106 if ( sd < 0 ) { sd = 0.0; }
1107 zi = p.z() + sd*v.z() ;
1108 xi = p.x() + sd*v.x() ;
1109 yi = p.y() + sd*v.y() ;
1116 rho2 = xi*xi + yi*yi ;
1117 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1118 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1121 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1150 if ( sd < 0 ) { sd = 0; }
1151 zi = p.z() + sd*v.z() ;
1152 xi = p.x() + sd*v.x() ;
1153 yi = p.y() + sd*v.y() ;
1160 xi = p.x() + sd*v.x() ;
1161 yi = p.y() + sd*v.y() ;
1162 rho2 = xi*xi + yi*yi ;
1163 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1164 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2)
1167 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1218 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi;
1223 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1225 safRMin =
fRMin- rho ;
1226 safRMax = rho -
fRMax ;
1240 if ( safRMin > safe ) { safe = safRMin; }
1241 if ( safRMax> safe ) { safe = safRMax; }
1251 if ( cosPsi < std::cos(
fDPhi*0.5) )
1263 if ( safePhi > safe ) { safe = safePhi; }
1266 if ( safe < 0 ) { safe = 0; }
1284 G4double deltaR, t1, t2, t3, b, c,
d2, roMin2 ;
1285 G4double distZLow,distZHigh,calfH,calfL;
1290 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1310 snxt = -distZHigh/calfH ;
1328 sz = -distZLow/calfL ;
1345 if((calfH<=0)&&(calfL<=0))
1361 t1 = 1.0 - v.z()*v.z() ;
1362 t2 = p.x()*v.x() + p.y()*v.y() ;
1363 t3 = p.x()*p.x() + p.y()*p.y() ;
1366 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }
1386 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); }
1405 roMin2 = t3 - t2*t2/t1 ;
1421 srd = c/(-b+std::sqrt(d2));
1426 if ( calcNorm ) { *validNorm =
false; }
1437 srd = -b + std::sqrt(d2) ;
1461 srd = -b + std::sqrt(d2) ;
1483 vphi = std::atan2(v.y(),v.x()) ;
1489 if ( p.x() || p.y() )
1512 sphi = pDistS/compS ;
1516 xi = p.x() + sphi*v.x() ;
1517 yi = p.y() + sphi*v.y() ;
1557 sphi2 = pDistE/compE ;
1563 xi = p.x() + sphi2*v.x() ;
1564 yi = p.y() + sphi2*v.y() ;
1575 else { sphi = 0.0 ; }
1586 else { sphi = 0.0 ; }
1632 xi = p.x() + snxt*v.x() ;
1633 yi = p.y() + snxt*v.y() ;
1639 *validNorm = false ;
1650 *validNorm = false ;
1662 *validNorm = false ;
1679 std::ostringstream message;
1680 G4int oldprc = message.precision(16);
1681 message <<
"Undefined side for valid surface normal to solid."
1683 <<
"Position:" << G4endl << G4endl
1684 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1685 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1686 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1687 <<
"Direction:" << G4endl << G4endl
1688 <<
"v.x() = " << v.x() << G4endl
1689 <<
"v.y() = " << v.y() << G4endl
1690 <<
"v.z() = " << v.z() << G4endl << G4endl
1691 <<
"Proposed distance :" << G4endl << G4endl
1692 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl ;
1693 message.precision(oldprc) ;
1694 G4Exception(
"G4CutTubs::DistanceToOut(p,v,..)",
"GeomSolids1002",
1709 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho;
1712 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1714 safRMin = rho -
fRMin ;
1715 safRMax =
fRMax - rho ;
1721 safZLow = std::fabs((p+vZ).dot(
fLowNorm));
1725 safZHigh = std::fabs((p-vZ).dot(
fHighNorm));
1728 if ( safRMin < safe ) { safe = safRMin; }
1729 if ( safRMax< safe ) { safe = safRMax; }
1743 if (safePhi < safe) { safe = safePhi ; }
1745 if ( safe < 0 ) { safe = 0; }
1764 G4double meshAngle, meshRMax, crossAngle,
1765 cosCrossAngle, sinCrossAngle, sAngle;
1766 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1767 G4int crossSection, noCrossSections;
1783 meshAngle =
fDPhi/(noCrossSections - 1) ;
1793 else { sAngle =
fSPhi ; }
1799 vertices->reserve(noCrossSections*4);
1800 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1804 crossAngle = sAngle + crossSection*meshAngle ;
1805 cosCrossAngle = std::cos(crossAngle) ;
1806 sinCrossAngle = std::sin(crossAngle) ;
1808 rMaxX = meshRMax*cosCrossAngle ;
1809 rMaxY = meshRMax*sinCrossAngle ;
1818 rMinX = meshRMin*cosCrossAngle ;
1819 rMinY = meshRMin*sinCrossAngle ;
1837 "Error in allocation of vertices. Out of memory !");
1866 G4int oldprc = os.precision(16);
1867 os <<
"-----------------------------------------------------------\n"
1868 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1869 <<
" ===================================================\n"
1870 <<
" Solid type: G4CutTubs\n"
1871 <<
" Parameters: \n"
1872 <<
" inner radius : " <<
fRMin/
mm <<
" mm \n"
1873 <<
" outer radius : " <<
fRMax/
mm <<
" mm \n"
1874 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1875 <<
" starting phi : " <<
fSPhi/
degree <<
" degrees \n"
1876 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1877 <<
" low Norm : " <<
fLowNorm <<
" \n"
1879 <<
"-----------------------------------------------------------\n";
1880 os.precision(oldprc);
1891 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1892 aOne, aTwo, aThr, aFou;
1901 cosphi = std::cos(phi);
1902 sinphi = std::sin(phi);
1910 if( (chose >=0) && (chose < aOne) )
1912 xRand = fRMax*cosphi;
1913 yRand = fRMax*sinphi;
1918 else if( (chose >= aOne) && (chose < aOne + aTwo) )
1920 xRand = fRMin*cosphi;
1921 yRand = fRMin*sinphi;
1926 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1928 xRand = rRand*cosphi;
1929 yRand = rRand*sinphi;
1933 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1935 xRand = rRand*cosphi;
1936 yRand = rRand*sinphi;
1940 else if( (chose >= aOne + aTwo + 2.*aThr)
1941 && (chose < aOne + aTwo + 2.*aThr + aFou) )
1943 xRand = rRand*std::cos(
fSPhi);
1944 yRand = rRand*std::sin(
fSPhi);
1971 typedef G4int G4int4[4];
1975 G4int nn=ph1->GetNoVertices();
1976 G4int nf=ph1->GetNoFacets();
1977 G4double3* xyz =
new G4double3[
nn];
1978 G4int4* faces =
new G4int4[nf] ;
1982 xyz[i][0]=ph1->GetVertex(i+1).x();
1983 xyz[i][1]=ph1->GetVertex(i+1).y();
1984 G4double tmpZ=ph1->GetVertex(i+1).z();
2001 for(
G4int i=0;i<nf;i++)
2003 ph1->GetFacet(i+1,n,iNodes,iEdge);
2006 faces[i][k]=iNodes[k];
2008 for(
G4int k=n;k<4;k++)
2013 ph->createPolyhedron(nn,nf,xyz,faces);
2030 G4double zXLow1,zXLow2,zYLow1,zYLow2;
2031 G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2;
2041 if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2)
2042 || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) {
return true; }
2083 G4bool in_range_low =
false;
2084 G4bool in_range_hi =
false;
2089 if (phiLow<0) { phiLow+=
twopi; }
2091 if (ddp<0) { ddp +=
twopi; }
2094 xc =
fRMin*std::cos(phiLow);
2095 yc =
fRMin*std::sin(phiLow);
2097 xc =
fRMax*std::cos(phiLow);
2098 yc =
fRMax*std::sin(phiLow);
2100 if (in_range_low) { zmin =
std::min(zmin, z1); }
2102 in_range_low =
true;
2109 if (phiHigh<0) { phiHigh+=
twopi; }
2111 if (ddp<0) { ddp +=
twopi; }
2114 xc =
fRMin*std::cos(phiHigh);
2115 yc =
fRMin*std::sin(phiHigh);
2117 xc =
fRMax*std::cos(phiHigh);
2118 yc =
fRMax*std::sin(phiHigh);
2120 if (in_range_hi) { zmax =
std::min(zmax, z1); }
2151 for (i = 1; i < 4; i++)
2153 if(z[i] < z[i-1])z1=z[i];
2165 for (i = 1; i < 4; i++)
2167 if(z[4+i] > z[4+i-1]) { z1=z[4+i]; }
2170 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
G4double GetRadialTolerance() const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
static const double twopi
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