31 double pRMin,
double pRMax,
33 double pSPhi,
double pDPhi)
34 :
VUSolid(pName.c_str()), fRMin(pRMin), fRMax(pRMax), fDz(pDz), fSPhi(0), fDPhi(0)
43 std::ostringstream message;
44 message <<
"Negative Z half-length (" << pDz <<
") in solid: " <<
GetName();
48 if ((pRMin >= pRMax) || (pRMin < 0))
50 std::ostringstream message;
51 message <<
"Invalid values for radii in solid: " <<
GetName()
53 <<
"pRMin = " << pRMin <<
", pRMax = " << pRMax;
69 :
VUSolid(
""), fCubicVolume(0.), fSurfaceArea(0.),
70 kRadTolerance(0.), kAngTolerance(0.),
71 fRMin(0.), fRMax(0.), fDz(0.), fSPhi(0.), fDPhi(0.),
72 fSinCPhi(0.), fCosCPhi(0.), fCosHDPhiOT(0.), fCosHDPhiIT(0.),
73 fSinSPhi(0.), fCosSPhi(0.), fSinEPhi(0.), fCosEPhi(0.),
74 fSinSPhiDPhi(0.), fCosSPhiDPhi(0.),
93 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
94 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance),
95 fRMin(rhs.fRMin), fRMax(rhs.fRMax), fDz(rhs.fDz),
96 fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
97 fSinCPhi(rhs.fSinCPhi), fCosCPhi(rhs.fCosCPhi),
98 fCosHDPhiOT(rhs.fCosHDPhiOT), fCosHDPhiIT(rhs.fCosHDPhiIT),
99 fSinSPhi(rhs.fSinSPhi), fCosSPhi(rhs.fCosSPhi),
100 fSinEPhi(rhs.fSinEPhi), fCosEPhi(rhs.fCosEPhi),
101 fSinSPhiDPhi(rhs.fSinSPhiDPhi), fCosSPhiDPhi(rhs.fCosSPhiDPhi),
102 fPhiFullTube(rhs.fPhiFullTube)
121 VUSolid::operator=(rhs);
157 double r2, pPhi, tolRMin, tolRMax;
163 if (std::fabs(p.
z()) <=
fDz - halfCarTolerance)
165 r2 = p.
x() * p.
x() + p.
y() * p.
y();
169 tolRMin =
fRMin + halfRadTolerance;
176 tolRMax =
fRMax - halfRadTolerance;
178 if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
189 if ((tolRMin == 0) && (std::fabs(p.
x()) <= halfCarTolerance)
190 && (std::fabs(p.
y()) <= halfCarTolerance))
196 pPhi = std::atan2(p.
y(), p.
x());
197 if (pPhi < -halfAngTolerance)
204 if ((std::fabs(pPhi) < halfAngTolerance)
209 if ((pPhi >=
fSPhi + halfAngTolerance)
210 && (pPhi <=
fSPhi +
fDPhi - halfAngTolerance))
214 else if ((pPhi >=
fSPhi - halfAngTolerance)
215 && (pPhi <=
fSPhi +
fDPhi + halfAngTolerance))
223 && (pPhi >=
fSPhi +
fDPhi + halfAngTolerance))
228 && (pPhi >=
fSPhi +
fDPhi - halfAngTolerance))
242 tolRMin =
fRMin - halfRadTolerance;
243 tolRMax =
fRMax + halfRadTolerance;
250 if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
252 if (
fPhiFullTube || (r2 <= halfRadTolerance * halfRadTolerance))
259 pPhi = std::atan2(p.
y(), p.
x());
261 if (pPhi < -halfAngTolerance)
267 if ((std::fabs(pPhi) < halfAngTolerance)
272 if ((pPhi >=
fSPhi - halfAngTolerance)
273 && (pPhi <=
fSPhi +
fDPhi + halfAngTolerance))
281 && (pPhi >=
fSPhi +
fDPhi + halfAngTolerance))
294 else if (std::fabs(p.
z()) <=
fDz + halfCarTolerance)
297 r2 = p.
x() * p.
x() + p.
y() * p.
y();
298 tolRMin =
fRMin - halfRadTolerance;
299 tolRMax =
fRMax + halfRadTolerance;
306 if ((r2 >= tolRMin * tolRMin) && (r2 <= tolRMax * tolRMax))
308 if (
fPhiFullTube || (r2 <= halfRadTolerance * halfRadTolerance))
315 pPhi = std::atan2(p.
y(), p.
x());
317 if (pPhi < -halfAngTolerance)
323 if ((std::fabs(pPhi) < halfAngTolerance)
328 if ((pPhi >=
fSPhi - halfAngTolerance)
329 && (pPhi <=
fSPhi +
fDPhi + halfAngTolerance))
337 && (pPhi >=
fSPhi +
fDPhi + halfAngTolerance))
362 double distZ, distRMin, distRMax;
372 rho = std::sqrt(p.
x() * p.
x() + p.
y() * p.
y());
374 distRMin = std::fabs(rho -
fRMin);
375 distRMax = std::fabs(rho -
fRMax);
376 distZ = std::fabs(std::fabs(p.
z()) -
fDz);
380 if (rho > halfCarTolerance)
382 pPhi = std::atan2(p.
y(), p.
x());
384 if (pPhi <
fSPhi - halfCarTolerance)
388 else if (pPhi >
fSPhi +
fDPhi + halfCarTolerance)
393 distSPhi = std::fabs(pPhi -
fSPhi);
404 if (rho > halfCarTolerance)
409 if (distRMax <= halfCarTolerance)
414 if (
fRMin && (distRMin <= halfCarTolerance))
421 if (distSPhi <= halfAngTolerance)
426 if (distEPhi <= halfAngTolerance)
432 if (distZ <= halfCarTolerance)
448 UWarning, 1,
"Point p is not on surface !?");
449 int oldprc = cout.precision(20);
450 cout <<
"UTubs::SN ( " << p.
x() <<
", " << p.
y() <<
", " << p.
z() <<
" ); "
451 << std::endl << std::endl;
452 cout.precision(oldprc);
456 else if (noSurfaces == 1)
462 norm = sumnorm.
Unit();
480 double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin;
482 rho = std::sqrt(p.
x() * p.
x() + p.
y() * p.
y());
484 distRMin = std::fabs(rho -
fRMin);
485 distRMax = std::fabs(rho -
fRMax);
486 distZ = std::fabs(std::fabs(p.
z()) -
fDz);
488 if (distRMin < distRMax)
490 if (distZ < distRMin)
503 if (distZ < distRMax)
516 phi = std::atan2(p.
y(), p.
x());
529 distSPhi = std::fabs(phi -
fSPhi) * rho;
531 distEPhi = std::fabs(phi -
fSPhi -
fDPhi) * rho;
533 if (distSPhi < distEPhi)
535 if (distSPhi < distMin)
542 if (distEPhi < distMin)
552 norm =
UVector3(-p.
x() / rho, -p.
y() / rho, 0);
587 "Undefined side for valid surface normal to solid.");
619 double tolORMin2, tolIRMax2;
620 double tolORMax2, tolIRMin2, tolODz, tolIDz;
621 const double dRmax = 100.*
fRMax;
628 double Dist, sd, xi, yi, zi, rho2, inum, iden, cosPsi, Comp;
629 double t1, t2, t3, b, c, d;
635 tolORMin2 = (
fRMin - halfRadTolerance) * (
fRMin - halfRadTolerance);
636 tolIRMin2 = (
fRMin + halfRadTolerance) * (
fRMin + halfRadTolerance);
643 tolORMax2 = (
fRMax + halfRadTolerance) * (
fRMax + halfRadTolerance);
644 tolIRMax2 = (
fRMax - halfRadTolerance) * (
fRMax - halfRadTolerance);
648 tolIDz =
fDz - halfCarTolerance;
649 tolODz =
fDz + halfCarTolerance;
651 if (std::fabs(p.
z()) >= tolIDz)
653 if (p.
z() * v.
z() < 0)
655 sd = (std::fabs(p.
z()) -
fDz) / std::fabs(v.
z());
662 xi = p.
x() + sd * v.
x();
663 yi = p.
y() + sd * v.
y();
664 rho2 = xi * xi + yi * yi;
668 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2))
675 iden = std::sqrt(rho2);
676 cosPsi = inum / iden;
690 if (snxt < halfCarTolerance)
710 t1 = 1.0 - v.
z() * v.
z();
711 t2 = p.
x() * v.
x() + p.
y() * v.
y();
712 t3 = p.
x() * p.
x() + p.
y() * p.
y();
718 if ((t3 >= tolORMax2) && (t2 < 0))
728 sd = c / (-b + std::sqrt(d));
734 double fTerm = sd - std::fmod(sd, dRmax);
739 zi = p.
z() + sd * v.
z();
740 if (std::fabs(zi) <= tolODz)
750 xi = p.
x() + sd * v.
x();
751 yi = p.
y() + sd * v.
y();
767 if ((t3 > tolIRMin2) && (t2 < 0) && (std::fabs(p.
z()) <= tolIDz))
774 iden = std::sqrt(t3);
775 cosPsi = inum / iden;
784 c = t3 - fRMax *
fRMax;
795 snxt = c / (-b + std::sqrt(d));
797 if (snxt < halfCarTolerance)
818 c = t3 - fRMax *
fRMax;
829 snxt = c / (-b + std::sqrt(d));
831 if (snxt < halfCarTolerance)
854 sd = (b > 0.) ? c / (-b - std::sqrt(d)) : (-b + std::sqrt(d));
855 if (sd >= -halfCarTolerance)
866 double fTerm = sd - std::fmod(sd, dRmax);
869 zi = p.
z() + sd * v.
z();
870 if (std::fabs(zi) <= tolODz)
880 xi = p.
x() + sd * v.
x();
881 yi = p.
y() + sd * v.
y();
916 if (Dist < halfCarTolerance)
926 zi = p.
z() + sd * v.
z();
927 if (std::fabs(zi) <= tolODz)
929 xi = p.
x() + sd * v.
x();
930 yi = p.
y() + sd * v.
y();
931 rho2 = xi * xi + yi * yi;
933 if (((rho2 >= tolIRMin2) && (rho2 <= tolIRMax2))
934 || ((rho2 > tolORMin2) && (rho2 < tolIRMin2)
937 || ((rho2 > tolIRMax2) && (rho2 < tolORMax2)
962 if (Dist < halfCarTolerance)
972 zi = p.
z() + sd * v.
z();
973 if (std::fabs(zi) <= tolODz)
975 xi = p.
x() + sd * v.
x();
976 yi = p.
y() + sd * v.
y();
977 rho2 = xi * xi + yi * yi;
978 if (((rho2 >= tolIRMin2) && (rho2 <= tolIRMax2))
979 || ((rho2 > tolORMin2) && (rho2 < tolIRMin2)
982 || ((rho2 > tolIRMax2) && (rho2 < tolORMax2)
999 if (snxt < halfCarTolerance)
1034 double safe = 0.0, rho, safe1, safe2, safe3;
1038 rho = std::sqrt(p.
x() * p.
x() + p.
y() * p.
y());
1039 safe1 =
fRMin - rho;
1040 safe2 = rho -
fRMax;
1041 safe3 = std::fabs(p.
z()) -
fDz;
1059 if ((outside) && (safePhi > safe))
1067 safe = 0;
return safe;
1069 if (!aAccurate)
return safe;
1074 safsq += safe1 * safe1;
1079 safsq += safe2 * safe2;
1084 safsq += safe3 * safe3;
1087 if (count == 1)
return safe;
1088 return std::sqrt(safsq);
1101 double deltaR, t1, t2, t3, b, c,
d2, roMin2;
1108 double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2;
1114 pdist =
fDz - p.
z();
1115 if (pdist > halfCarTolerance)
1117 snxt = pdist / v.
z();
1129 pdist =
fDz + p.
z();
1131 if (pdist > halfCarTolerance)
1133 snxt = -pdist / v.
z();
1160 t1 = 1.0 - v.
z() * v.
z();
1161 t2 = p.
x() * v.
x() + p.
y() * v.
y();
1162 t3 = p.
x() * p.
x() + p.
y() * p.
y();
1170 roi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3;
1193 srd = c / (-b - std::sqrt(d2));
1212 roMin2 = t3 - t2 * t2 / t1;
1228 srd = c / (-b + std::sqrt(d2));
1245 srd = -b + std::sqrt(d2);
1266 srd = -b + std::sqrt(d2);
1286 vphi = std::atan2(v.
y(), v.
x());
1288 if (vphi <
fSPhi - halfAngTolerance)
1292 else if (vphi >
fSPhi +
fDPhi + halfAngTolerance)
1313 && (pDistE <= halfCarTolerance)))
1315 && (pDistE > halfCarTolerance))))
1321 sphi = pDistS / compS;
1323 if (sphi >= -halfCarTolerance)
1325 xi = p.
x() + sphi * v.
x();
1326 yi = p.
y() + sphi * v.
y();
1334 if (((
fSPhi - halfAngTolerance) <= vphi)
1335 && ((
fSPhi +
fDPhi + halfAngTolerance) >= vphi))
1347 if (pDistS > -halfCarTolerance)
1365 sphi2 = pDistE / compE;
1369 if ((sphi2 > -halfCarTolerance) && (sphi2 < sphi))
1371 xi = p.
x() + sphi2 * v.
x();
1372 yi = p.
y() + sphi2 * v.
y();
1378 if (!((
fSPhi - halfAngTolerance <= vphi)
1379 && (
fSPhi +
fDPhi + halfAngTolerance >= vphi)))
1382 if (pDistE <= -halfCarTolerance)
1399 if (pDistE <= -halfCarTolerance)
1421 if ((
fSPhi - halfAngTolerance <= vphi)
1422 && (vphi <=
fSPhi +
fDPhi + halfAngTolerance))
1452 xi = p.
x() + snxt * v.
x();
1453 yi = p.
y() + snxt * v.
y();
1459 xi = p.
x() + snxt * v.
x();
1460 yi = p.
y() + snxt * v.
y();
1504 std::ostringstream message;
1505 int oldprc = message.precision(16);
1506 message <<
"Undefined side for valid surface normal to solid."
1508 <<
"Position:" << std::endl << std::endl
1509 <<
"p.x = " << p.
x() <<
" mm" << std::endl
1510 <<
"p.y = " << p.
y() <<
" mm" << std::endl
1511 <<
"p.z = " << p.
z() <<
" mm" << std::endl << std::endl
1512 <<
"Direction:" << std::endl << std::endl
1513 <<
"v.x = " << v.
x() << std::endl
1514 <<
"v.y = " << v.
y() << std::endl
1515 <<
"v.z = " << v.
z() << std::endl << std::endl
1516 <<
"Proposed distance :" << std::endl << std::endl
1517 <<
"snxt = " << snxt <<
" mm" << std::endl;
1518 message.precision(oldprc);
1520 UWarning, 1, message.str().c_str());
1524 if (snxt < halfCarTolerance)
1538 double safe = 0.0, rho, safeZ;
1539 rho = std::sqrt(p.
x() * p.
x() + p.
y() * p.
y());
1544 int oldprc = cout.precision(16);
1547 cout <<
"Position:" << std::endl << std::endl;
1548 cout <<
"p.x = " << p.
x() <<
" mm" << std::endl;
1549 cout <<
"p.y = " << p.
y() <<
" mm" << std::endl;
1550 cout <<
"p.z = " << p.
z() <<
" mm" << std::endl << std::endl;
1551 cout.precision(oldprc);
1553 UWarning, 1,
"Point p is outside !?");
1557 safeZ =
fDz - std::fabs(p.
z());
1578 return std::string(
"Tubs");
1587 return new UTubs(*
this);
1596 int oldprc = os.precision(16);
1597 os <<
"-----------------------------------------------------------\n"
1598 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1599 <<
" ===================================================\n"
1600 <<
" Solid type: UTubs\n"
1601 <<
" Parameters: \n"
1602 <<
" inner radius : " <<
fRMin <<
" mm \n"
1603 <<
" outer radius : " <<
fRMax <<
" mm \n"
1604 <<
" half length Z: " <<
fDz <<
" mm \n"
1607 <<
"-----------------------------------------------------------\n";
1608 os.precision(oldprc);
1619 double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1620 aOne, aTwo, aThr, aFou;
1625 aThr = 0.5 *
fDPhi * (fRMax * fRMax - fRMin *
fRMin);
1629 cosphi = std::cos(phi);
1630 sinphi = std::sin(phi);
1641 if ((chose >= 0) && (chose < aOne))
1643 xRand = fRMax * cosphi;
1644 yRand = fRMax * sinphi;
1646 return UVector3(xRand, yRand, zRand);
1648 else if ((chose >= aOne) && (chose < aOne + aTwo))
1650 xRand = fRMin * cosphi;
1651 yRand = fRMin * sinphi;
1653 return UVector3(xRand, yRand, zRand);
1655 else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr))
1657 xRand = rRand * cosphi;
1658 yRand = rRand * sinphi;
1660 return UVector3(xRand, yRand, zRand);
1662 else if ((chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr))
1664 xRand = rRand * cosphi;
1665 yRand = rRand * sinphi;
1667 return UVector3(xRand, yRand, zRand);
1669 else if ((chose >= aOne + aTwo + 2.*aThr)
1670 && (chose < aOne + aTwo + 2.*aThr + aFou))
1675 return UVector3(xRand, yRand, zRand);
1682 return UVector3(xRand, yRand, zRand);
double GetStartPhiAngle() const
static double frTolerance
double GetZHalfLength() const
double SafetyFromOutside(const UVector3 &p, bool precise=false) const
double SafetyFromInside(const UVector3 &p, bool precise=false) const
virtual void GetParametersList(int, double *) const
const std::string & GetName() const
double GetInnerRadius() const
VUSolid::EnumInside Inside(const UVector3 &p) const
static double Tolerance()
void CheckPhiAngles(double sPhi, double dPhi)
static const double kInfinity
double DistanceToOut(const UVector3 &p, const UVector3 &v, UVector3 &n, bool &validNorm, double aPstep=UUtils::kInfinity) const
double GetDeltaPhiAngle() const
UGeometryType GetEntityType() const
bool Normal(const UVector3 &p, UVector3 &normal) const
virtual UVector3 ApproxSurfaceNormal(const UVector3 &p) const
UTubs & operator=(const UTubs &rhs)
static double faTolerance
UVector3 GetPointOnSurface() const
void Extent(UVector3 &aMin, UVector3 &aMax) const
double GetRadiusInRing(double rmin, double rmax)
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
std::ostream & StreamInfo(std::ostream &os) const
double GetOuterRadius() const
std::string UGeometryType
double SafetyToPhi(const UVector3 &p, const double rho, bool &outside) const
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
double Random(double min=0.0, double max=1.0)
double SafetyFromInsideR(const UVector3 &p, const double rho, bool precise=false) const