Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-09 09:27:33 

0001 /// @file GenericPolyconeImplementation.h
0002 /// @author Raman Sehgal (raman.sehgal@cern.ch)
0003 
0004 #ifndef VECGEOM_VOLUMES_KERNEL_GENERICPOLYCONEIMPLEMENTATION_H_
0005 #define VECGEOM_VOLUMES_KERNEL_GENERICPOLYCONEIMPLEMENTATION_H_
0006 
0007 #include "VecGeom/base/Vector3D.h"
0008 #include "VecGeom/volumes/GenericPolyconeStruct.h"
0009 #include "VecGeom/volumes/kernel/GenericKernels.h"
0010 #include <VecCore/VecCore>
0011 #include "VecGeom/volumes/kernel/CoaxialConesImplementation.h"
0012 #include "VecGeom/volumes/kernel/shapetypes/ConeTypes.h"
0013 
0014 #include <cstdio>
0015 
0016 namespace vecgeom {
0017 
0018 VECGEOM_DEVICE_FORWARD_DECLARE(struct GenericPolyconeImplementation;);
0019 VECGEOM_DEVICE_DECLARE_CONV(struct, GenericPolyconeImplementation);
0020 
0021 inline namespace VECGEOM_IMPL_NAMESPACE {
0022 
0023 class PlacedGenericPolycone;
0024 template <typename T>
0025 struct GenericPolyconeStruct;
0026 class UnplacedGenericPolycone;
0027 
0028 struct GenericPolyconeImplementation {
0029 
0030   using PlacedShape_t    = PlacedGenericPolycone;
0031   using UnplacedStruct_t = GenericPolyconeStruct<Precision>;
0032   using UnplacedVolume_t = UnplacedGenericPolycone;
0033 
0034   template <typename Real_v, bool ForInside>
0035   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForASection(
0036       UnplacedStruct_t const &unplaced, int isect, Vector3D<Real_v> const &polyconePoint,
0037       typename vecCore::Mask_v<Real_v> &secFullyInside, typename vecCore::Mask_v<Real_v> &secFullyOutside)
0038   {
0039 
0040     using namespace ConeTypes;
0041 
0042     if (isect < 0) {
0043       secFullyInside  = false;
0044       secFullyOutside = true;
0045       return;
0046     }
0047 
0048     GenericPolyconeSection const &sec = unplaced.GetSection(isect);
0049     Vector3D<Precision> secLocalp     = polyconePoint - Vector3D<Precision>(0, 0, sec.fShift);
0050 #ifdef POLYCONEDEBUG
0051     std::cerr << " isect=" << isect << "/" << unplaced.GetNSections() << " secLocalP=" << secLocalp
0052               << ", secShift=" << sec.fShift << " sec.fSolid=" << sec.fSolid << std::endl;
0053     if (sec.fSolid) sec.fSolid->Print();
0054 #endif
0055 
0056     CoaxialConesImplementation::template GenericKernelForContainsAndInside<Real_v, typename vecCore::Mask_v<Real_v>,
0057                                                                            ForInside>(*sec.fCoaxialCones, secLocalp,
0058                                                                                       secFullyInside, secFullyOutside);
0059   }
0060 
0061   template <typename Real_v, typename Bool_v>
0062   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Contains(UnplacedStruct_t const &genericPolycone,
0063                                                                     Vector3D<Real_v> const &point, Bool_v &inside)
0064   {
0065     Bool_v unused(false), outside(false);
0066     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(genericPolycone, point, unused, outside);
0067     inside = !outside;
0068   }
0069 
0070   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0071   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0072   template <typename Real_v, typename Inside_t>
0073   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void Inside(UnplacedStruct_t const &genericPolycone,
0074                                                                   Vector3D<Real_v> const &point, Inside_t &inside)
0075   {
0076 
0077     using Bool_v       = vecCore::Mask_v<Real_v>;
0078     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0079     Bool_v completelyinside(false), completelyoutside(false);
0080     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(genericPolycone, point, completelyinside,
0081                                                             completelyoutside);
0082     inside = EInside::kSurface;
0083     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0084     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0085   }
0086 
0087   template <typename Real_v, typename Bool_v, bool ForInside>
0088   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void GenericKernelForContainsAndInside(
0089       UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &localPoint, Bool_v &completelyInside,
0090       Bool_v &completelyOutside)
0091   {
0092     /* TODO : Logic to check where the point is inside or not.
0093     **
0094     ** if ForInside is false then it will only check if the point is outside,
0095     ** and is used by Contains function
0096     **
0097     ** if ForInside is true then it will check whether the point is inside or outside,
0098     ** and if neither inside nor outside then it is on the surface.
0099     ** and is used by Inside function
0100     */
0101 
0102     typedef Bool_v Bool_t;
0103 
0104     int indexLow  = unplaced.GetSectionIndex(localPoint.z() - kTolerance);
0105     int indexHigh = unplaced.GetSectionIndex(localPoint.z() + kTolerance);
0106     if (indexLow < 0 && indexHigh < 0) {
0107       completelyOutside = true;
0108       return;
0109     }
0110     if (indexLow < 0 && indexHigh == 0) {
0111       // Check location in section 0 and return
0112       GenericKernelForASection<Real_v, ForInside>(unplaced, 0, localPoint, completelyInside, completelyOutside);
0113       return;
0114     }
0115     if (indexHigh < 0 && indexLow == (unplaced.GetNSections() - 1)) {
0116       // Check location in section N-1 and return
0117       GenericKernelForASection<Real_v, ForInside>(unplaced, (unplaced.GetNSections() - 1), localPoint, completelyInside,
0118                                                   completelyOutside);
0119 
0120       return;
0121     }
0122     if (indexLow >= 0 && indexHigh >= 0) {
0123       if (indexLow == indexHigh) {
0124         // Check location in section indexLow and return
0125         GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, completelyInside,
0126                                                     completelyOutside);
0127 
0128         return;
0129       } else {
0130 
0131         Bool_t secInLow = false, secOutLow = false;
0132         Bool_t secInHigh = false, secOutHigh = false;
0133 
0134         GenericPolyconeSection const &sectionLow  = unplaced.GetSection(indexLow);
0135         GenericPolyconeSection const &sectionHigh = unplaced.GetSection(indexHigh);
0136 
0137         Bool_t onRing(false);
0138         onRing |= (CoaxialConesImplementation::template IsOnRing<Real_v, false>(
0139                        *sectionLow.fCoaxialCones, localPoint - Vector3D<Precision>(0, 0, sectionLow.fShift)) ||
0140                    CoaxialConesImplementation::template IsOnRing<Real_v, true>(
0141                        *sectionHigh.fCoaxialCones, localPoint - Vector3D<Precision>(0, 0, sectionHigh.fShift)));
0142 
0143         GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, secInLow, secOutLow);
0144         GenericKernelForASection<Real_v, ForInside>(unplaced, indexHigh, localPoint, secInHigh, secOutHigh);
0145         Bool_t surfLow  = !secInLow && !secOutLow;
0146         Bool_t surfHigh = !secInHigh && !secOutHigh;
0147 
0148         if (surfLow && surfHigh) {
0149           completelyInside = !onRing; // true;
0150           return;
0151         } else {
0152           // else if point is on surface of only one of the two sections then point is actually on surface , the default
0153           // case,
0154           // so no need to check
0155 
0156           // What needs to check is if it is outside both ie. Outside indexLow section and Outside indexHigh section
0157           // then it is certainly outside
0158           if (secOutLow && secOutHigh) {
0159             completelyOutside = true;
0160             return;
0161           }
0162         }
0163       }
0164     }
0165   }
0166 
0167   template <typename Real_v>
0168   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToIn(UnplacedStruct_t const &polycone,
0169                                                                         Vector3D<Real_v> const &point,
0170                                                                         Vector3D<Real_v> const &direction,
0171                                                                         Real_v const &stepMax, Real_v &distance)
0172   {
0173     // using namespace PolyconeTypes;
0174     Vector3D<Real_v> p = point;
0175     Vector3D<Real_v> v = direction;
0176 
0177 #ifdef POLYCONEDEBUG
0178     std::cerr << "Polycone::DistToIn() (spot 1): point=" << point << ", dir=" << direction << ", localPoint=" << p
0179               << ", localDir=" << v << "\n";
0180 #endif
0181 
0182     // TODO: add bounding box check maybe??
0183 
0184     distance      = kInfLength;
0185     int increment = (v.z() > 0) ? 1 : -1;
0186     if (std::fabs(v.z()) < kTolerance) increment = 0;
0187     int index = polycone.GetSectionIndex(p.z());
0188     if (index == -1) index = 0;
0189     if (index == -2) index = polycone.GetNSections() - 1;
0190 
0191     do {
0192       // now we have to find a section
0193       GenericPolyconeSection const &sec = polycone.GetSection(index);
0194 
0195 #ifdef POLYCONEDEBUG
0196       std::cerr << "Polycone::DistToIn() (spot 2):"
0197                 << " index=" << index << " NSec=" << polycone.GetNSections() << " &sec=" << &sec << " - secPars:"
0198                 << " secOffset=" << sec.fShift << " Dz=" << sec.fSolid->GetDz() << " Rmin1=" << sec.fSolid->GetRmin1()
0199                 << " Rmin2=" << sec.fSolid->GetRmin2() << " Rmax1=" << sec.fSolid->GetRmax1()
0200                 << " Rmax2=" << sec.fSolid->GetRmax2() << " -- calling Cone::DistToIn()...\n";
0201 #endif
0202 
0203       CoaxialConesImplementation::template DistanceToIn<Real_v>(
0204           *sec.fCoaxialCones, p - Vector3D<Precision>(0, 0, sec.fShift), v, stepMax, distance);
0205 
0206 #ifdef POLYCONEDEBUG
0207       std::cerr << "Polycone::DistToIn() (spot 3):"
0208                 << " distToIn() = " << distance << "\n";
0209 #endif
0210 
0211       if (distance < kInfLength || !increment) break;
0212       index += increment;
0213     } while (index >= 0 && index < polycone.GetNSections());
0214     return;
0215   }
0216 
0217   template <typename Real_v>
0218   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void DistanceToOut(UnplacedStruct_t const &polycone,
0219                                                                          Vector3D<Real_v> const &point,
0220                                                                          Vector3D<Real_v> const &dir,
0221                                                                          Real_v const &stepMax, Real_v &distance)
0222   {
0223     distance            = kInfLength;
0224     Vector3D<Real_v> pn = point;
0225     Precision dist      = 0.;
0226     int increment       = (dir.z() > 0) ? 1 : -1;
0227 
0228     // specialization for N==1??? It should be a cone in the first place
0229     if (polycone.GetNSections() == 1) {
0230       const GenericPolyconeSection &section = polycone.GetSection(0);
0231 
0232       CoaxialConesImplementation::template DistanceToOut<Real_v>(
0233           *section.fCoaxialCones, point - Vector3D<Precision>(0, 0, section.fShift), dir, stepMax, distance);
0234 
0235       return;
0236     }
0237 
0238     int indexLow  = polycone.GetSectionIndex(point.z() - kTolerance);
0239     int indexHigh = polycone.GetSectionIndex(point.z() + kTolerance);
0240     int index     = 0;
0241     if (indexLow < 0 && indexHigh < 0) {
0242       distance = -1;
0243       return;
0244     }
0245 
0246     if (indexLow < 0 && indexHigh >= 0 && dir.z() < 0.) {
0247       index                                 = indexHigh;
0248       const GenericPolyconeSection &section = polycone.GetSection(index);
0249       CoaxialConesImplementation::template DistanceToOut<Real_v>(*section.fCoaxialCones, pn, dir, stepMax, dist);
0250       distance = dist;
0251       return;
0252     }
0253 
0254     if (indexLow > 0 && indexHigh < 0 && dir.z() > 0.) {
0255       index                                 = indexLow;
0256       const GenericPolyconeSection &section = polycone.GetSection(index);
0257       CoaxialConesImplementation::template DistanceToOut<Real_v>(*section.fCoaxialCones, pn, dir, stepMax, dist);
0258       distance = dist;
0259       return;
0260     }
0261 
0262     Inside_t inside;
0263     Precision totalDistance = 0.;
0264     int count               = 0;
0265     do {
0266 
0267       if (indexLow >= 0 && indexHigh >= 0) {
0268         count++;
0269         index                                 = indexLow;
0270         const GenericPolyconeSection &section = polycone.GetSection(index);
0271         pn.z() -= section.fShift;
0272         CoaxialConesImplementation::template Inside<Real_v>(*section.fCoaxialCones, pn, inside);
0273         if (inside == EInside::kOutside) {
0274           if (count == 1) {
0275             distance = -1;
0276             return;
0277           } else {
0278             distance = totalDistance;
0279             return;
0280           }
0281         } else {
0282           CoaxialConesImplementation::template DistanceToOut<Real_v>(*section.fCoaxialCones, pn, dir, stepMax, dist);
0283           if (dist < 0.) break;
0284           totalDistance += dist;
0285           pn += dir * dist;
0286           pn.z() += section.fShift;
0287         }
0288         indexLow += increment;
0289         indexHigh += increment;
0290       }
0291     } while (indexLow > -1 && indexLow < polycone.GetNSections()); // end of do-while
0292     distance = totalDistance;
0293     return;
0294   }
0295 
0296   template <typename Real_v>
0297   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToIn(UnplacedStruct_t const &polycone,
0298                                                                       Vector3D<Real_v> const &point, Real_v &safety)
0299   {
0300 
0301     Vector3D<Real_v> p = point;
0302     int index          = polycone.GetSectionIndex(p.z());
0303 
0304     bool needZ = false;
0305     if (index < 0) {
0306       needZ = true;
0307       if (index == -1) index = 0;
0308       if (index == -2) index = polycone.GetNSections() - 1;
0309     }
0310     Precision minSafety               = 0; //= SafetyFromOutsideSection(index, p);
0311     GenericPolyconeSection const &sec = polycone.GetSection(index);
0312     // safety to current segment
0313     if (needZ) {
0314       CoaxialConesImplementation::template SafetyToIn<Real_v>(*sec.fCoaxialCones,
0315                                                               p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0316     } else
0317 
0318       CoaxialConesImplementation::template SafetyToIn<Real_v>(*sec.fCoaxialCones,
0319                                                               p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0320 
0321     if (safety < kTolerance) return;
0322     minSafety       = safety;
0323     Precision zbase = polycone.fZs[index + 1];
0324     // going right
0325     for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0326       Precision dz = polycone.fZs[i] - zbase;
0327       if (dz >= minSafety) break;
0328 
0329       GenericPolyconeSection const &sect = polycone.GetSection(i);
0330       CoaxialConesImplementation::template SafetyToIn<Real_v>(*sect.fCoaxialCones,
0331                                                               p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0332       if (safety < minSafety) minSafety = safety;
0333     }
0334 
0335     // going left if this is possible
0336     if (index > 0) {
0337       zbase = polycone.fZs[index - 1];
0338       for (int i = index - 1; i >= 0; --i) {
0339         Precision dz = zbase - polycone.fZs[i];
0340         if (dz >= minSafety) break;
0341         GenericPolyconeSection const &sect = polycone.GetSection(i);
0342 
0343         CoaxialConesImplementation::template SafetyToIn<Real_v>(*sect.fCoaxialCones,
0344                                                                 p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0345 
0346         if (safety < minSafety) minSafety = safety;
0347       }
0348     }
0349     safety = minSafety;
0350 
0351     return;
0352   }
0353 
0354   template <typename Real_v>
0355   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut(UnplacedStruct_t const &polycone,
0356                                                                        Vector3D<Real_v> const &point, Real_v &safety)
0357   {
0358     typedef typename vecCore::Mask_v<Real_v> Bool_v;
0359     Bool_v compIn(false), compOut(false);
0360     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, compIn, compOut);
0361     if (compOut) {
0362       safety = -1;
0363       return;
0364     }
0365 
0366     if (!compIn && !compOut) {
0367       safety = 0.;
0368       return;
0369     }
0370 
0371     int index = polycone.GetSectionIndex(point.z());
0372     if (index < 0) {
0373       safety = -1;
0374       return;
0375     }
0376 
0377     GenericPolyconeSection const &sec = polycone.GetSection(index);
0378 
0379     Vector3D<Real_v> p = point - Vector3D<Precision>(0, 0, sec.fShift);
0380     CoaxialConesImplementation::template SafetyToOut<Real_v>(*sec.fCoaxialCones, p, safety);
0381 
0382     Precision minSafety = safety;
0383     if (minSafety == kInfLength) {
0384       safety = 0.;
0385       return;
0386     }
0387     if (minSafety < kTolerance) {
0388       safety = 0.;
0389       return;
0390     }
0391 
0392     Precision zbase = polycone.fZs[index + 1];
0393     for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0394       Precision dz = polycone.fZs[i] - zbase;
0395       if (dz >= minSafety) break;
0396       GenericPolyconeSection const &sect = polycone.GetSection(i);
0397       p                                  = point - Vector3D<Precision>(0, 0, sect.fShift);
0398 
0399       CoaxialConesImplementation::template SafetyToIn<Real_v>(*sect.fCoaxialCones, p, safety);
0400 
0401       if (safety < minSafety) minSafety = safety;
0402     }
0403 
0404     if (index > 0) {
0405       zbase = polycone.fZs[index - 1];
0406       for (int i = index - 1; i >= 0; --i) {
0407         Precision dz = zbase - polycone.fZs[i];
0408         if (dz >= minSafety) break;
0409         GenericPolyconeSection const &sect = polycone.GetSection(i);
0410         p                                  = point - Vector3D<Precision>(0, 0, sect.fShift);
0411 
0412         CoaxialConesImplementation::template SafetyToIn<Real_v>(*sect.fCoaxialCones, p, safety);
0413 
0414         if (safety < minSafety) minSafety = safety;
0415       }
0416     }
0417 
0418     safety = minSafety;
0419     return;
0420   }
0421 };
0422 } // namespace VECGEOM_IMPL_NAMESPACE
0423 } // namespace vecgeom
0424 
0425 #endif // VECGEOM_VOLUMES_KERNEL_GENERICPOLYCONEIMPLEMENTATION_H_
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!