Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 10:14:01 

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