File indexing completed on 2026-09-26 09:15:49
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #ifndef VECGEOM_VOLUMES_KERNEL_POLYCONEIMPLEMENTATION_H_
0013 #define VECGEOM_VOLUMES_KERNEL_POLYCONEIMPLEMENTATION_H_
0014
0015 #include "VecGeom/base/Vector3D.h"
0016 #include "VecGeom/volumes/PolyconeStruct.h"
0017 #include "VecGeom/volumes/kernel/GenericKernels.h"
0018 #include <VecCore/VecCore>
0019 #include "VecGeom/volumes/kernel/ConeImplementation.h"
0020 #include "VecGeom/volumes/kernel/shapetypes/ConeTypes.h"
0021
0022 #include <cstdio>
0023
0024 namespace vecgeom {
0025
0026 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(struct, PolyconeImplementation, typename);
0027
0028 inline namespace VECGEOM_IMPL_NAMESPACE {
0029
0030 template <typename T>
0031 class SPlacedPolycone;
0032 template <typename T>
0033 class SUnplacedPolycone;
0034
0035 template <typename polyconeTypeT>
0036 struct PolyconeImplementation {
0037
0038 using UnplacedStruct_t = PolyconeStruct<Precision>;
0039 using UnplacedVolume_t = SUnplacedPolycone<polyconeTypeT>;
0040 using PlacedShape_t = SPlacedPolycone<UnplacedVolume_t>;
0041
0042 template <typename Real_v, bool ForInside>
0043 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForASection(
0044 UnplacedStruct_t const &unplaced, int isect, Vector3D<Real_v> const &polyconePoint,
0045 typename vecCore::Mask_v<Real_v> &secFullyInside, typename vecCore::Mask_v<Real_v> &secFullyOutside)
0046 {
0047
0048
0049 using namespace ConeTypes;
0050
0051 if (isect < 0) {
0052 secFullyInside = false;
0053 secFullyOutside = true;
0054 return;
0055 }
0056
0057 PolyconeSection const &sec = unplaced.GetSection(isect);
0058 Vector3D<Precision> secLocalp = polyconePoint - Vector3D<Precision>(0, 0, sec.fShift);
0059 #ifdef POLYCONEDEBUG
0060 std::cerr << " isect=" << isect << "/" << unplaced.GetNSections() << " secLocalP=" << secLocalp
0061 << ", secShift=" << sec.fShift << " sec.fSolid=" << sec.fSolid << std::endl;
0062 if (sec.fSolid) sec.fSolid->Print();
0063 #endif
0064
0065 ConeHelpers<Real_v, polyconeTypeT>::template GenericKernelForContainsAndInside<ForInside>(
0066 sec.fSolid, secLocalp, secFullyInside, secFullyOutside);
0067 }
0068
0069 template <typename Real_v, typename Bool_v>
0070 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &polycone,
0071 Vector3D<Real_v> const &point, Bool_v &inside)
0072 {
0073
0074 Bool_v unused(false), outside(false);
0075 GenericKernelForContainsAndInside<Real_v, Bool_v, false>(polycone, point, unused, outside);
0076 inside = !outside;
0077 }
0078
0079
0080
0081 template <typename Real_v, typename Inside_t>
0082 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &polycone,
0083 Vector3D<Real_v> const &point, Inside_t &inside)
0084 {
0085
0086 using Bool_v = vecCore::Mask_v<Real_v>;
0087 using InsideBool_v = vecCore::Mask_v<Inside_t>;
0088 Bool_v completelyinside(false), completelyoutside(false);
0089 GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, completelyinside, completelyoutside);
0090 inside = EInside::kSurface;
0091 vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0092 vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0093 }
0094
0095 template <typename Real_v, typename Bool_v, bool ForInside>
0096 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0097 UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &localPoint, Bool_v &completelyInside,
0098 Bool_v &completelyOutside)
0099 {
0100
0101 typedef Bool_v Bool_t;
0102
0103 int indexLow = unplaced.GetSectionIndex(localPoint.z() - kTolerance);
0104 int indexHigh = unplaced.GetSectionIndex(localPoint.z() + kTolerance);
0105 if (indexLow < 0 && indexHigh < 0) {
0106 completelyOutside = true;
0107 return;
0108 }
0109 if (indexLow < 0 && indexHigh == 0) {
0110
0111 GenericKernelForASection<Real_v, ForInside>(unplaced, 0, localPoint, completelyInside, completelyOutside);
0112 return;
0113 }
0114 if (indexHigh < 0 && indexLow == (unplaced.GetNSections() - 1)) {
0115
0116 GenericKernelForASection<Real_v, ForInside>(unplaced, (unplaced.GetNSections() - 1), localPoint, completelyInside,
0117 completelyOutside);
0118
0119 return;
0120 }
0121 if (indexLow >= 0 && indexHigh >= 0) {
0122 if (indexLow == indexHigh) {
0123
0124 GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, completelyInside,
0125 completelyOutside);
0126
0127 return;
0128 } else {
0129
0130 Bool_t secInLow = false, secOutLow = false;
0131 Bool_t secInHigh = false, secOutHigh = false;
0132 GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, secInLow, secOutLow);
0133 GenericKernelForASection<Real_v, ForInside>(unplaced, indexHigh, localPoint, secInHigh, secOutHigh);
0134 Bool_t surfLow = !secInLow && !secOutLow;
0135 Bool_t surfHigh = !secInHigh && !secOutHigh;
0136
0137 if (surfLow && surfHigh) {
0138 completelyInside = true;
0139 return;
0140 } else {
0141
0142
0143
0144
0145
0146
0147 if (secOutLow && secOutHigh) {
0148 completelyOutside = true;
0149 return;
0150 }
0151 }
0152 }
0153 }
0154 }
0155
0156 template <typename Real_v>
0157 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &polycone,
0158 Vector3D<Real_v> const &point,
0159 Vector3D<Real_v> const &direction,
0160 Real_v const &stepMax, Real_v &distance)
0161 {
0162
0163 Vector3D<Real_v> p = point;
0164 Vector3D<Real_v> v = direction;
0165
0166 #ifdef POLYCONEDEBUG
0167 std::cerr << "Polycone::DistToIn() (spot 1): point=" << point << ", dir=" << direction << ", localPoint=" << p
0168 << ", localDir=" << v << "\n";
0169 #endif
0170
0171
0172
0173 distance = kInfLength;
0174 int increment = (v.z() > 0) ? 1 : -1;
0175 if (std::fabs(v.z()) < kTolerance) increment = 0;
0176 int index = polycone.GetSectionIndex(p.z());
0177 if (index == -1) index = 0;
0178 if (index == -2) index = polycone.GetNSections() - 1;
0179
0180 do {
0181
0182 PolyconeSection const &sec = polycone.GetSection(index);
0183
0184 #ifdef POLYCONEDEBUG
0185 std::cerr << "Polycone::DistToIn() (spot 2):"
0186 << " index=" << index << " NSec=" << polycone.GetNSections() << " &sec=" << &sec << " - secPars:"
0187 << " secOffset=" << sec.fShift << " Dz=" << sec.fSolid.GetDz() << " Rmin1=" << sec.fSolid.GetRmin1()
0188 << " Rmin2=" << sec.fSolid.GetRmin2() << " Rmax1=" << sec.fSolid.GetRmax1()
0189 << " Rmax2=" << sec.fSolid.GetRmax2() << " -- calling Cone::DistToIn()...\n";
0190 #endif
0191
0192 ConeImplementation<polyconeTypeT>::template DistanceToIn<Real_v>(
0193 sec.fSolid, p - Vector3D<Precision>(0, 0, sec.fShift), v, stepMax, distance);
0194
0195 #ifdef POLYCONEDEBUG
0196 std::cerr << "Polycone::DistToIn() (spot 3):"
0197 << " distToIn() = " << distance << "\n";
0198 #endif
0199
0200 if (distance < kInfLength || !increment) break;
0201 index += increment;
0202 } while (index >= 0 && index < polycone.GetNSections());
0203 return;
0204 }
0205
0206 template <typename Real_v>
0207 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &polycone,
0208 Vector3D<Real_v> const &point,
0209 Vector3D<Real_v> const &dir,
0210 Real_v const &stepMax, Real_v &distance)
0211 {
0212
0213 distance = kInfLength;
0214 Vector3D<Real_v> pn = point;
0215
0216
0217 if (polycone.GetNSections() == 1) {
0218 const PolyconeSection §ion = polycone.GetSection(0);
0219
0220 ConeImplementation<polyconeTypeT>::template DistanceToOut<Real_v>(
0221 section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), dir, stepMax, distance);
0222
0223 return;
0224 }
0225
0226 int indexLow = polycone.GetSectionIndex(point.z() - kTolerance);
0227 int indexHigh = polycone.GetSectionIndex(point.z() + kTolerance);
0228 int index = 0;
0229
0230
0231
0232
0233 if (indexLow < 0 && indexHigh < 0) {
0234 distance = -1;
0235 return;
0236 } else if (indexLow < 0 && indexHigh >= 0) {
0237 index = indexHigh;
0238 const PolyconeSection §ion = polycone.GetSection(index);
0239
0240 Inside_t inside;
0241
0242
0243 ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0244 section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0245 if (inside == EInside::kOutside) {
0246 distance = -1;
0247 return;
0248 }
0249 } else if (indexLow != indexHigh && (indexLow >= 0)) {
0250
0251 const PolyconeSection §ion = polycone.GetSection(indexLow);
0252
0253 Inside_t inside;
0254
0255
0256
0257 ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0258 section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0259
0260 if (inside == EInside::kOutside) {
0261 index = indexHigh;
0262 } else {
0263 index = indexLow;
0264 }
0265 } else {
0266 index = indexLow;
0267 if (index < 0) index = polycone.GetSectionIndex(point.z());
0268 }
0269 if (index < 0) {
0270 distance = 0.;
0271 return;
0272 }
0273
0274 else {
0275 const PolyconeSection §ion = polycone.GetSection(index);
0276
0277 Inside_t inside;
0278
0279
0280 ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0281 section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0282 if (inside == EInside::kOutside) {
0283 distance = -1;
0284 return;
0285 }
0286 }
0287
0288 Precision totalDistance = 0.;
0289 Precision dist;
0290 int increment = (dir.z() > 0) ? 1 : -1;
0291 if (std::fabs(dir.z()) < kTolerance) increment = 0;
0292
0293
0294 int istep = 0;
0295 do {
0296 const PolyconeSection §ion = polycone.GetSection(index);
0297
0298 if ((totalDistance != 0) || (istep < 2)) {
0299 pn = point + totalDistance * dir;
0300 pn.z() -= section.fShift;
0301 Inside_t inside;
0302
0303 ConeImplementation<polyconeTypeT>::template Inside<Real_v>(section.fSolid, pn, inside);
0304
0305 if (inside == EInside::kOutside) {
0306 break;
0307 }
0308 } else
0309 pn.z() -= section.fShift;
0310
0311 istep++;
0312
0313
0314 ConeImplementation<polyconeTypeT>::template DistanceToOut<Real_v>(section.fSolid, pn, dir, stepMax, dist);
0315 if (dist == -1) return;
0316
0317
0318 if (std::fabs(dist) < 0.5 * kTolerance) {
0319 int index1 = index;
0320 if ((index > 0) && (index < polycone.GetNSections() - 1)) {
0321 index1 += increment;
0322 } else {
0323 if ((index == 0) && (increment > 0)) index1 += increment;
0324 if ((index == polycone.GetNSections() - 1) && (increment < 0)) index1 += increment;
0325 }
0326
0327 Vector3D<Precision> pte = point + (totalDistance + dist) * dir;
0328 const PolyconeSection §ion1 = polycone.GetSection(index1);
0329 pte.z() -= section1.fShift;
0330 Vector3D<Precision> localp;
0331 Inside_t inside22;
0332
0333 ConeImplementation<polyconeTypeT>::template Inside<Real_v>(section1.fSolid, pte, inside22);
0334 if (inside22 == 3 || (increment == 0)) {
0335 break;
0336 }
0337 }
0338
0339 totalDistance += dist;
0340 index += increment;
0341 } while (increment != 0 && index >= 0 && index < polycone.GetNSections());
0342
0343 distance = totalDistance;
0344
0345 return;
0346 }
0347
0348 template <typename Real_v>
0349 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &polycone,
0350 Vector3D<Real_v> const &point, Real_v &safety)
0351 {
0352
0353 Vector3D<Real_v> p = point;
0354 int index = polycone.GetSectionIndex(p.z());
0355
0356 bool needZ = false;
0357 if (index < 0) {
0358 needZ = true;
0359 if (index == -1) index = 0;
0360 if (index == -2) index = polycone.GetNSections() - 1;
0361 }
0362 Precision minSafety = 0;
0363 PolyconeSection const &sec = polycone.GetSection(index);
0364
0365 if (needZ) {
0366
0367
0368 ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sec.fSolid,
0369 p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0370 } else {
0371
0372
0373
0374 ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sec.fSolid,
0375 p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0376
0377 if (safety < kTolerance) return;
0378 minSafety = safety;
0379 Precision zbase = polycone.fZs[index + 1];
0380
0381 for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0382 Precision dz = polycone.fZs[i] - zbase;
0383 if (dz >= minSafety) break;
0384
0385 PolyconeSection const § = polycone.GetSection(i);
0386
0387
0388
0389
0390 ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(
0391 sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0392
0393 if (safety < minSafety) minSafety = safety;
0394 }
0395
0396
0397 if (index > 0) {
0398 zbase = polycone.fZs[index - 1];
0399 for (int i = index - 1; i >= 0; --i) {
0400 Precision dz = zbase - polycone.fZs[i];
0401 if (dz >= minSafety) break;
0402 PolyconeSection const § = polycone.GetSection(i);
0403
0404
0405
0406
0407 ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(
0408 sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0409
0410 if (safety < minSafety) minSafety = safety;
0411 }
0412 }
0413 safety = minSafety;
0414 }
0415 return;
0416 }
0417
0418 template <typename Real_v>
0419 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Real_v DistanceToSeg(Vector3D<Real_v> const &point,
0420 Vector3D<Precision> segment_start,
0421 Vector3D<Precision> segment_end)
0422 {
0423 Vector3D<Real_v> segment_vector = segment_end - segment_start;
0424 Vector3D<Real_v> point_vector = point - segment_start;
0425 Real_v projection_scalar = point_vector.Dot(segment_vector) / segment_vector.Mag2();
0426 Vector3D<Real_v> projection_point = segment_start + projection_scalar * segment_vector;
0427 vecCore__MaskedAssignFunc(projection_point, projection_scalar < Real_v(0.), Vector3D<Real_v>(segment_start));
0428 vecCore__MaskedAssignFunc(projection_point, projection_scalar > Real_v(1.), Vector3D<Real_v>(segment_end));
0429 Vector3D<Real_v> distVec = point - projection_point;
0430 return distVec.Mag();
0431 }
0432
0433
0434
0435
0436 template <typename Real_v>
0437 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &polycone,
0438 Vector3D<Real_v> const &point, Real_v &safety)
0439 {
0440 typedef typename vecCore::Mask_v<Real_v> Bool_v;
0441 safety = Real_v(kInfLength);
0442 Bool_v compIn(false), compOut(false), done(false);
0443 GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, compIn, compOut);
0444 vecCore__MaskedAssignFunc(safety, compOut && !done, Real_v(-1.));
0445 done |= compOut;
0446 if (vecCore::MaskFull(done)) return;
0447
0448 Vector3D<Real_v> twoDPoint = Vector3D<Real_v>(point.Perp(), point.z(), 0.);
0449 for (unsigned int currSegIndex = 0; currSegIndex < polycone.fTwoDVec.size(); currSegIndex++) {
0450 unsigned int nextSegIndex = currSegIndex + 1;
0451 if (currSegIndex == polycone.fTwoDVec.size() - 1) {
0452 nextSegIndex = 0;
0453 }
0454
0455 Real_v distance(kInfLength);
0456 if ((polycone.fTwoDVec[currSegIndex].x() == 0 && polycone.fTwoDVec[nextSegIndex].x() == 0) ||
0457 (polycone.fTwoDVec[currSegIndex].x() == polycone.fTwoDVec[nextSegIndex].x() &&
0458 polycone.fTwoDVec[currSegIndex].y() == polycone.fTwoDVec[nextSegIndex].y()))
0459 distance = Real_v(kInfLength);
0460 else
0461 distance = DistanceToSeg(twoDPoint, polycone.fTwoDVec[currSegIndex], polycone.fTwoDVec[nextSegIndex]);
0462
0463 vecCore__MaskedAssignFunc(safety, (distance < safety) && !done, distance);
0464 }
0465
0466 if (polycone.fDeltaPhi < 2 * kPi) {
0467 Real_v safetyPhi = polycone.fPhiWedge.SafetyToOut<Real_v>(point);
0468 vecCore__MaskedAssignFunc(safety, safetyPhi < safety, safetyPhi);
0469 }
0470 }
0471
0472 #if (0)
0473
0474 template <typename Real_v>
0475 VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut_Old(UnplacedStruct_t const &polycone,
0476 Vector3D<Real_v> const &point,
0477 Real_v &safety)
0478 {
0479 typedef typename vecCore::Mask_v<Real_v> Bool_v;
0480 Bool_v compIn(false), compOut(false);
0481 GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, compIn, compOut);
0482 if (compOut) {
0483 safety = -1;
0484 return;
0485 }
0486
0487 int index = polycone.GetSectionIndex(point.z());
0488 if (index < 0) {
0489 safety = -1;
0490 return;
0491 }
0492
0493 PolyconeSection const &sec = polycone.GetSection(index);
0494
0495 Vector3D<Real_v> p = point - Vector3D<Precision>(0, 0, sec.fShift);
0496
0497 ConeImplementation<polyconeTypeT>::template SafetyToOut<Real_v>(sec.fSolid, p, safety);
0498
0499 Precision minSafety = safety;
0500 if (minSafety == kInfLength) {
0501 safety = 0.;
0502 return;
0503 }
0504 if (minSafety < kTolerance) {
0505 safety = 0.;
0506 return;
0507 }
0508
0509 Precision zbase = polycone.fZs[index + 1];
0510 for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0511 Precision dz = polycone.fZs[i] - zbase;
0512 if (dz >= minSafety) break;
0513 PolyconeSection const § = polycone.GetSection(i);
0514 p = point - Vector3D<Precision>(0, 0, sect.fShift);
0515
0516
0517 ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sect.fSolid, p, safety);
0518
0519 if (safety < minSafety) minSafety = safety;
0520 }
0521
0522 if (index > 0) {
0523 zbase = polycone.fZs[index - 1];
0524 for (int i = index - 1; i >= 0; --i) {
0525 Precision dz = zbase - polycone.fZs[i];
0526 if (dz >= minSafety) break;
0527 PolyconeSection const § = polycone.GetSection(i);
0528 p = point - Vector3D<Precision>(0, 0, sect.fShift);
0529
0530
0531 ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(sect.fSolid, p, safety);
0532
0533 if (safety < minSafety) minSafety = safety;
0534 }
0535 }
0536
0537 safety = minSafety;
0538 return;
0539 }
0540 #endif
0541 };
0542 }
0543 }
0544
0545 #endif