File indexing completed on 2025-01-18 10:14:07
0001
0002
0003
0004
0005
0006
0007
0008 #ifndef VOLUMES_HYPEUTILITIES_H_
0009 #define VOLUMES_HYPEUTILITIES_H_
0010
0011 #include "VecGeom/base/Global.h"
0012 #include "VecGeom/volumes/Wedge_Evolution.h"
0013 #include "VecGeom/base/Vector3D.h"
0014 #include "VecGeom/volumes/kernel/GenericKernels.h"
0015 #include "VecGeom/volumes/kernel/shapetypes/HypeTypes.h"
0016 #include <cstdio>
0017
0018 namespace vecgeom {
0019
0020 inline namespace VECGEOM_IMPL_NAMESPACE {
0021
0022 class UnplacedHype;
0023 template <typename T>
0024 struct HypeStruct;
0025
0026 namespace HypeUtilities {
0027 using UnplacedStruct_t = HypeStruct<Precision>;
0028 template <typename Real_v, typename hypeType>
0029 VECCORE_ATT_HOST_DEVICE
0030 typename vecCore::Mask_v<Real_v> IsCompletelyOutside(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point)
0031 {
0032 using namespace ::vecgeom::HypeTypes;
0033 using Bool_v = typename vecCore::Mask_v<Real_v>;
0034 Real_v r2 = point.Perp2();
0035 Real_v oRad2 = (hype.fRmax2 + hype.fTOut2 * point.z() * point.z());
0036
0037 Bool_v completelyoutside = (Abs(point.z()) > (hype.fDz + hype.zToleranceLevel));
0038 if (vecCore::MaskFull(completelyoutside)) return completelyoutside;
0039 completelyoutside |= (r2 > oRad2 + hype.outerRadToleranceLevel);
0040 if (vecCore::MaskFull(completelyoutside)) return completelyoutside;
0041
0042
0043 if (checkInnerSurfaceTreatment<hypeType>(hype)) {
0044 Real_v iRad2 = (hype.fRmin2 + hype.fTIn2 * point.z() * point.z());
0045 completelyoutside |= (r2 < (iRad2 - hype.innerRadToleranceLevel));
0046 }
0047 return completelyoutside;
0048 }
0049
0050 template <typename Real_v, typename hypeType>
0051 VECCORE_ATT_HOST_DEVICE
0052 typename vecCore::Mask_v<Real_v> IsCompletelyInside(UnplacedStruct_t const &hype, Vector3D<Real_v> const &point)
0053 {
0054 using namespace ::vecgeom::HypeTypes;
0055 using Bool_v = typename vecCore::Mask_v<Real_v>;
0056 Real_v r2 = point.Perp2();
0057 Real_v oRad2 = (hype.fRmax2 + hype.fTOut2 * point.z() * point.z());
0058
0059 Bool_v completelyinside =
0060 (Abs(point.z()) < (hype.fDz - hype.zToleranceLevel)) && (r2 < oRad2 - hype.outerRadToleranceLevel);
0061
0062 if (checkInnerSurfaceTreatment<hypeType>(hype)) {
0063 Real_v iRad2 = (hype.fRmin2 + hype.fTIn2 * point.z() * point.z());
0064 completelyinside &= (r2 > (iRad2 + hype.innerRadToleranceLevel));
0065 }
0066 return completelyinside;
0067 }
0068
0069 template <typename Real_v, bool ForInnerRad>
0070 VECCORE_ATT_HOST_DEVICE
0071 Real_v RadiusHypeSq(UnplacedStruct_t const &hype, Real_v z)
0072 {
0073
0074 if (ForInnerRad)
0075 return (hype.fRmin2 + hype.fTIn2 * z * z);
0076 else
0077 return (hype.fRmax2 + hype.fTOut2 * z * z);
0078 }
0079
0080 template <typename Real_v>
0081 VECCORE_ATT_HOST_DEVICE
0082 typename vecCore::Mask_v<Real_v> IsPointMovingInsideOuterSurface(UnplacedStruct_t const &hype,
0083 Vector3D<Real_v> const &point,
0084 Vector3D<Real_v> const &direction)
0085 {
0086 Real_v pz = point.z();
0087 Real_v vz = direction.z();
0088 vecCore__MaskedAssignFunc(vz, pz < Real_v(0.), -vz);
0089 vecCore__MaskedAssignFunc(pz, pz < Real_v(0.), -pz);
0090 return ((point.x() * direction.x() + point.y() * direction.y() - pz * hype.fTOut2 * vz) < Real_v(0.));
0091 }
0092
0093 template <typename Real_v>
0094 VECCORE_ATT_HOST_DEVICE
0095 typename vecCore::Mask_v<Real_v> IsPointMovingInsideInnerSurface(UnplacedStruct_t const &hype,
0096 Vector3D<Real_v> const &point,
0097 Vector3D<Real_v> const &direction)
0098 {
0099 Real_v pz = point.z();
0100 Real_v vz = direction.z();
0101
0102 vecCore__MaskedAssignFunc(vz, pz < Real_v(0.), -vz);
0103 vecCore__MaskedAssignFunc(pz, pz < Real_v(0.), -pz);
0104
0105
0106 return ((point.x() * direction.x() + point.y() * direction.y() - pz * hype.fTIn2 * vz) > Real_v(0.));
0107 }
0108
0109 template <typename Real_v, typename hypeType>
0110 VECCORE_ATT_HOST_DEVICE
0111 typename vecCore::Mask_v<Real_v> IsPointOnSurfaceAndMovingInside(UnplacedStruct_t const &hype,
0112 Vector3D<Real_v> const &point,
0113 Vector3D<Real_v> const &direction)
0114 {
0115 using namespace ::vecgeom::HypeTypes;
0116 using Bool_v = typename vecCore::Mask_v<Real_v>;
0117 Bool_v innerHypeSurf(false), outerHypeSurf(false), zSurf(false);
0118 Bool_v done(false);
0119 Real_v rho2 = point.Perp2();
0120 Real_v radI2 = RadiusHypeSq<Real_v, true>(hype, point.z());
0121 Real_v radO2 = RadiusHypeSq<Real_v, false>(hype, point.z());
0122
0123 Bool_v in(false);
0124 zSurf = ((rho2 - hype.fEndOuterRadius2) < kTolerance) && ((hype.fEndInnerRadius2 - rho2) < kTolerance) &&
0125 (Abs(Abs(point.z()) - hype.fDz) < kTolerance);
0126 in |= (zSurf && (point.z() * direction.z() < Real_v(0.)));
0127
0128 done |= zSurf;
0129 if (vecCore::MaskFull(done)) return in;
0130
0131 outerHypeSurf |= (!zSurf && (Abs((radO2) - (rho2)) < hype.outerRadToleranceLevel));
0132 in |= (!done && outerHypeSurf && IsPointMovingInsideOuterSurface<Real_v>(hype, point, direction));
0133
0134
0135 if (checkInnerSurfaceTreatment<hypeType>(hype)) {
0136 done |= (!zSurf && outerHypeSurf);
0137 if (vecCore::MaskFull(done)) return in;
0138
0139 innerHypeSurf |= (!zSurf && !outerHypeSurf && (Abs((radI2) - (rho2)) < hype.innerRadToleranceLevel));
0140 in |= (!done && !zSurf && innerHypeSurf && IsPointMovingInsideInnerSurface<Real_v>(hype, point, direction));
0141 done |= (!zSurf && innerHypeSurf);
0142 if (vecCore::MaskFull(done)) return in;
0143 }
0144 return in;
0145 }
0146
0147 template <typename Real_v, typename hypeType, bool ForDistToIn>
0148 VECCORE_ATT_HOST_DEVICE
0149 typename vecCore::Mask_v<Real_v> GetPointOfIntersectionWithZPlane(UnplacedStruct_t const &hype,
0150 Vector3D<Real_v> const &point,
0151 Vector3D<Real_v> const &direction, Real_v &zDist)
0152 {
0153 using namespace ::vecgeom::HypeTypes;
0154 zDist = (Sign(ForDistToIn ? point.z() : direction.z()) * hype.fDz - point.z()) / NonZero(direction.z());
0155
0156 auto r2 = (point + zDist * direction).Perp2();
0157
0158 if (!checkInnerSurfaceTreatment<hypeType>(hype))
0159 return (r2 < hype.fEndOuterRadius2);
0160 else
0161 return ((r2 < hype.fEndOuterRadius2) && (r2 > hype.fEndInnerRadius2));
0162 }
0163
0164 template <typename Real_v>
0165 VECCORE_ATT_HOST_DEVICE
0166 typename vecCore::Mask_v<Real_v> IsPointMovingOutsideOuterSurface(UnplacedStruct_t const &hype,
0167 Vector3D<Real_v> const &point,
0168 Vector3D<Real_v> const &direction)
0169 {
0170
0171 using Bool_v = typename vecCore::Mask_v<Real_v>;
0172 Bool_v out(false);
0173
0174 Real_v pz = point.z();
0175 Real_v vz = direction.z();
0176 vecCore__MaskedAssignFunc(pz, vz < Real_v(0.), -pz);
0177 vecCore__MaskedAssignFunc(vz, vz < Real_v(0.), -vz);
0178 Vector3D<Real_v> normHere(point.x(), point.y(), -point.z() * hype.fTOut2);
0179 out = (normHere.Dot(direction) > Real_v(0.));
0180 return out;
0181 }
0182
0183 template <typename Real_v>
0184 VECCORE_ATT_HOST_DEVICE
0185 typename vecCore::Mask_v<Real_v> IsPointMovingOutsideInnerSurface(UnplacedStruct_t const &hype,
0186 Vector3D<Real_v> const &point,
0187 Vector3D<Real_v> const &direction)
0188 {
0189
0190 Real_v pz = point.z();
0191 Real_v vz = direction.z();
0192 vecCore__MaskedAssignFunc(pz, vz < Real_v(0.), -pz);
0193 vecCore__MaskedAssignFunc(vz, vz < Real_v(0.), -vz);
0194 Vector3D<Real_v> normHere(-point.x(), -point.y(), point.z() * hype.fTIn2);
0195 return (normHere.Dot(direction) > Real_v(0.));
0196 }
0197
0198 template <typename Real_v>
0199 VECCORE_ATT_HOST_DEVICE
0200 typename vecCore::Mask_v<Real_v> IsPointOnOuterSurfaceAndMovingOutside(UnplacedStruct_t const &hype,
0201 Vector3D<Real_v> const &point,
0202 Vector3D<Real_v> const &direction)
0203 {
0204
0205 using Bool_v = typename vecCore::Mask_v<Real_v>;
0206 Real_v rho2 = point.x() * point.x() + point.y() * point.y();
0207 Real_v absZ = Abs(point.z());
0208 Real_v radO2 = RadiusHypeSq<Real_v, false>(hype, point.z());
0209 Bool_v out(false), outerHypeSurf(false);
0210 outerHypeSurf = (Abs((radO2) - (rho2)) < hype.outerRadToleranceLevel) && (absZ >= Real_v(0.)) && (absZ < hype.fDz);
0211 out = outerHypeSurf && IsPointMovingOutsideOuterSurface<Real_v>(hype, point, direction);
0212 return out;
0213 }
0214
0215 template <typename Real_v, typename hypeType>
0216 VECCORE_ATT_HOST_DEVICE
0217 typename vecCore::Mask_v<Real_v> IsPointOnInnerSurfaceAndMovingOutside(UnplacedStruct_t const &hype,
0218 Vector3D<Real_v> const &point,
0219 Vector3D<Real_v> const &direction)
0220 {
0221 using namespace ::vecgeom::HypeTypes;
0222 using Bool_v = typename vecCore::Mask_v<Real_v>;
0223 Real_v rho2 = point.x() * point.x() + point.y() * point.y();
0224 Real_v absZ = Abs(point.z());
0225 Real_v radI2 = RadiusHypeSq<Real_v, true>(hype, point.z());
0226 Bool_v out(false), innerHypeSurf(false);
0227
0228 if (checkInnerSurfaceTreatment<hypeType>(hype)) {
0229 innerHypeSurf = (Abs((radI2) - (rho2)) < hype.innerRadToleranceLevel) && (absZ >= Real_v(0.)) && (absZ < hype.fDz);
0230 out = innerHypeSurf && HypeUtilities::IsPointMovingOutsideInnerSurface<Real_v>(hype, point, direction);
0231 }
0232 return out;
0233 }
0234
0235 template <typename Real_v, typename hypeType>
0236 VECCORE_ATT_HOST_DEVICE
0237 typename vecCore::Mask_v<Real_v> IsPointOnSurfaceAndMovingOutside(UnplacedStruct_t const &hype,
0238 Vector3D<Real_v> const &point,
0239 Vector3D<Real_v> const &direction)
0240 {
0241 using namespace ::vecgeom::HypeTypes;
0242 using Bool_v = typename vecCore::Mask_v<Real_v>;
0243 Bool_v innerHypeSurf(false), outerHypeSurf(false), zSurf(false);
0244 Bool_v done(false);
0245
0246 Real_v rho2 = point.x() * point.x() + point.y() * point.y();
0247 Real_v radI2 = RadiusHypeSq<Real_v, true>(hype, point.z());
0248 Real_v radO2 = RadiusHypeSq<Real_v, false>(hype, point.z());
0249
0250 Bool_v out(false);
0251
0252
0253 zSurf = ((rho2 - hype.fEndOuterRadius2) < kTolerance) && ((hype.fEndInnerRadius2 - rho2) < kTolerance) &&
0254 (Abs(Abs(point.z()) - hype.fDz) < kTolerance);
0255
0256 out |= (zSurf && (point.z() * direction.z() > Real_v(0.)));
0257
0258 done = out;
0259 if (vecCore::MaskFull(done)) return out;
0260
0261 outerHypeSurf |= !done && (Abs(radO2 - rho2) < hype.outerRadToleranceLevel);
0262
0263 out |= (outerHypeSurf && HypeUtilities::IsPointMovingOutsideOuterSurface<Real_v>(hype, point, direction));
0264
0265
0266 done |= out;
0267 if (vecCore::MaskFull(done)) return out;
0268
0269
0270 if (checkInnerSurfaceTreatment<hypeType>(hype)) {
0271
0272 innerHypeSurf |= (!done && (Abs(radI2 - rho2) < hype.innerRadToleranceLevel));
0273
0274 out |= (innerHypeSurf && HypeUtilities::IsPointMovingOutsideInnerSurface<Real_v>(hype, point, direction));
0275
0276 done |= out;
0277 if (vecCore::MaskFull(done)) return out;
0278 }
0279
0280 return out;
0281 }
0282
0283 template <typename Real_v>
0284 VECCORE_ATT_HOST_DEVICE
0285 Real_v ApproxDistOutside(Real_v pr, Real_v pz, Precision r0, Precision tanPhi)
0286 {
0287 Real_v r1 = Sqrt(r0 * r0 + tanPhi * tanPhi * pz * pz);
0288 Real_v z1 = pz;
0289 Real_v r2 = pr;
0290 Real_v z2 = Sqrt((pr * pr - r0 * r0) / (tanPhi * tanPhi));
0291 Real_v dz = z2 - z1;
0292 Real_v dr = r2 - r1;
0293 Real_v r3 = Sqrt(dr * dr + dz * dz);
0294 auto mask = r3 < vecCore::NumericLimits<Real_v>::Min();
0295 return vecCore::Blend(mask, (r2 - r1), (r2 - r1) * dz / r3);
0296 }
0297
0298 template <typename Real_v>
0299 VECCORE_ATT_HOST_DEVICE
0300 Real_v ApproxDistInside(Real_v pr, Real_v pz, Precision r0, Precision tan2Phi)
0301 {
0302 using Bool_v = typename vecCore::Mask_v<Real_v>;
0303 Bool_v done(false);
0304 Real_v ret(0.);
0305 Real_v tan2Phi_v(tan2Phi);
0306 vecCore__MaskedAssignFunc(ret, (tan2Phi_v < vecCore::NumericLimits<Real_v>::Min()), r0 - pr);
0307 done |= (tan2Phi_v < vecCore::NumericLimits<Real_v>::Min());
0308 if (vecCore::MaskFull(done)) return ret;
0309
0310 Real_v rh = Sqrt(r0 * r0 + pz * pz * tan2Phi_v);
0311 Real_v dr = -rh;
0312 Real_v dz = pz * tan2Phi_v;
0313 Real_v len = Sqrt(dr * dr + dz * dz);
0314
0315 vecCore__MaskedAssignFunc(ret, !done, Abs((pr - rh) * dr) / len);
0316 return ret;
0317 }
0318
0319 }
0320
0321
0322
0323
0324 template <class Real_v, bool ForDistToIn, bool ForInnerSurface>
0325 class HypeHelpers {
0326
0327 public:
0328 HypeHelpers() {}
0329 ~HypeHelpers() {}
0330
0331 VECCORE_ATT_HOST_DEVICE
0332 static typename vecCore::Mask_v<Real_v> GetPointOfIntersectionWithHyperbolicSurface(HypeStruct<Precision> const &hype,
0333 Vector3D<Real_v> const &point,
0334 Vector3D<Real_v> const &direction,
0335 Real_v &dist)
0336 {
0337
0338 using Bool_v = typename vecCore::Mask_v<Real_v>;
0339
0340 if (ForInnerSurface) {
0341 Real_v a = direction.Perp2() - hype.fTIn2 * direction.z() * direction.z();
0342 Real_v b = (direction.x() * point.x() + direction.y() * point.y() - hype.fTIn2 * direction.z() * point.z());
0343 Real_v c = point.Perp2() - hype.fTIn2 * point.z() * point.z() - hype.fRmin2;
0344 Bool_v exist = (b * b - a * c > Real_v(0.));
0345 if (ForDistToIn) {
0346 vecCore__MaskedAssignFunc(dist, exist && b < Real_v(0.), ((-b + Sqrt(b * b - a * c)) / (a)));
0347 vecCore__MaskedAssignFunc(dist, exist && b >= Real_v(0.), ((c) / (-b - Sqrt(b * b - a * c))));
0348
0349 } else {
0350 vecCore__MaskedAssignFunc(dist, exist && b > Real_v(0.), ((-b - Sqrt(b * b - a * c)) / (a)));
0351 vecCore__MaskedAssignFunc(dist, exist && b <= Real_v(0.), ((c) / (-b + Sqrt(b * b - a * c))));
0352 }
0353
0354 } else {
0355 Real_v a = direction.Perp2() - hype.fTOut2 * direction.z() * direction.z();
0356 Real_v b = (direction.x() * point.x() + direction.y() * point.y() - hype.fTOut2 * direction.z() * point.z());
0357 Real_v c = point.Perp2() - hype.fTOut2 * point.z() * point.z() - hype.fRmax2;
0358 Bool_v exist = (b * b - a * c > Real_v(0.));
0359 if (ForDistToIn) {
0360 vecCore__MaskedAssignFunc(dist, exist && b >= Real_v(0.), ((-b - Sqrt(b * b - a * c)) / (a)));
0361 vecCore__MaskedAssignFunc(dist, exist && b < Real_v(0.), ((c) / (-b + Sqrt(b * b - a * c))));
0362 } else {
0363 vecCore__MaskedAssignFunc(dist, exist && b < Real_v(0.), ((-b + Sqrt(b * b - a * c)) / (a)));
0364 vecCore__MaskedAssignFunc(dist, exist && b >= Real_v(0.), ((c) / (-b - Sqrt(b * b - a * c))));
0365 }
0366 }
0367
0368 vecCore__MaskedAssignFunc(dist, dist < Real_v(0.), InfinityLength<Real_v>());
0369
0370 Real_v newPtZ = point.z() + dist * direction.z();
0371
0372 return (Abs(newPtZ) <= hype.fDz);
0373 }
0374 };
0375
0376 template <bool ForDistToIn, bool ForInnerSurface>
0377 class HypeHelpers<Precision, ForDistToIn, ForInnerSurface> {
0378 public:
0379 HypeHelpers() {}
0380 ~HypeHelpers() {}
0381
0382 VECCORE_ATT_HOST_DEVICE
0383 static bool GetPointOfIntersectionWithHyperbolicSurface(HypeStruct<Precision> const &hype,
0384 Vector3D<Precision> const &point,
0385 Vector3D<Precision> const &direction, Precision &dist)
0386 {
0387 if (ForInnerSurface) {
0388 Precision a = direction.Perp2() - hype.fTIn2 * direction.z() * direction.z();
0389 Precision b = (direction.x() * point.x() + direction.y() * point.y() - hype.fTIn2 * direction.z() * point.z());
0390 Precision c = point.Perp2() - hype.fTIn2 * point.z() * point.z() - hype.fRmin2;
0391 bool exist = (b * b - a * c > 0.);
0392 if (exist) {
0393 if (ForDistToIn) {
0394 if (b < 0.)
0395 dist = ((-b + Sqrt(b * b - a * c)) / (a));
0396 else
0397 dist = ((c) / (-b - Sqrt(b * b - a * c)));
0398 } else {
0399
0400 if (b > 0.)
0401 dist = ((-b - Sqrt(b * b - a * c)) / (a));
0402 else
0403 dist = ((c) / (-b + Sqrt(b * b - a * c)));
0404 }
0405 } else
0406 return false;
0407
0408 } else {
0409 Precision a = direction.Perp2() - hype.fTOut2 * direction.z() * direction.z();
0410 Precision b = (direction.x() * point.x() + direction.y() * point.y() - hype.fTOut2 * direction.z() * point.z());
0411 Precision c = point.Perp2() - hype.fTOut2 * point.z() * point.z() - hype.fRmax2;
0412 bool exist = (b * b - a * c > 0.);
0413
0414 if (exist) {
0415 if (ForDistToIn) {
0416 if (b >= 0.)
0417 dist = ((-b - Sqrt(b * b - a * c)) / (a));
0418 else
0419 dist = ((c) / (-b + Sqrt(b * b - a * c)));
0420 } else {
0421
0422 if (b < 0.)
0423 dist = ((-b + Sqrt(b * b - a * c)) / a);
0424 else
0425 dist = (c / (-b - Sqrt(b * b - a * c)));
0426 }
0427 } else
0428 return false;
0429 }
0430 if (dist < 0.) dist = kInfLength;
0431
0432
0433 Precision newPtZ = point.z() + dist * direction.z();
0434
0435 return (Abs(newPtZ) <= hype.fDz);
0436 }
0437 };
0438
0439 }
0440 }
0441
0442 #endif
