62 #if !(defined(G4GEOM_USE_UTORUS) && defined(G4GEOM_USE_SYS_USOLIDS))
80 using namespace CLHEP;
128 message <<
"Invalid swept radius for Solid: " <<
GetName() <<
G4endl
129 <<
" pRtor = " << pRtor <<
", pRmax = " << pRmax;
136 if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
138 if (pRmin >= 1.e2*kCarTolerance) {
fRmin = pRmin ; }
139 else {
fRmin = 0.0 ; }
145 message <<
"Invalid values of radii for Solid: " <<
GetName() <<
G4endl
146 <<
" pRmin = " << pRmin <<
", pRmax = " << pRmax;
162 if (pDPhi > 0) {
fDPhi = pDPhi ; }
166 message <<
"Invalid Z delta-Phi for Solid: " <<
GetName() <<
G4endl
167 <<
" pDPhi = " << pDPhi;
189 :
G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
190 fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
191 kRadTolerance(0.), kAngTolerance(0.),
192 halfCarTolerance(0.), halfAngTolerance(0.)
208 :
G4CSGSolid(rhs), fRmin(rhs.fRmin),fRmax(rhs.fRmax),
209 fRtor(rhs.fRtor),fSPhi(rhs.fSPhi),fDPhi(rhs.fDPhi),
210 fRminTolerance(rhs.fRminTolerance), fRmaxTolerance(rhs.fRmaxTolerance),
211 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
212 halfCarTolerance(rhs.halfCarTolerance),
213 halfAngTolerance(rhs.halfAngTolerance)
225 if (
this == &rhs) {
return *
this; }
265 std::vector<G4double>& roots )
const
279 c[2] = 2*( (d + 2*pDotV*pDotV - r2) + 2*Rtor2*v.
z()*v.
z());
280 c[3] = 4*(pDotV*(d - r2) + 2*Rtor2*p.
z()*v.
z()) ;
281 c[4] = (d-r2)*(d-r2) +4*Rtor2*(p.
z()*p.
z()-r2);
285 num = torusEq.
FindRoots( c, 4, srd, si );
287 for ( i = 0; i < num; i++ )
289 if( si[i] == 0. ) { roots.push_back(srd[i]) ; }
292 std::sort(roots.begin() , roots.end() ) ;
305 G4bool IsDistanceToIn )
const
314 std::vector<G4double> roots ;
315 std::vector<G4double> rootsrefined ;
322 for (
size_t k = 0 ; k<roots.size() ; k++ )
332 if ( rootsrefined.size()==roots.size() )
334 t = t + rootsrefined[k] ;
340 G4double theta = std::atan2(ptmp.
y(),ptmp.
x());
362 if ( IsDistanceToIn ==
true )
370 p.
y()*(1-
fRtor/std::hypot(p.
x(),p.
y())),
375 if ( r ==
GetRmin() ) { scal = -scal ; }
376 if ( scal < 0 ) {
return 0.0 ; }
383 if ( IsDistanceToIn ==
false )
390 p.
y()*(1-
fRtor/std::hypot(p.
x(),p.
y())),
395 if ( r ==
GetRmin() ) { scal = -scal ; }
396 if ( scal > 0 ) {
return 0.0 ; }
429 pMin.
set(-rext,-rext,-dz);
430 pMax.
set( rext, rext, dz);
439 pMin.
set(vmin.
x(),vmin.
y(),-dz);
440 pMax.
set(vmax.
x(),vmax.
y(), dz);
445 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
448 message <<
"Bad bounding box (min >= max) for solid: "
450 <<
"\npMin = " << pMin
451 <<
"\npMax = " << pMax;
452 G4Exception(
"G4Torus::BoundingLimits()",
"GeomMgt0001",
476 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
480 return exist = (pMin < pMax) ?
true :
false;
497 static const G4int NPHI = 24;
498 static const G4int NDISK = 16;
499 static const G4double sinHalfDisk = std::sin(
pi/NDISK);
500 static const G4double cosHalfDisk = std::cos(
pi/NDISK);
501 static const G4double sinStepDisk = 2.*sinHalfDisk*cosHalfDisk;
502 static const G4double cosStepDisk = 1. - 2.*sinHalfDisk*sinHalfDisk;
505 G4int kphi = (dphi <= astep) ? 1 : (
G4int)((dphi-
deg)/astep) + 1;
508 G4double sinHalf = std::sin(0.5*ang);
509 G4double cosHalf = std::cos(0.5*ang);
510 G4double sinStep = 2.*sinHalf*cosHalf;
511 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
515 for (
G4int k=0; k<NDISK+1; ++k) pols[k].resize(4);
517 std::vector<const G4ThreeVectorList *> polygons;
518 polygons.resize(NDISK+1);
519 for (
G4int k=0; k<NDISK+1; ++k) polygons[k] = &pols[k];
525 if ((rtor-rmin*sinHalfDisk)/cosHalf > (rtor+rmin*sinHalfDisk)) rmin = 0;
529 for (
G4int k=0; k<NDISK; ++k)
531 G4double rmincur = rtor + rmin*cosCurDisk;
532 if (cosCurDisk < 0 && rmin > 0) rmincur /= cosHalf;
533 rzmin[k].
set(rmincur,rmin*sinCurDisk);
535 G4double rmaxcur = rtor + rmax*cosCurDisk;
536 if (cosCurDisk > 0) rmaxcur /= cosHalf;
537 rzmax[k].
set(rmaxcur,rmax*sinCurDisk);
540 sinCurDisk = sinCurDisk*cosStepDisk + cosCurDisk*sinStepDisk;
541 cosCurDisk = cosCurDisk*cosStepDisk - sinTmpDisk*sinStepDisk;
550 G4double sinCur1 = 0, cosCur1 = 0, sinCur2 = 0, cosCur2 = 0;
551 for (
G4int i=0; i<kphi+1; ++i)
557 sinCur2 = sinCur1*cosHalf + cosCur1*sinHalf;
558 cosCur2 = cosCur1*cosHalf - sinCur1*sinHalf;
564 sinCur2 = (i == kphi) ? sinEnd : sinCur1*cosStep + cosCur1*sinStep;
565 cosCur2 = (i == kphi) ? cosEnd : cosCur1*cosStep - sinCur1*sinStep;
567 for (
G4int k=0; k<NDISK; ++k)
569 G4double r1 = rzmin[k].
x(), r2 = rzmax[k].
x();
570 G4double z1 = rzmin[k].
y(), z2 = rzmax[k].
y();
571 pols[k][0].set(r1*cosCur1,r1*sinCur1,z1);
572 pols[k][1].set(r2*cosCur1,r2*sinCur1,z2);
573 pols[k][2].set(r2*cosCur2,r2*sinCur2,z2);
574 pols[k][3].set(r1*cosCur2,r1*sinCur2,z1);
576 pols[NDISK] = pols[0];
581 DiskExtent(rint,rext,sinCur1,cosCur1,sinCur2,cosCur2,vmin,vmax);
588 if (!benv.
CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax))
continue;
589 if (emin < pMin) pMin = emin;
590 if (emax > pMax) pMax =
emax;
591 if (eminlim > pMin && emaxlim < pMax)
break;
593 return (pMin < pMax);
602 G4double r, pt2, pPhi, tolRMin, tolRMax ;
608 r = std::hypot(p.
x(),p.
y());
616 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
627 pPhi = std::atan2(p.
y(),p.
x()) ;
664 if (tolRMin < 0 ) { tolRMin = 0 ; }
666 if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
674 pPhi = std::atan2(p.
y(),p.
x()) ;
713 G4int noSurfaces = 0;
727 rho = std::hypot(p.
x(),p.
y());
728 pt = std::hypot(p.
z(),rho-
fRtor);
733 if( rho > delta && pt != 0.0 )
746 pPhi = std::atan2(p.
y(),p.
x());
751 distSPhi = std::fabs( pPhi -
fSPhi );
757 if( distRMax <= delta )
762 else if(
fRmin && (distRMin <= delta) )
773 if (distSPhi <= dAngle)
778 if (distEPhi <= dAngle)
784 if ( noSurfaces == 0 )
794 ed <<
" ERROR> Surface Normal was called for Torus,"
795 <<
" with point not on surface." <<
G4endl;
799 ed <<
" ERROR> Surface Normal has not found a surface, "
800 <<
" despite the point being on the surface. " <<
G4endl;
811 ed <<
" Coordinates of point : " << p <<
G4endl;
812 ed <<
" Parameters of solid : " << G4endl << *
this <<
G4endl;
816 G4Exception(
"G4Torus::SurfaceNormal(p)",
"GeomSolids1002",
818 "Failing to find normal, even though point is on surface!");
822 static const char* NameInside[3]= {
"Inside",
"Surface",
"Outside" };
823 ed <<
" The point is " << NameInside[inIt] <<
" the solid. "<<
G4endl;
824 G4Exception(
"G4Torus::SurfaceNormal(p)",
"GeomSolids1002",
830 else if ( noSurfaces == 1 ) { norm = sumnorm; }
831 else { norm = sumnorm.
unit(); }
846 G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
848 rho = std::hypot(p.
x(),p.
y());
849 pt = std::hypot(p.
z(),rho-
fRtor);
852 G4cout <<
" G4Torus::ApproximateSurfaceNormal called for point " << p
856 distRMax = std::fabs(pt -
fRmax) ;
860 distRMin = std::fabs(pt -
fRmin) ;
862 if (distRMin < distRMax)
880 phi = std::atan2(p.
y(),p.
x()) ;
882 if (phi < 0) { phi +=
twopi ; }
885 else { distSPhi = std::fabs(phi-
fSPhi)*rho ; }
887 distEPhi = std::fabs(phi -
fSPhi -
fDPhi)*rho ;
889 if (distSPhi < distEPhi)
891 if (distSPhi<distMin) side =
kNSPhi ;
895 if (distEPhi < distMin) { side =
kNEPhi ; }
920 "Undefined side for valid surface normal to solid.");
961 G4double cPhi,sinCPhi=0.,cosCPhi=0.;
978 cPhi =
fSPhi + hDPhi ;
979 sinCPhi = std::sin(cPhi) ;
980 cosCPhi = std::cos(cPhi) ;
1004 if ( sd[0] < snxt ) { snxt = sd[0] ; }
1019 sinSPhi = std::sin(
fSPhi) ;
1020 cosSPhi = std::cos(
fSPhi) ;
1021 Comp = v.
x()*sinSPhi - v.
y()*cosSPhi ;
1025 Dist = (p.
y()*cosSPhi - p.
x()*sinSPhi) ;
1032 if ( sphi < 0 ) { sphi = 0 ; }
1034 xi = p.
x() + sphi*v.
x() ;
1035 yi = p.
y() + sphi*v.
y() ;
1036 zi = p.
z() + sphi*v.
z() ;
1037 rhoi = std::hypot(xi,yi);
1040 if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
1045 if ((yi*cosCPhi-xi*sinCPhi)<=0) { snxt=sphi; }
1051 sinEPhi=std::sin(ePhi);
1052 cosEPhi=std::cos(ePhi);
1053 Comp=-(v.
x()*sinEPhi-v.
y()*cosEPhi);
1057 Dist = -(p.
y()*cosEPhi - p.
x()*sinEPhi) ;
1065 if (sphi < 0 ) { sphi = 0 ; }
1067 xi = p.
x() + sphi*v.
x() ;
1068 yi = p.
y() + sphi*v.
y() ;
1069 zi = p.
z() + sphi*v.
z() ;
1070 rhoi = std::hypot(xi,yi);
1073 if (it2 >= tolORMin2 && it2 <= tolORMax2)
1078 if ((yi*cosCPhi-xi*sinCPhi)>=0) { snxt=sphi; }
1099 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1102 rho = std::hypot(p.
x(),p.
y());
1103 pt = std::hypot(p.
z(),rho-
fRtor);
1105 safe2 = pt -
fRmax ;
1107 if (safe1 > safe2) { safe = safe1; }
1108 else { safe = safe2; }
1113 cosPhiC = std::cos(phiC) ;
1114 sinPhiC = std::sin(phiC) ;
1115 cosPsi = (p.
x()*cosPhiC + p.
y()*sinPhiC)/rho ;
1117 if (cosPsi < std::cos(
fDPhi*0.5) )
1119 if ((p.
y()*cosPhiC - p.
x()*sinPhiC) <= 0 )
1121 safePhi = std::fabs(p.
x()*std::sin(
fSPhi) - p.
y()*std::cos(
fSPhi)) ;
1126 safePhi = std::fabs(p.
x()*std::sin(ePhi) - p.
y()*std::cos(ePhi)) ;
1128 if (safePhi > safe) { safe = safePhi ; }
1131 if (safe < 0 ) { safe = 0 ; }
1152 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
1154 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
1177 if( (pt*pt > tolRMax*tolRMax) && (vDotNmax >= 0) )
1183 if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) )
1186 p.
y()*(1 -
fRtor/rho)/pt,
1203 if ( (pt*pt < tolRMin*tolRMin) && (vDotNmax < 0) )
1205 if (calcNorm) { *validNorm = false ; }
1235 if ( calcNorm && (snxt == 0.0) )
1237 *validNorm = false ;
1245 sinSPhi = std::sin(
fSPhi) ;
1246 cosSPhi = std::cos(
fSPhi) ;
1248 sinEPhi = std::sin(ePhi) ;
1249 cosEPhi = std::cos(ePhi) ;
1250 cPhi =
fSPhi + fDPhi*0.5 ;
1251 sinCPhi = std::sin(cPhi) ;
1252 cosCPhi = std::cos(cPhi) ;
1257 vphi = std::atan2(v.
y(),v.
x()) ;
1262 if ( p.
x() || p.
y() )
1264 pDistS = p.
x()*sinSPhi - p.
y()*cosSPhi ;
1265 pDistE = -p.
x()*sinEPhi + p.
y()*cosEPhi ;
1269 compS = -sinSPhi*v.
x() + cosSPhi*v.
y() ;
1270 compE = sinEPhi*v.
x() - cosEPhi*v.
y() ;
1282 sphi = pDistS/compS ;
1286 xi = p.
x() + sphi*v.
x() ;
1287 yi = p.
y() + sphi*v.
y() ;
1302 else if ( yi*cosCPhi-xi*sinCPhi >=0 )
1323 sphi2 = pDistE/compE ;
1329 xi = p.
x() + sphi2*v.
x() ;
1330 yi = p.
y() + sphi2*v.
y() ;
1346 if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1368 vphi = std::atan2(v.
y(),v.
x());
1400 xi = p.
x() + snxt*v.
x() ;
1401 yi = p.
y() + snxt*v.
y() ;
1402 zi = p.
z() + snxt*v.
z() ;
1403 rhoi = std::hypot(xi,yi);
1404 it = hypot(zi,rhoi-
fRtor);
1406 iDotxyNmax = (1-
fRtor/rhoi) ;
1407 if(iDotxyNmax >= -2.*fRmaxTolerance)
1410 yi*(1-
fRtor/rhoi)/it,
1416 *validNorm = false ;
1421 *validNorm = false ;
1432 *validNorm = false ;
1444 *validNorm = false ;
1455 G4int oldprc = message.precision(16);
1456 message <<
"Undefined side for valid surface normal to solid."
1458 <<
"Position:" << G4endl << G4endl
1459 <<
"p.x() = " << p.
x()/
mm <<
" mm" << G4endl
1460 <<
"p.y() = " << p.
y()/
mm <<
" mm" << G4endl
1461 <<
"p.z() = " << p.
z()/
mm <<
" mm" << G4endl << G4endl
1462 <<
"Direction:" << G4endl << G4endl
1463 <<
"v.x() = " << v.
x() << G4endl
1464 <<
"v.y() = " << v.
y() << G4endl
1465 <<
"v.z() = " << v.
z() << G4endl << G4endl
1466 <<
"Proposed distance :" << G4endl << G4endl
1467 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl;
1468 message.precision(oldprc);
1487 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
1489 rho = std::hypot(p.
x(),p.
y());
1490 pt = std::hypot(p.
z(),rho-
fRtor);
1502 G4cout.precision(oldprc);
1503 G4Exception(
"G4Torus::DistanceToOut(p)",
"GeomSolids1002",
1510 safeR1 = pt -
fRmin ;
1513 if (safeR1 < safeR2) { safe = safeR1 ; }
1514 else { safe = safeR2 ; }
1526 cosPhiC = std::cos(phiC) ;
1527 sinPhiC = std::sin(phiC) ;
1529 if ((p.
y()*cosPhiC-p.
x()*sinPhiC)<=0)
1531 safePhi = -(p.
x()*std::sin(
fSPhi) - p.
y()*std::cos(
fSPhi)) ;
1536 safePhi = (p.
x()*std::sin(ePhi) - p.
y()*std::cos(ePhi)) ;
1538 if (safePhi < safe) { safe = safePhi ; }
1540 if (safe < 0) { safe = 0 ; }
1568 G4int oldprc = os.precision(16);
1569 os <<
"-----------------------------------------------------------\n"
1570 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1571 <<
" ===================================================\n"
1572 <<
" Solid type: G4Torus\n"
1573 <<
" Parameters: \n"
1574 <<
" inner radius: " <<
fRmin/
mm <<
" mm \n"
1575 <<
" outer radius: " <<
fRmax/
mm <<
" mm \n"
1576 <<
" swept radius: " <<
fRtor/
mm <<
" mm \n"
1577 <<
" starting phi: " <<
fSPhi/
degree <<
" degrees \n"
1578 <<
" delta phi : " <<
fDPhi/
degree <<
" degrees \n"
1579 <<
"-----------------------------------------------------------\n";
1580 os.precision(oldprc);
1591 G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
1596 cosu = std::cos(phi); sinu = std::sin(phi);
1597 cosv = std::cos(theta); sinv = std::sin(theta);
1613 else if( (chose >= aOut) && (chose < aOut + aIn) )
1618 else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
1622 (
fRtor+rRand*cosv)*std::sin(
fSPhi), rRand*sinv);
1647 #endif // !defined(G4GEOM_USE_TORUS) || !defined(G4GEOM_USE_SYS_USOLIDS)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void SetAllParameters(G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
void set(double x, double y, double z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetCosEndPhi() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetMaxExtent(const EAxis pAxis) const
std::ostringstream G4ExceptionDescription
std::vector< ExP01TrackerHit * > a
void TorusRootsJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, std::vector< G4double > &roots) const
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
static constexpr double mm
EInside Inside(const G4ThreeVector &p) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4GeometryType GetEntityType() const
void message(RunManager *runmanager)
G4double GetRadiusInRing(G4double rmin, G4double rmax) const
G4double GetAngularTolerance() const
G4double GetRadialTolerance() const
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
static const G4double emax
G4double halfCarTolerance
G4Polyhedron * CreatePolyhedron() const
G4double SolveNumericJT(const G4ThreeVector &p, const G4ThreeVector &v, G4double r, G4bool IsDistanceToIn) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
std::vector< G4ThreeVector > G4ThreeVectorList
G4double halfAngTolerance
G4Torus(const G4String &pName, G4double pRmin, G4double pRmax, G4double pRtor, G4double pSPhi, G4double pDPhi)
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
G4double GetMinExtent(const EAxis pAxis) const
static constexpr double deg
static constexpr double twopi
G4ThreeVector GetPointOnSurface() const
std::ostream & StreamInfo(std::ostream &os) const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4double GetSinStartPhi() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double GetSinEndPhi() const
void set(double x, double y)
void DescribeYourselfTo(G4VGraphicsScene &scene) const
static constexpr double degree
G4GLOB_DLL std::ostream G4cout
static G4GeometryTolerance * GetInstance()
static constexpr double pi
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=G4bool(false), G4bool *validNorm=0, G4ThreeVector *n=0) const
G4bool fRebuildPolyhedron
G4double GetCosStartPhi() const
G4Torus & operator=(const G4Torus &rhs)