Back to home page

EIC code displayed by LXR
 
 

    

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

0001 /*
0002  * PolyconeImplementation.h
0003  *
0004  *  Created on: Dec 8, 2014
0005  *      Author: swenzel
0006  */
0007 
0008 /// History notes:
0009 /// Jan-March 2017: revision + moving to use new Cone Kernels (Raman Sehgal)
0010 /// May-June 2017: revision + moving to new Structure (Raman Sehgal)
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   VECCORE_ATT_HOST_DEVICE
0043   static void PrintType()
0044   {
0045     //  printf("SpecializedPolycone<%i, %i>", transCodeT, rotCodeT);
0046   }
0047 
0048   template <typename Stream>
0049   static void PrintType(Stream &st, int transCodeT = translation::kGeneric, int rotCodeT = rotation::kGeneric)
0050   {
0051     st << "SpecializedPolycone<" << transCodeT << "," << rotCodeT << ">";
0052   }
0053 
0054   template <typename Stream>
0055   static void PrintImplementationType(Stream &st)
0056   {
0057     (void)st;
0058     // st << "PolyconeImplementation<" << transCodeT << "," << rotCodeT << ">";
0059   }
0060 
0061   template <typename Stream>
0062   static void PrintUnplacedType(Stream &st)
0063   {
0064     (void)st;
0065     // TODO: this is wrong
0066     // st << "UnplacedPolycone";
0067   }
0068 
0069   template <typename Real_v, bool ForInside>
0070   VECGEOM_FORCE_INLINE
0071   VECCORE_ATT_HOST_DEVICE
0072   static void GenericKernelForASection(UnplacedStruct_t const &unplaced, int isect,
0073                                        Vector3D<Real_v> const &polyconePoint,
0074                                        typename vecCore::Mask_v<Real_v> &secFullyInside,
0075                                        typename vecCore::Mask_v<Real_v> &secFullyOutside)
0076   {
0077 
0078     // using namespace PolyconeTypes;
0079     using namespace ConeTypes;
0080 
0081     if (isect < 0) {
0082       secFullyInside  = false;
0083       secFullyOutside = true;
0084       return;
0085     }
0086 
0087     PolyconeSection const &sec    = unplaced.GetSection(isect);
0088     Vector3D<Precision> secLocalp = polyconePoint - Vector3D<Precision>(0, 0, sec.fShift);
0089 #ifdef POLYCONEDEBUG
0090     std::cout << " isect=" << isect << "/" << unplaced.GetNSections() << " secLocalP=" << secLocalp
0091               << ", secShift=" << sec.fShift << " sec.fSolid=" << sec.fSolid << std::endl;
0092     if (sec.fSolid) sec.fSolid->Print();
0093 #endif
0094 
0095     ConeHelpers<Real_v, polyconeTypeT>::template GenericKernelForContainsAndInside<ForInside>(
0096         *sec.fSolid, secLocalp, secFullyInside, secFullyOutside);
0097   }
0098 
0099   template <typename Real_v, typename Bool_v>
0100   VECGEOM_FORCE_INLINE
0101   VECCORE_ATT_HOST_DEVICE
0102   static void Contains(UnplacedStruct_t const &polycone, Vector3D<Real_v> const &point, Bool_v &inside)
0103   {
0104 
0105     Bool_v unused(false), outside(false);
0106     GenericKernelForContainsAndInside<Real_v, Bool_v, false>(polycone, point, unused, outside);
0107     inside = !outside;
0108   }
0109 
0110   // BIG QUESTION: DO WE WANT TO GIVE ALL 3 TEMPLATE PARAMETERS
0111   // -- OR -- DO WE WANT TO DEDUCE Bool_v, Index_t from Real_v???
0112   template <typename Real_v, typename Inside_t>
0113   VECGEOM_FORCE_INLINE
0114   VECCORE_ATT_HOST_DEVICE
0115   static void Inside(UnplacedStruct_t const &polycone, Vector3D<Real_v> const &point, Inside_t &inside)
0116   {
0117 
0118     using Bool_v       = vecCore::Mask_v<Real_v>;
0119     using InsideBool_v = vecCore::Mask_v<Inside_t>;
0120     Bool_v completelyinside(false), completelyoutside(false);
0121     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, completelyinside, completelyoutside);
0122     inside = EInside::kSurface;
0123     vecCore::MaskedAssign(inside, (InsideBool_v)completelyoutside, Inside_t(EInside::kOutside));
0124     vecCore::MaskedAssign(inside, (InsideBool_v)completelyinside, Inside_t(EInside::kInside));
0125   }
0126 
0127   template <typename Real_v, typename Bool_v, bool ForInside>
0128   VECGEOM_FORCE_INLINE
0129   VECCORE_ATT_HOST_DEVICE
0130   static void GenericKernelForContainsAndInside(UnplacedStruct_t const &unplaced, Vector3D<Real_v> const &localPoint,
0131                                                 Bool_v &completelyInside, Bool_v &completelyOutside)
0132   {
0133 
0134     typedef Bool_v Bool_t;
0135 
0136     int indexLow  = unplaced.GetSectionIndex(localPoint.z() - kTolerance);
0137     int indexHigh = unplaced.GetSectionIndex(localPoint.z() + kTolerance);
0138     if (indexLow < 0 && indexHigh < 0) {
0139       completelyOutside = true;
0140       return;
0141     }
0142     if (indexLow < 0 && indexHigh == 0) {
0143       // Check location in section 0 and return
0144       GenericKernelForASection<Real_v, ForInside>(unplaced, 0, localPoint, completelyInside, completelyOutside);
0145       return;
0146     }
0147     if (indexHigh < 0 && indexLow == (unplaced.GetNSections() - 1)) {
0148       // Check location in section N-1 and return
0149       GenericKernelForASection<Real_v, ForInside>(unplaced, (unplaced.GetNSections() - 1), localPoint, completelyInside,
0150                                                   completelyOutside);
0151 
0152       return;
0153     }
0154     if (indexLow >= 0 && indexHigh >= 0) {
0155       if (indexLow == indexHigh) {
0156         // Check location in section indexLow and return
0157         GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, completelyInside,
0158                                                     completelyOutside);
0159 
0160         return;
0161       } else {
0162 
0163         Bool_t secInLow = false, secOutLow = false;
0164         Bool_t secInHigh = false, secOutHigh = false;
0165         GenericKernelForASection<Real_v, ForInside>(unplaced, indexLow, localPoint, secInLow, secOutLow);
0166         GenericKernelForASection<Real_v, ForInside>(unplaced, indexHigh, localPoint, secInHigh, secOutHigh);
0167         Bool_t surfLow  = !secInLow && !secOutLow;
0168         Bool_t surfHigh = !secInHigh && !secOutHigh;
0169 
0170         if (surfLow && surfHigh) {
0171           completelyInside = true;
0172           return;
0173         } else {
0174           // else if point is on surface of only one of the two sections then point is actually on surface , the default
0175           // case,
0176           // so no need to check
0177 
0178           // What needs to check is if it is outside both ie. Outside indexLow section and Outside indexHigh section
0179           // then it is certainly outside
0180           if (secOutLow && secOutHigh) {
0181             completelyOutside = true;
0182             return;
0183           }
0184         }
0185       }
0186     }
0187   }
0188 
0189   template <typename Real_v>
0190   VECGEOM_FORCE_INLINE
0191   VECCORE_ATT_HOST_DEVICE
0192   static void DistanceToIn(UnplacedStruct_t const &polycone, Vector3D<Real_v> const &point,
0193                            Vector3D<Real_v> const &direction, Real_v const &stepMax, Real_v &distance)
0194   {
0195     // using namespace PolyconeTypes;
0196     Vector3D<Real_v> p = point;
0197     Vector3D<Real_v> v = direction;
0198 
0199 #ifdef POLYCONEDEBUG
0200     std::cout << "Polycone::DistToIn() (spot 1): point=" << point << ", dir=" << direction << ", localPoint=" << p
0201               << ", localDir=" << v << "\n";
0202 #endif
0203 
0204     // TODO: add bounding box check maybe??
0205 
0206     distance      = kInfLength;
0207     int increment = (v.z() > 0) ? 1 : -1;
0208     if (std::fabs(v.z()) < kTolerance) increment = 0;
0209     int index = polycone.GetSectionIndex(p.z());
0210     if (index == -1) index = 0;
0211     if (index == -2) index = polycone.GetNSections() - 1;
0212 
0213     do {
0214       // now we have to find a section
0215       PolyconeSection const &sec = polycone.GetSection(index);
0216 
0217 #ifdef POLYCONEDEBUG
0218       std::cout << "Polycone::DistToIn() (spot 2):"
0219                 << " index=" << index << " NSec=" << polycone.GetNSections() << " &sec=" << &sec << " - secPars:"
0220                 << " secOffset=" << sec.fShift << " Dz=" << sec.fSolid->GetDz() << " Rmin1=" << sec.fSolid->GetRmin1()
0221                 << " Rmin2=" << sec.fSolid->GetRmin2() << " Rmax1=" << sec.fSolid->GetRmax1()
0222                 << " Rmax2=" << sec.fSolid->GetRmax2() << " -- calling Cone::DistToIn()...\n";
0223 #endif
0224 
0225       ConeImplementation<polyconeTypeT>::template DistanceToIn<Real_v>(
0226           *sec.fSolid, p - Vector3D<Precision>(0, 0, sec.fShift), v, stepMax, distance);
0227 
0228 #ifdef POLYCONEDEBUG
0229       std::cout << "Polycone::DistToIn() (spot 3):"
0230                 << " distToIn() = " << distance << "\n";
0231 #endif
0232 
0233       if (distance < kInfLength || !increment) break;
0234       index += increment;
0235     } while (index >= 0 && index < polycone.GetNSections());
0236     return;
0237   }
0238 
0239   template <typename Real_v>
0240   VECGEOM_FORCE_INLINE
0241   VECCORE_ATT_HOST_DEVICE
0242   static void DistanceToOut(UnplacedStruct_t const &polycone, Vector3D<Real_v> const &point,
0243                             Vector3D<Real_v> const &dir, Real_v const &stepMax, Real_v &distance)
0244   {
0245     // using namespace PolyconeTypes;
0246     distance            = kInfLength;
0247     Vector3D<Real_v> pn = point;
0248 
0249     // specialization for N==1??? It should be a cone in the first place
0250     if (polycone.GetNSections() == 1) {
0251       const PolyconeSection &section = polycone.GetSection(0);
0252 
0253       ConeImplementation<polyconeTypeT>::template DistanceToOut<Real_v>(
0254           *section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), dir, stepMax, distance);
0255 
0256       return;
0257     }
0258 
0259     int indexLow  = polycone.GetSectionIndex(point.z() - kTolerance);
0260     int indexHigh = polycone.GetSectionIndex(point.z() + kTolerance);
0261     int index     = 0;
0262 
0263     // section index is -1 when out of left-end
0264     // section index is -2 when beyond right-end
0265 
0266     if (indexLow < 0 && indexHigh < 0) {
0267       distance = -1;
0268       return;
0269     } else if (indexLow < 0 && indexHigh >= 0) {
0270       index                          = indexHigh;
0271       const PolyconeSection &section = polycone.GetSection(index);
0272 
0273       Inside_t inside;
0274       //      ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(
0275       //          *section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0276       ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0277           *section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0278       if (inside == EInside::kOutside) {
0279         distance = -1;
0280         return;
0281       }
0282     } else if (indexLow != indexHigh && (indexLow >= 0)) {
0283       // we are close to an intermediate Surface, section has to be identified
0284       const PolyconeSection &section = polycone.GetSection(indexLow);
0285 
0286       Inside_t inside;
0287       //      ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(
0288       //        *section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0289 
0290       ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0291           *section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0292 
0293       if (inside == EInside::kOutside) {
0294         index = indexHigh;
0295       } else {
0296         index = indexLow;
0297       }
0298     } else {
0299       index = indexLow;
0300       if (index < 0) index = polycone.GetSectionIndex(point.z());
0301     }
0302     if (index < 0) {
0303       distance = 0.;
0304       return;
0305     }
0306     // Added
0307     else {
0308       const PolyconeSection &section = polycone.GetSection(index);
0309 
0310       Inside_t inside;
0311       //      ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(
0312       //          *section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0313       ConeImplementation<polyconeTypeT>::template Inside<Real_v>(
0314           *section.fSolid, point - Vector3D<Precision>(0, 0, section.fShift), inside);
0315       if (inside == EInside::kOutside) {
0316         distance = -1;
0317         return;
0318       }
0319     }
0320 
0321     Precision totalDistance = 0.;
0322     Precision dist;
0323     int increment = (dir.z() > 0) ? 1 : -1;
0324     if (std::fabs(dir.z()) < kTolerance) increment = 0;
0325 
0326     // What is the relevance of istep?
0327     int istep = 0;
0328     do {
0329       const PolyconeSection &section = polycone.GetSection(index);
0330 
0331       if ((totalDistance != 0) || (istep < 2)) {
0332         pn = point + totalDistance * dir; // point must be shifted, so it could eventually get into another solid
0333         pn.z() -= section.fShift;
0334         Inside_t inside;
0335         //        ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(*section.fSolid, pn, inside);
0336         ConeImplementation<polyconeTypeT>::template Inside<Real_v>(*section.fSolid, pn, inside);
0337 
0338         if (inside == EInside::kOutside) {
0339           break;
0340         }
0341       } else
0342         pn.z() -= section.fShift;
0343 
0344       istep++;
0345 
0346       // ConeImplementation<ConeTypes::UniversalCone>::DistanceToOut<Real_v>(*section.fSolid, pn, dir, stepMax, dist);
0347       ConeImplementation<polyconeTypeT>::template DistanceToOut<Real_v>(*section.fSolid, pn, dir, stepMax, dist);
0348       if (dist == -1) return;
0349 
0350       // Section Surface case
0351       if (std::fabs(dist) < 0.5 * kTolerance) {
0352         int index1 = index;
0353         if ((index > 0) && (index < polycone.GetNSections() - 1)) {
0354           index1 += increment;
0355         } else {
0356           if ((index == 0) && (increment > 0)) index1 += increment;
0357           if ((index == polycone.GetNSections() - 1) && (increment < 0)) index1 += increment;
0358         }
0359 
0360         Vector3D<Precision> pte         = point + (totalDistance + dist) * dir;
0361         const PolyconeSection &section1 = polycone.GetSection(index1);
0362         pte.z() -= section1.fShift;
0363         Vector3D<Precision> localp;
0364         Inside_t inside22;
0365         // ConeImplementation<ConeTypes::UniversalCone>::Inside<Real_v>(*section1.fSolid, pte, inside22);
0366         ConeImplementation<polyconeTypeT>::template Inside<Real_v>(*section1.fSolid, pte, inside22);
0367         if (inside22 == 3 || (increment == 0)) {
0368           break;
0369         }
0370       } // end if surface case
0371 
0372       totalDistance += dist;
0373       index += increment;
0374     } while (increment != 0 && index >= 0 && index < polycone.GetNSections());
0375 
0376     distance = totalDistance;
0377 
0378     return;
0379   }
0380 
0381   template <typename Real_v>
0382   VECGEOM_FORCE_INLINE
0383   VECCORE_ATT_HOST_DEVICE
0384   static void SafetyToIn(UnplacedStruct_t const &polycone, Vector3D<Real_v> const &point, Real_v &safety)
0385   {
0386 
0387     Vector3D<Real_v> p = point;
0388     int index          = polycone.GetSectionIndex(p.z());
0389 
0390     bool needZ = false;
0391     if (index < 0) {
0392       needZ = true;
0393       if (index == -1) index = 0;
0394       if (index == -2) index = polycone.GetNSections() - 1;
0395     }
0396     Precision minSafety        = 0; //= SafetyFromOutsideSection(index, p);
0397     PolyconeSection const &sec = polycone.GetSection(index);
0398     // safety to current segment
0399     if (needZ) {
0400       //      ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0401       //          *sec.fSolid, p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0402       ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(*sec.fSolid,
0403                                                                      p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0404     } else {
0405 
0406       //      ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0407       //          *sec.fSolid, p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0408       ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(*sec.fSolid,
0409                                                                      p - Vector3D<Precision>(0, 0, sec.fShift), safety);
0410 
0411       if (safety < kTolerance) return;
0412       minSafety       = safety;
0413       Precision zbase = polycone.fZs[index + 1];
0414       // going right
0415       for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0416         Precision dz = polycone.fZs[i] - zbase;
0417         if (dz >= minSafety) break;
0418 
0419         PolyconeSection const &sect = polycone.GetSection(i);
0420 
0421         //      ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0422         //          *sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0423 
0424         ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(
0425             *sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0426 
0427         if (safety < minSafety) minSafety = safety;
0428       }
0429 
0430       // going left if this is possible
0431       if (index > 0) {
0432         zbase = polycone.fZs[index - 1];
0433         for (int i = index - 1; i >= 0; --i) {
0434           Precision dz = zbase - polycone.fZs[i];
0435           if (dz >= minSafety) break;
0436           PolyconeSection const &sect = polycone.GetSection(i);
0437 
0438           //        ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(
0439           //            *sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0440 
0441           ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(
0442               *sect.fSolid, p - Vector3D<Precision>(0, 0, sect.fShift), safety);
0443 
0444           if (safety < minSafety) minSafety = safety;
0445         }
0446       }
0447       safety = minSafety;
0448     }
0449     return;
0450   }
0451   /* A function to calculation the shortest distance of a point from a segment */
0452   template <typename Real_v>
0453   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static Real_v DistanceToSeg(Vector3D<Real_v> const &point,
0454                                                                            Vector3D<Precision> segment_start,
0455                                                                            Vector3D<Precision> segment_end)
0456   {
0457     Vector3D<Real_v> segment_vector   = segment_end - segment_start;
0458     Vector3D<Real_v> point_vector     = point - segment_start;
0459     Real_v projection_scalar          = point_vector.Dot(segment_vector) / segment_vector.Mag2();
0460     Vector3D<Real_v> projection_point = segment_start + projection_scalar * segment_vector;
0461     vecCore__MaskedAssignFunc(projection_point, projection_scalar < Real_v(0.), Vector3D<Real_v>(segment_start));
0462     vecCore__MaskedAssignFunc(projection_point, projection_scalar > Real_v(1.), Vector3D<Real_v>(segment_end));
0463     Vector3D<Real_v> distVec = point - projection_point;
0464     return distVec.Mag();
0465   }
0466 
0467   /*
0468   ** New definition that uses Segments and Single Wedge, instead of SafetyToOut from ConeImplementation
0469   */
0470   template <typename Real_v>
0471   VECGEOM_FORCE_INLINE
0472   VECCORE_ATT_HOST_DEVICE
0473   static void SafetyToOut(UnplacedStruct_t const &polycone, Vector3D<Real_v> const &point, Real_v &safety)
0474   {
0475     typedef typename vecCore::Mask_v<Real_v> Bool_v;
0476     safety = Real_v(kInfLength);
0477     Bool_v compIn(false), compOut(false), done(false);
0478     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, compIn, compOut);
0479     vecCore__MaskedAssignFunc(safety, compOut && !done, Real_v(-1.));
0480     done |= compOut;
0481     if (vecCore::MaskFull(done)) return;
0482 
0483     Vector3D<Real_v> twoDPoint = Vector3D<Real_v>(point.Perp(), point.z(), 0.);
0484     for (unsigned int currSegIndex = 0; currSegIndex < polycone.fTwoDVec.size(); currSegIndex++) {
0485       unsigned int nextSegIndex = currSegIndex + 1;
0486       if (currSegIndex == polycone.fTwoDVec.size() - 1) {
0487         nextSegIndex = 0;
0488       }
0489 
0490       Real_v distance(kInfLength);
0491       if ((polycone.fTwoDVec[currSegIndex].x() == 0 && polycone.fTwoDVec[nextSegIndex].x() == 0) ||
0492           (polycone.fTwoDVec[currSegIndex].x() == polycone.fTwoDVec[nextSegIndex].x() &&
0493            polycone.fTwoDVec[currSegIndex].y() == polycone.fTwoDVec[nextSegIndex].y()))
0494         distance = Real_v(kInfLength);
0495       else
0496         distance = DistanceToSeg(twoDPoint, polycone.fTwoDVec[currSegIndex], polycone.fTwoDVec[nextSegIndex]);
0497 
0498       vecCore__MaskedAssignFunc(safety, (distance < safety) && !done, distance);
0499     }
0500 
0501     if (polycone.fDeltaPhi < 2 * kPi) {
0502       Real_v safetyPhi = polycone.fPhiWedge.SafetyToOut<Real_v>(point);
0503       vecCore__MaskedAssignFunc(safety, safetyPhi < safety, safetyPhi);
0504     }
0505   }
0506 
0507 #if (0)
0508   /* Retaining the old definition for the time being */
0509   template <typename Real_v>
0510   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE static void SafetyToOut_Old(UnplacedStruct_t const &polycone,
0511                                                                            Vector3D<Real_v> const &point,
0512                                                                            Real_v &safety)
0513   {
0514     typedef typename vecCore::Mask_v<Real_v> Bool_v;
0515     Bool_v compIn(false), compOut(false);
0516     GenericKernelForContainsAndInside<Real_v, Bool_v, true>(polycone, point, compIn, compOut);
0517     if (compOut) {
0518       safety = -1;
0519       return;
0520     }
0521 
0522     int index = polycone.GetSectionIndex(point.z());
0523     if (index < 0) {
0524       safety = -1;
0525       return;
0526     }
0527 
0528     PolyconeSection const &sec = polycone.GetSection(index);
0529 
0530     Vector3D<Real_v> p = point - Vector3D<Precision>(0, 0, sec.fShift);
0531     // ConeImplementation<ConeTypes::UniversalCone>::SafetyToOut<Real_v>(*sec.fSolid, p, safety);
0532     ConeImplementation<polyconeTypeT>::template SafetyToOut<Real_v>(*sec.fSolid, p, safety);
0533 
0534     Precision minSafety = safety;
0535     if (minSafety == kInfLength) {
0536       safety = 0.;
0537       return;
0538     }
0539     if (minSafety < kTolerance) {
0540       safety = 0.;
0541       return;
0542     }
0543 
0544     Precision zbase = polycone.fZs[index + 1];
0545     for (int i = index + 1; i < polycone.GetNSections(); ++i) {
0546       Precision dz = polycone.fZs[i] - zbase;
0547       if (dz >= minSafety) break;
0548       PolyconeSection const &sect = polycone.GetSection(i);
0549       p                           = point - Vector3D<Precision>(0, 0, sect.fShift);
0550 
0551       // ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(*sect.fSolid, p, safety);
0552       ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(*sect.fSolid, p, safety);
0553 
0554       if (safety < minSafety) minSafety = safety;
0555     }
0556 
0557     if (index > 0) {
0558       zbase = polycone.fZs[index - 1];
0559       for (int i = index - 1; i >= 0; --i) {
0560         Precision dz = zbase - polycone.fZs[i];
0561         if (dz >= minSafety) break;
0562         PolyconeSection const &sect = polycone.GetSection(i);
0563         p                           = point - Vector3D<Precision>(0, 0, sect.fShift);
0564 
0565         // ConeImplementation<ConeTypes::UniversalCone>::SafetyToIn<Real_v>(*sect.fSolid, p, safety);
0566         ConeImplementation<polyconeTypeT>::template SafetyToIn<Real_v>(*sect.fSolid, p, safety);
0567 
0568         if (safety < minSafety) minSafety = safety;
0569       }
0570     }
0571 
0572     safety = minSafety;
0573     return;
0574   }
0575 #endif
0576 };
0577 } // namespace VECGEOM_IMPL_NAMESPACE
0578 } // namespace vecgeom
0579 
0580 #endif // VECGEOM_VOLUMES_KERNEL_polyconeIMPLEMENTATION_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!