VecGeom/navigation/VNavigator.h

0001 /*
0002  * VNavigator.h
0003  *
0004  *  Created on: 17.09.2015
0005  *      Author: swenzel
0006  */
0007
0008 #ifndef NAVIGATION_VNAVIGATOR_H_
0009 #define NAVIGATION_VNAVIGATOR_H_
0010
0011 #include "VecGeom/base/Global.h"
0012 #include "VecGeom/base/Vector3D.h"
0013 #include "VecGeom/base/SOA3D.h"
0014 #include "VecGeom/base/Transformation3D.h"
0015 #include "VecGeom/navigation/NavigationState.h"
0016 #include "VecGeom/navigation/GlobalLocator.h"
0017 #include "VecGeom/volumes/PlacedVolume.h"
0018 #include "VecGeom/volumes/LogicalVolume.h"
0019 #include "VecGeom/navigation/VSafetyEstimator.h"
0020 #include "VecGeom/navigation/NavStateFwd.h"
0021
0022 namespace vecgeom {
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 // some forward declarations
0026 template <typename T>
0027 class Vector3D;
0028 // class NavigationState;
0029 class LogicalVolume;
0030 class Transformation3D;
0031 class VPlacedVolume;
0032
0033 //! base class defining basic interface for navigation ( hit-detection )
0034 //! sub classes implement optimized algorithms for logical volumes
0035 class VNavigator {
0036
0037 public:
0038   VNavigator() : fSafetyEstimator(nullptr) {}
0039
0040   VECCORE_ATT_HOST_DEVICE
0041   VSafetyEstimator const *GetSafetyEstimator() const { return fSafetyEstimator; }
0042
0043   //! computes the step (distance) to the next object in the geometry hierarchy obtained
0044   //! by propagating with step along the ray
0045   //! the next object could be after a boundary and does not necessarily coincide with the object
0046   //! hit by the ray
0047
0048   //! this methods transforms the global coordinates into local ones usually calls more specialized methods
0049   //! like the hit detection on local coordinates
0050   VECCORE_ATT_HOST_DEVICE
0051   virtual Precision ComputeStepAndPropagatedState(Vector3D<Precision> const & /*globalpoint*/,
0052                                                   Vector3D<Precision> const & /*globaldir*/,
0053                                                   Precision /*(physics) step limit */,
0054                                                   NavigationState const & /*in_state*/,
0055                                                   NavigationState & /*out_state*/) const = 0;
0056
0057   //! computes the step (distance) to the next object in the geometry hierarchy obtained
0058   //! by propagating with step along the ray
0059   //! updates out_state to contain information about the next hitting boundary:
0060   //!   - if a daugher is hit: out_state.Top() will be daughter
0061   //!   - if ray leaves volume: out_state.Top() will point to current volume
0062   //!   - if step limit > step: out_state == in_state
0063   //!
0064   //! This function is essentialy equal to ComputeStepAndPropagatedState without
0065   //! the relocation part
0066   virtual Precision ComputeStep(Vector3D<Precision> const & /*globalpoint*/, Vector3D<Precision> const & /*globaldir*/,
0067                                 Precision /*(physics) step limit */, NavigationState const & /*in_state*/,
0068                                 NavigationState & /*out_state*/) const = 0;
0069
0070   //! as above ... also returns the safety ... does not give_back an out_state
0071   //! but the in_state might be modified to contain the next daughter when
0072   //! user specifies indicateDaughterHit = true
0073   VECCORE_ATT_HOST_DEVICE
0074   virtual Precision ComputeStepAndSafety(Vector3D<Precision> const & /*globalpoint*/,
0075                                          Vector3D<Precision> const & /*globaldir*/, Precision /*(physics) step limit */,
0076                                          NavigationState & /*in_state*/, bool /*calcsafety*/, Precision & /*safety*/,
0077                                          bool indicateDaughterHit = false) const = 0;
0078
0079   // an alias interface ( using TGeo name )
0080   VECCORE_ATT_HOST_DEVICE
0081   void FindNextBoundaryAndStep(Vector3D<Precision> const &globalpoint, Vector3D<Precision> const &globaldir,
0082                                NavigationState const &in_state, NavigationState &out_state, Precision step_limit,
0083                                Precision &step) const
0084   {
0085     step = ComputeStepAndPropagatedState(globalpoint, globaldir, step_limit, in_state, out_state);
0086   }
0087
0088   // an alias interface ( using TGeo name )
0089   void FindNextBoundaryAndStepAndSafety(Vector3D<Precision> const &globalpoint, Vector3D<Precision> const &globaldir,
0090                                         NavigationState const &in_state, NavigationState &out_state,
0091                                         Precision step_limit, Precision &step, bool calcsafety, Precision &safety) const
0092   {
0093     step = ComputeStepAndSafetyAndPropagatedState(globalpoint, globaldir, step_limit, in_state, out_state, calcsafety,
0094                                                   safety);
0095   }
0096
0097   // a similar interface, in addition also returning the safety as a result
0098   VECCORE_ATT_HOST_DEVICE
0099   virtual Precision ComputeStepAndSafetyAndPropagatedState(Vector3D<Precision> const & /*globalpoint*/,
0100                                                            Vector3D<Precision> const & /*globaldir*/,
0101                                                            Precision /*(physics) step limit */,
0102                                                            NavigationState const & /*in_state*/,
0103                                                            NavigationState & /*out_state*/, bool /*calcsafefty*/,
0104                                                            Precision & /*safety_out*/) const = 0;
0105
0106   // the bool return type indicates if out_state was already modified; this may happen in assemblies;
0107   // in this case we don't need to copy the in_state to the out state later on
0108   // NavigationState a pointer since we might want to pass nullptr
0109   VECCORE_ATT_HOST_DEVICE
0110   virtual bool CheckDaughterIntersections(LogicalVolume const * /*lvol*/, Vector3D<Precision> const & /*localpoint*/,
0111                                           Vector3D<Precision> const & /*localdir*/,
0112                                           NavigationState const * /*in_state*/, NavigationState * /*out_state*/,
0113                                           Precision & /*step*/, VPlacedVolume const *& /*hitcandidate*/) const = 0;
0114
0115   /// check if a ray given by localpoint, localdir intersects with any daughter. Possibility
0116   /// to pass a volume which is blocked/should be ignored in the query. Updates the step as well as the hitcandidate
0117   /// volume. (This version is useful for G4; assemblies not supported)
0118   VECCORE_ATT_HOST_DEVICE
0119   virtual bool CheckDaughterIntersections(LogicalVolume const * /*lvol*/, Vector3D<Precision> const & /*localpoint*/,
0120                                           Vector3D<Precision> const & /*localdir*/, VPlacedVolume const * /*blocked*/,
0121                                           Precision & /*step*/, VPlacedVolume const *& /*hitcandidate*/) const
0122   {
0123     assert(false); // Not implemented --- notify of failure !!
0124     return false;
0125   }
0126
0127   // interfaces for vector/basket navigation
0128   virtual void ComputeStepsAndPropagatedStates(SOA3D<Precision> const & /*globalpoints*/,
0129                                                SOA3D<Precision> const & /*globaldirs*/,
0130                                                Precision const * /*(physics) step limits */,
0131                                                NavigationState const *const * /*in_states*/,
0132                                                NavigationState ** /*out_states*/, Precision * /*out_steps*/) const = 0;
0133
0134   // for vector navigation
0135   VECCORE_ATT_HOST_DEVICE
0136   virtual void ComputeStepsAndSafetiesAndPropagatedStates(SOA3D<Precision> const &points, SOA3D<Precision> const &dirs,
0137                                                           Precision const *psteps, NavStatePool const &instates,
0138                                                           NavStatePool &outstates, Precision *outsteps,
0139                                                           bool const *calcsafety, Precision *safeties) const = 0;
0140
0141   VECCORE_ATT_HOST_DEVICE
0142   virtual void ComputeStepsAndSafetiesAndPropagatedStates(SOA3D<Precision> const & /*globalpoints*/,
0143                                                           SOA3D<Precision> const & /*globaldirs*/,
0144                                                           Precision const * /*(physics) step limits */,
0145                                                           NavigationState const *const * /*in_states*/,
0146                                                           NavigationState ** /*out_states*/, Precision * /*out_steps*/,
0147                                                           bool const * /*calcsafety*/,
0148                                                           Precision * /*out_safeties*/) const = 0;
0149
0150   // for vector navigation (not doing relocation) -- interface for GeantV
0151   VECCORE_ATT_HOST_DEVICE
0152   virtual void ComputeStepsAndSafeties(SOA3D<Precision> const & /*globalpoints*/,
0153                                        SOA3D<Precision> const & /*globaldirs*/,
0154                                        Precision const * /*(physics) step limits */,
0155                                        NavigationState const *const * /*in_states*/, Precision * /*out_steps*/,
0156                                        bool const * /*calcsafety*/, Precision * /*out_safeties*/) const = 0;
0157
0158 protected:
0159   // a common relocate method ( to calculate propagated states after the boundary )
0160   VECCORE_ATT_HOST_DEVICE
0161   virtual void Relocate(Vector3D<Precision> const & /*localpoint*/, NavigationState const &__restrict__ /*in_state*/,
0162                         NavigationState &__restrict__ /*out_state*/) const = 0;
0163
0164   // a common function to be used by all navigators to ensure consistency in transporting points
0165   // after a boundary
0166   VECGEOM_FORCE_INLINE
0167   VECCORE_ATT_HOST_DEVICE
0168   static Vector3D<Precision> MovePointAfterBoundary(Vector3D<Precision> const &localpoint,
0169                                                     Vector3D<Precision> const &dir, Precision step)
0170   {
0171     const Precision extra = 1E-6; // TODO: to be revisited (potentially going for a more relative approach)
0172     return localpoint + (step + extra) * dir;
0173   }
0174
0175 public:
0176   VECCORE_ATT_DEVICE
0177   virtual ~VNavigator(){};
0178
0179   // get name of implementing class
0180   virtual const char *GetName() const = 0;
0181
0182   typedef VSafetyEstimator SafetyEstimator_t;
0183
0184 protected:
0185   VECCORE_ATT_HOST_DEVICE
0186   VNavigator(VSafetyEstimator *s) : fSafetyEstimator(s) {}
0187   VSafetyEstimator *fSafetyEstimator; // a pointer to the safetyEstimator which can be used by the Navigator
0188
0189   // some common code to prepare the outstate
0190   VECGEOM_FORCE_INLINE
0191   VECCORE_ATT_HOST_DEVICE
0192   static Precision PrepareOutState(NavigationState const &__restrict__ in_state,
0193                                    NavigationState &__restrict__ out_state, Precision geom_step, Precision step_limit,
0194                                    VPlacedVolume const *hitcandidate, bool &doneafterthisstep)
0195   {
0196     // now we have the candidates and we prepare the out_state
0197     in_state.CopyTo(&out_state);
0198     doneafterthisstep = false;
0199
0200     // if the following is the case we are in the wrong volume;
0201     // assuming that DistanceToIn returns negative number when point is inside
0202     // do nothing (step=0) and retry one level higher
0203
0204     // TODO: put diagnostic code here ( like in original SimpleNavigator )
0205     if (geom_step == kInfLength && step_limit > 0.) {
0206       geom_step = vecgeom::kTolerance;
0207       out_state.SetBoundaryState(true);
0208       do {
0209         out_state.Pop();
0210       } while (out_state.Top()->GetLogicalVolume()->GetUnplacedVolume()->IsAssembly());
0211       doneafterthisstep = true;
0212       return geom_step;
0213     }
0214
0215     // is geometry further away than physics step?
0216     // this is a physics step
0217     if (geom_step > step_limit) {
0218       // don't need to do anything
0219       geom_step = step_limit;
0220       out_state.SetBoundaryState(false);
0221       return geom_step;
0222     }
0223
0224     // otherwise it is a geometry step
0225     out_state.SetBoundaryState(true);
0226     out_state.SetLastExited();
0227     if (hitcandidate) out_state.Push(hitcandidate);
0228
0229     if (geom_step < 0.) {
0230       // std::cerr << "WARNING: STEP NEGATIVE; NEXTVOLUME " << nexthitvolume << std::endl;
0231       // InspectEnvironmentForPointAndDirection( globalpoint, globaldir, currentstate );
0232       geom_step = 0.;
0233     }
0234     return geom_step;
0235   }
0236
0237   // kernel to be used with both scalar and vector types
0238   template <typename T>
0239   VECGEOM_FORCE_INLINE
0240   VECCORE_ATT_HOST_DEVICE
0241   static T TreatDistanceToMother(VPlacedVolume const *pvol, Vector3D<T> const &localpoint, Vector3D<T> const &localdir,
0242                                  T step_limit)
0243   {
0244     T step;
0245     assert(pvol != nullptr && "currentvolume is null in navigation");
0246     step = pvol->DistanceToOut(localpoint, localdir, step_limit);
0247     vecCore::MaskedAssign(step, step < T(0.), T(0.));
0248     return step;
0249   }
0250
0251   // default static function doing the global to local transformation
0252   // may be redefined in concrete implementations ( for instance in cases where we know the form of the global matrix
0253   // a-priori )
0254   // input
0255   // TODO: think about how we can have scalar + SIMD version
0256   // note: the last argument is a trick to pass information across function calls ( only exploited in specialized
0257   // navigators )
0258   template <typename T>
0259   VECGEOM_FORCE_INLINE
0260   VECCORE_ATT_HOST_DEVICE
0261   static void DoGlobalToLocalTransformation(NavigationState const &in_state, Vector3D<T> const &globalpoint,
0262                                             Vector3D<T> const &globaldir, Vector3D<T> &localpoint,
0263                                             Vector3D<T> &localdir)
0264   {
0265     // calculate local point/dir from global point/dir
0266     Transformation3D m;
0267     in_state.TopMatrix(m);
0268     localpoint = m.Transform(globalpoint);
0269     localdir   = m.TransformDirection(globaldir);
0270   }
0271
0272   // version used for SIMD processing
0273   // should be specialized in Impl for faster treatment
0274   template <typename T, unsigned int ChunkSize>
0275   VECGEOM_FORCE_INLINE
0276   VECCORE_ATT_HOST_DEVICE
0277   static void DoGlobalToLocalTransformations(NavigationState const *const *in_states,
0278                                              SOA3D<Precision> const &globalpoints, SOA3D<Precision> const &globaldirs,
0279                                              unsigned int from_index, Vector3D<T> &localpoint, Vector3D<T> &localdir)
0280   {
0281     for (unsigned int i = 0; i < ChunkSize; ++i) {
0282       unsigned int trackid = from_index + i;
0283       Transformation3D m;
0284       //  assert(in_states[trackid]->Top()->GetLogicalVolume() == lvol &&
0285       //         "not all states in same logical volume"); // the logical volume of all the states should be the same
0286       in_states[trackid]->TopMatrix(m);              // could benefit from interal vec
0287       auto tmp = m.Transform(globalpoints[trackid]); // could benefit from internal vec
0288
0289       using vecCore::AssignLane;
0290       AssignLane(localpoint.x(), i, tmp.x());
0291       AssignLane(localpoint.y(), i, tmp.y());
0292       AssignLane(localpoint.z(), i, tmp.z());
0293       tmp = m.TransformDirection(globaldirs[trackid]); // could benefit from internal vec
0294       AssignLane(localdir.x(), i, tmp.x());
0295       AssignLane(localdir.y(), i, tmp.y());
0296       AssignLane(localdir.z(), i, tmp.z());
0297     }
0298   }
0299
0300   // version used for SIMD processing
0301   // should be specialized in Impl for faster treatment
0302   template <typename T, unsigned int ChunkSize>
0303   VECGEOM_FORCE_INLINE
0304   VECCORE_ATT_HOST_DEVICE
0305   static void DoGlobalToLocalTransformations(NavStatePool const &in_states, SOA3D<Precision> const &globalpoints,
0306                                              SOA3D<Precision> const &globaldirs, unsigned int from_index,
0307                                              Vector3D<T> &localpoint, Vector3D<T> &localdir)
0308   {
0309     for (unsigned int i = 0; i < ChunkSize; ++i) {
0310       unsigned int trackid = from_index + i;
0311       Transformation3D m;
0312       //  assert(in_states[trackid]->Top()->GetLogicalVolume() == lvol &&
0313       //         "not all states in same logical volume"); // the logical volume of all the states should be the same
0314       in_states[trackid]->TopMatrix(m);              // could benefit from internal vec
0315       auto tmp = m.Transform(globalpoints[trackid]); // could benefit from internal vec
0316
0317       using vecCore::AssignLane;
0318       AssignLane(localpoint.x(), i, tmp.x());
0319       AssignLane(localpoint.y(), i, tmp.y());
0320       AssignLane(localpoint.z(), i, tmp.z());
0321       tmp = m.TransformDirection(globaldirs[trackid]); // could benefit from internal vec
0322       AssignLane(localdir.x(), i, tmp.x());
0323       AssignLane(localdir.y(), i, tmp.y());
0324       AssignLane(localdir.z(), i, tmp.z());
0325     }
0326   }
0327 };
0328
0329 //! template class providing a standard implementation for
0330 //! some interfaces in VNavigator (using the CRT pattern)
0331 template <typename Impl, bool MotherIsConvex = false>
0332 class VNavigatorHelper : public VNavigator {
0333 protected:
0334   using VNavigator::VNavigator;
0335
0336 public:
0337   // the default implementation for hit detection with daughters for a chunk of data
0338   // is to loop over the implementation for the scalar case
0339   // this static function may be overridden by the specialized implementations (such as done in NewSimpleNavigator)
0340   // the from_index, to_index indicate which states from the NavigationState ** are actually treated
0341   // in the worst case, we might have to implement this stuff over there
0342   template <typename T, unsigned int ChunkSize>
0343   VECCORE_ATT_HOST_DEVICE
0344   static void DaughterIntersectionsLooper(VNavigator const *nav, LogicalVolume const *lvol,
0345                                           Vector3D<T> const &localpoint, Vector3D<T> const &localdir,
0346                                           NavigationState const *const *in_states, NavigationState **out_states,
0347                                           unsigned int from_index, Precision *out_steps,
0348                                           VPlacedVolume const *hitcandidates[ChunkSize])
0349   {
0350     // dispatch to ordinary implementation ( which itself might be vectorized )
0351     using vecCore::LaneAt;
0352     for (unsigned int i = 0; i < ChunkSize; ++i) {
0353       unsigned int trackid = from_index + i;
0354       ((Impl *)nav)
0355           ->Impl::CheckDaughterIntersections(
0356               lvol,
0357               Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i), LaneAt(localpoint.z(), i)),
0358               Vector3D<Precision>(LaneAt(localdir.x(), i), LaneAt(localdir.y(), i), LaneAt(localdir.z(), i)),
0359               in_states[trackid], out_states[trackid], out_steps[trackid], hitcandidates[i]);
0360     }
0361   }
0362
0363   // the default implementation for hit detection with daughters for a chunk of data
0364   // is to loop over the implementation for the scalar case
0365   // this static function may be overridden by the specialized implementations (such as done in NewSimpleNavigator)
0366   // the from_index, to_index indicate which states from the NavigationState ** are actually treated
0367   // in the worst case, we might have to implement this stuff over there
0368   template <typename T, unsigned int ChunkSize>
0369   VECCORE_ATT_HOST_DEVICE
0370   static void DaughterIntersectionsLooper(VNavigator const *nav, LogicalVolume const *lvol,
0371                                           Vector3D<T> const &localpoint, Vector3D<T> const &localdir,
0372                                           NavStatePool const &in_states, NavStatePool &out_states,
0373                                           unsigned int from_index, Precision *out_steps,
0374                                           VPlacedVolume const *hitcandidates[ChunkSize])
0375   {
0376     // dispatch to ordinary implementation ( which itself might be vectorized )
0377     using vecCore::LaneAt;
0378     for (unsigned int i = 0; i < ChunkSize; ++i) {
0379       unsigned int trackid = from_index + i;
0380       ((Impl *)nav)
0381           ->Impl::CheckDaughterIntersections(
0382               lvol,
0383               Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i), LaneAt(localpoint.z(), i)),
0384               Vector3D<Precision>(LaneAt(localdir.x(), i), LaneAt(localdir.y(), i), LaneAt(localdir.z(), i)),
0385               in_states[trackid], out_states[trackid], out_steps[trackid], hitcandidates[i]);
0386     }
0387   }
0388
0389   // the default implementation for hit detection with daughters for a chunk of data
0390   // is to loop over the implementation for the scalar case
0391   // similar to previous version; does not care about out_state
0392   template <typename T, unsigned int ChunkSize>
0393   VECCORE_ATT_HOST_DEVICE
0394   static void DaughterIntersectionsLooper(VNavigator const *nav, LogicalVolume const *lvol,
0395                                           Vector3D<T> const &localpoint, Vector3D<T> const &localdir,
0396                                           NavigationState const *const *in_states, unsigned int from_index,
0397                                           Precision *out_steps, VPlacedVolume const *hitcandidates[ChunkSize])
0398   {
0399     // dispatch to ordinary implementation ( which itself might be vectorized )
0400     using vecCore::LaneAt;
0401     for (unsigned int i = 0; i < ChunkSize; ++i) {
0402       unsigned int trackid = from_index + i;
0403       ((Impl *)nav)
0404           ->Impl::CheckDaughterIntersections(
0405               lvol,
0406               Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i), LaneAt(localpoint.z(), i)),
0407               Vector3D<Precision>(LaneAt(localdir.x(), i), LaneAt(localdir.y(), i), LaneAt(localdir.z(), i)),
0408               in_states[trackid], nullptr, out_steps[trackid], hitcandidates[i]);
0409     }
0410   }
0411
0412   // the default implementation for safety calculation for a chunk of data
0413   template <typename T, unsigned int ChunkSize> // we may go to Backend as template parameter in future
0414   VECCORE_ATT_HOST_DEVICE static void SafetyLooper(VNavigator const *nav, VPlacedVolume const *pvol,
0415                                                    Vector3D<T> const &localpoint, unsigned int from_index,
0416                                                    bool const *calcsafeties, Precision *out_safeties)
0417   {
0418     // dispatch to ordinary implementation ( which itself might be vectorized )
0419     // TODO: check calcsafeties if we need to do this at all
0420     using SafetyE_t = typename Impl::SafetyEstimator_t;
0421     using Bool_v    = vecCore::Mask_v<T>;
0422     Bool_v m;
0423     for (unsigned int i = 0; i < ChunkSize; ++i) {
0424       vecCore::AssignMaskLane(m, i, calcsafeties[from_index + i]);
0425     }
0426
0427     T safety(0.);
0428     if (!vecCore::MaskEmpty(m)) {
0429       safety = ((SafetyE_t *)nav->GetSafetyEstimator())->ComputeSafetyForLocalPoint(localpoint, pvol, m);
0430     }
0431     vecCore::Store(safety, &out_safeties[from_index]);
0432   }
0433
0434 public:
0435   VECCORE_ATT_HOST_DEVICE
0436   virtual Precision ComputeStepAndPropagatedState(Vector3D<Precision> const &globalpoint,
0437                                                   Vector3D<Precision> const &globaldir, Precision step_limit,
0438                                                   NavigationState const &in_state,
0439                                                   NavigationState &out_state) const override
0440   {
0441 #ifdef DEBUGNAV
0442     static size_t counter = 0;
0443     counter++;
0444 #endif
0445     // calculate local point/dir from global point/dir
0446     // call the static function for this provided/specialized by the Impl
0447     Vector3D<Precision> localpoint;
0448     Vector3D<Precision> localdir;
0449     Impl::DoGlobalToLocalTransformation(in_state, globalpoint, globaldir, localpoint, localdir);
0450
0451     VPlacedVolume const *hitcandidate = nullptr;
0452     auto pvol                         = in_state.Top();
0453     auto lvol                         = pvol->GetLogicalVolume();
0454     Precision step                    = Impl::TreatDistanceToMother(pvol, localpoint, localdir, step_limit);
0455     // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0456     if (lvol->GetDaughters().size() > 0)
0457       ((Impl *)this)
0458           ->Impl::CheckDaughterIntersections(lvol, localpoint, localdir, &in_state, &out_state, step, hitcandidate);
0459
0460     // fix state
0461     bool done;
0462     step = Impl::PrepareOutState(in_state, out_state, step, step_limit, hitcandidate, done);
0463     if (done) {
0464       if (out_state.Top() != nullptr) {
0465         assert(!out_state.Top()->GetLogicalVolume()->GetUnplacedVolume()->IsAssembly());
0466       }
0467       return step;
0468     }
0469     // step was physics limited
0470     if (!out_state.IsOnBoundary()) return step;
0471
0472     // otherwise if necessary do a relocation
0473     // try relocation to refine out_state to correct location after the boundary
0474     ((Impl *)this)->Impl::Relocate(MovePointAfterBoundary(localpoint, localdir, step), in_state, out_state);
0475     if (out_state.Top() != nullptr) {
0476       while (out_state.Top()->IsAssembly()) {
0477         out_state.Pop();
0478       }
0479       assert(!out_state.Top()->GetLogicalVolume()->GetUnplacedVolume()->IsAssembly());
0480     }
0481     return step;
0482   }
0483
0484   virtual Precision ComputeStep(Vector3D<Precision> const &globalpoint, Vector3D<Precision> const &globaldir,
0485                                 Precision step_limit, NavigationState const &in_state,
0486                                 NavigationState &out_state) const override
0487   {
0488 #ifdef DEBUGNAV
0489     static size_t counter = 0;
0490     counter++;
0491 #endif
0492     // calculate local point/dir from global point/dir
0493     // call the static function for this provided/specialized by the Impl
0494     Vector3D<Precision> localpoint;
0495     Vector3D<Precision> localdir;
0496     Impl::DoGlobalToLocalTransformation(in_state, globalpoint, globaldir, localpoint, localdir);
0497
0498     VPlacedVolume const *hitcandidate = nullptr;
0499     auto pvol                         = in_state.Top();
0500     auto lvol                         = pvol->GetLogicalVolume();
0501     Precision step                    = Impl::TreatDistanceToMother(pvol, localpoint, localdir, step_limit);
0502     // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0503     if (lvol->GetDaughters().size() > 0)
0504       ((Impl *)this)
0505           ->Impl::CheckDaughterIntersections(lvol, localpoint, localdir, &in_state, &out_state, step, hitcandidate);
0506
0507     // fix state
0508     bool done;
0509     step = Impl::PrepareOutState(in_state, out_state, step, step_limit, hitcandidate, done);
0510     if (done) {
0511       if (out_state.Top() != nullptr) {
0512         assert(!out_state.Top()->GetLogicalVolume()->GetUnplacedVolume()->IsAssembly());
0513       }
0514       return step;
0515     }
0516     // step was physics limited
0517     if (!out_state.IsOnBoundary()) return step;
0518
0519     return step;
0520   }
0521
0522   VECCORE_ATT_HOST_DEVICE
0523   virtual Precision ComputeStepAndSafety(Vector3D<Precision> const &globalpoint, Vector3D<Precision> const &globaldir,
0524                                          Precision step_limit, NavigationState &in_state, bool calcsafety,
0525                                          Precision &safety, bool indicateDaughterHit = false) const override
0526   {
0527 // FIXME: combine this kernel and the one for ComputeStep() into one generic function
0528 #ifdef DEBUGNAV
0529     static size_t counter = 0;
0530     counter++;
0531 #endif
0532     // calculate local point/dir from global point/dir
0533     // call the static function for this provided/specialized by the Impl
0534     Vector3D<Precision> localpoint;
0535     Vector3D<Precision> localdir;
0536     NavigationState *out_state = nullptr;
0537
0538     Impl::DoGlobalToLocalTransformation(in_state, globalpoint, globaldir, localpoint, localdir);
0539
0540     // get safety first ( the benefit here is that we reuse the local points )
0541     using SafetyE_t = typename Impl::SafetyEstimator_t;
0542     if (calcsafety) {
0543       // call the appropriate safety Estimator
0544       safety = ((SafetyE_t *)fSafetyEstimator)->SafetyE_t::ComputeSafetyForLocalPoint(localpoint, in_state.Top());
0545     }
0546
0547     VPlacedVolume const *hitcandidate = nullptr;
0548     auto pvol                         = in_state.Top();
0549     auto lvol                         = pvol->GetLogicalVolume();
0550     Precision step                    = step_limit;
0551
0552     // is the next object certainly further away than the safety
0553     bool safetydone = calcsafety && safety >= step;
0554
0555     if (!safetydone) {
0556       step = Impl::TreatDistanceToMother(pvol, localpoint, localdir, step_limit);
0557       // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0558       if (lvol->GetDaughters().size() > 0)
0559         ((Impl *)this)
0560             ->Impl::CheckDaughterIntersections(lvol, localpoint, localdir, &in_state, out_state, step, hitcandidate);
0561     }
0562     if (indicateDaughterHit && hitcandidate) in_state.Push(hitcandidate);
0563     return Min(step, step_limit);
0564   }
0565
0566   VECCORE_ATT_HOST_DEVICE
0567   virtual void ComputeStepsAndSafeties(SOA3D<Precision> const &globalpoints, SOA3D<Precision> const &globaldirs,
0568                                        Precision const *step_limits, NavigationState const *const *in_states,
0569                                        Precision *out_steps, bool const *calcsafeties,
0570                                        Precision *out_safeties) const override
0571   {
0572     // process SIMD part and TAIL part
0573     using Real_v = vecgeom::VectorBackend::Real_v;
0574     (void)(Real_v *)(nullptr); // Avoid spurrious warning: unused type alias 'Real_v' [-Wunused-local-typedef]
0575     const auto size = globalpoints.size();
0576     auto pvol       = in_states[0]->Top();
0577     auto lvol       = pvol->GetLogicalVolume();
0578     // loop over all tracks in chunks
0579     auto i             = decltype(size){0};
0580     constexpr auto kVS = vecCore::VectorSize<Real_v>();
0581     for (; i < (size - (kVS - 1)); i += kVS) {
0582       NavigateAChunkNoReloc<Real_v, kVS>(this, pvol, lvol, globalpoints, globaldirs, step_limits, in_states, out_steps,
0583                                          calcsafeties, out_safeties, i);
0584     }
0585
0586     // fall back to scalar interface for tail treatment
0587     for (; i < size; ++i) {
0588       out_steps[i] = ((Impl *)this)
0589                          ->Impl::ComputeStepAndSafety(globalpoints[i], globaldirs[i], step_limits[i],
0590                                                       *const_cast<NavigationState *>(in_states[i]), calcsafeties[i],
0591                                                       out_safeties[i]);
0592     }
0593     i = decltype(size){0};
0594     for (; i < size; ++i)
0595       out_steps[i] = vecCore::math::Min(out_steps[i], step_limits[i]);
0596   }
0597
0598   // this kernel is a generic implementation to navigate with chunks of data
0599   // can be used also for the scalar implementation
0600   template <typename T, unsigned int ChunkSize>
0601   static void NavigateAChunk(VNavigator const *__restrict__ nav, VPlacedVolume const *__restrict__ pvol,
0602                              LogicalVolume const *__restrict__ lvol, SOA3D<Precision> const &__restrict__ globalpoints,
0603                              SOA3D<Precision> const &__restrict__ globaldirs, Precision const *__restrict__ step_limits,
0604                              NavigationState const *const *__restrict__ in_states,
0605                              NavigationState **__restrict__ out_states, Precision *__restrict__ out_steps,
0606                              unsigned int from_index)
0607   {
0608
0609     VPlacedVolume const *hitcandidates[ChunkSize] = {}; // initialize all to nullptr
0610
0611     Vector3D<T> localpoint, localdir;
0612     Impl::template DoGlobalToLocalTransformations<T, ChunkSize>(in_states, globalpoints, globaldirs, from_index,
0613                                                                 localpoint, localdir);
0614
0615     T slimit(vecCore::FromPtr<T>(step_limits + from_index));
0616     // need to calc DistanceToOut first
0617     T step = Impl::template TreatDistanceToMother<T>(pvol, localpoint, localdir, slimit);
0618     vecCore::Store(step, out_steps + from_index);
0619
0620     // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0621     Impl::template DaughterIntersectionsLooper<T, ChunkSize>(nav, lvol, localpoint, localdir, in_states, out_states,
0622                                                               from_index, out_steps, hitcandidates);
0623
0624     // fix state ( seems to be serial so we iterate over indices )
0625     for (unsigned int i = 0; i < ChunkSize; ++i) {
0626       unsigned int trackid = from_index + i;
0627       bool done;
0628       out_steps[trackid] = Impl::PrepareOutState(*in_states[trackid], *out_states[trackid], out_steps[trackid],
0629                                                  vecCore::LaneAt(slimit, i), hitcandidates[i], done);
0630       if (done) continue;
0631       // step was physics limited
0632       if (!out_states[trackid]->IsOnBoundary()) continue;
0633       // otherwise if necessary do a relocation
0634       // try relocation to refine out_state to correct location after the boundary
0635       using vecCore::LaneAt;
0636       ((Impl *)nav)
0637           ->Impl::Relocate(
0638               MovePointAfterBoundary(
0639                   Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i), LaneAt(localpoint.z(), i)),
0640                   Vector3D<Precision>(LaneAt(localdir.x(), i), LaneAt(localdir.y(), i), LaneAt(localdir.z(), i)),
0641                   out_steps[trackid]),
0642               *in_states[trackid], *out_states[trackid]);
0643     }
0644   }
0645
0646   // this kernel is a generic implementation to navigate with chunks of data
0647   // can be used also for the scalar imple
0648   // same as above but including safety treatment
0649   // TODO: remerge with above to avoid code duplication
0650   template <typename T, unsigned int ChunkSize>
0651   VECCORE_ATT_HOST_DEVICE
0652   static void NavigateAChunk(VNavigator const *__restrict__ nav, VPlacedVolume const *__restrict__ pvol,
0653                              LogicalVolume const *__restrict__ lvol, SOA3D<Precision> const &__restrict__ globalpoints,
0654                              SOA3D<Precision> const &__restrict__ globaldirs, Precision const *__restrict__ step_limits,
0655                              NavigationState const *const *__restrict__ in_states,
0656                              NavigationState **__restrict__ out_states, Precision *__restrict__ out_steps,
0657                              bool const *__restrict__ calcsafeties, Precision *__restrict__ out_safeties,
0658                              unsigned int from_index)
0659   {
0660
0661     VPlacedVolume const *hitcandidates[ChunkSize] = {}; // initialize all to nullptr
0662
0663     Vector3D<T> localpoint, localdir;
0664     Impl::template DoGlobalToLocalTransformations<T, ChunkSize>(in_states, globalpoints, globaldirs, from_index,
0665                                                                 localpoint, localdir);
0666
0667     // safety part
0668     Impl::template SafetyLooper<T, ChunkSize>(nav, pvol, localpoint, from_index, calcsafeties, out_safeties);
0669
0670     T slimit(vecCore::FromPtr<T>(step_limits + from_index)); // will only work with new ScalarWrapper
0671     // need to calc DistanceToOut first
0672     T step = Impl::template TreatDistanceToMother<T>(pvol, localpoint, localdir, slimit);
0673     vecCore::Store(step, out_steps + from_index);
0674
0675     // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0676     Impl::template DaughterIntersectionsLooper<T, ChunkSize>(nav, lvol, localpoint, localdir, in_states, out_states,
0677                                                               from_index, out_steps, hitcandidates);
0678
0679     using vecCore::LaneAt;
0680
0681     // fix state ( seems to be serial so we iterate over indices )
0682     for (unsigned int i = 0; i < ChunkSize; ++i) {
0683       unsigned int trackid = from_index + i;
0684       bool done;
0685       out_steps[trackid] = Impl::PrepareOutState(*in_states[trackid], *out_states[trackid], out_steps[trackid],
0686                                                  LaneAt(slimit, i), hitcandidates[i], done);
0687       if (done) continue;
0688       // step was physics limited
0689       if (!out_states[trackid]->IsOnBoundary()) continue;
0690       // otherwise if necessary do a relocation
0691       // try relocation to refine out_state to correct location after the boundary
0692       ((Impl *)nav)
0693           ->Impl::Relocate(
0694               MovePointAfterBoundary(
0695                   Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i), LaneAt(localpoint.z(), i)),
0696                   Vector3D<Precision>(LaneAt(localdir.x(), i), LaneAt(localdir.y(), i), LaneAt(localdir.z(), i)),
0697                   out_steps[trackid]),
0698               *in_states[trackid], *out_states[trackid]);
0699     }
0700   }
0701
0702   // this kernel is a generic implementation to navigate with chunks of data
0703   // can be used also for the scalar implementation
0704   // same as above but including safety treatment
0705   // TODO: remerge with above to avoid code duplication
0706   template <typename T, unsigned int ChunkSize>
0707   VECCORE_ATT_HOST_DEVICE
0708   static void NavigateAChunk(VNavigator const *__restrict__ nav, VPlacedVolume const *__restrict__ pvol,
0709                              LogicalVolume const *__restrict__ lvol, SOA3D<Precision> const &__restrict__ globalpoints,
0710                              SOA3D<Precision> const &__restrict__ globaldirs, Precision const *__restrict__ step_limits,
0711                              NavStatePool const &in_states, NavStatePool &out_states, Precision *__restrict__ out_steps,
0712                              bool const *__restrict__ calcsafeties, Precision *__restrict__ out_safeties,
0713                              unsigned int from_index)
0714   {
0715     VPlacedVolume const *hitcandidates[ChunkSize] = {}; // initialize all to nullptr
0716
0717     Vector3D<T> localpoint, localdir;
0718     Impl::template DoGlobalToLocalTransformations<T, ChunkSize>(in_states, globalpoints, globaldirs, from_index,
0719                                                                 localpoint, localdir);
0720
0721     // safety part
0722     Impl::template SafetyLooper<T, ChunkSize>(nav, pvol, localpoint, from_index, calcsafeties, out_safeties);
0723
0724     T slimit(vecCore::FromPtr<T>(step_limits + from_index)); // will only work with new ScalarWrapper
0725     // need to calc DistanceToOut first
0726     T step = Impl::template TreatDistanceToMother<T>(pvol, localpoint, localdir, slimit);
0727     vecCore::Store(step, out_steps + from_index);
0728
0729     // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0730     Impl::template DaughterIntersectionsLooper<T, ChunkSize>(nav, lvol, localpoint, localdir, in_states, out_states,
0731                                                               from_index, out_steps, hitcandidates);
0732
0733     using vecCore::LaneAt;
0734
0735     // fix state ( seems to be serial so we iterate over indices )
0736     for (unsigned int i = 0; i < ChunkSize; ++i) {
0737       unsigned int trackid = from_index + i;
0738       bool done;
0739       out_steps[trackid] = Impl::PrepareOutState(*in_states[trackid], *out_states[trackid], out_steps[trackid],
0740                                                  LaneAt(slimit, i), hitcandidates[i], done);
0741       if (done) continue;
0742       // step was physics limited
0743       if (!out_states[trackid]->IsOnBoundary()) continue;
0744       // otherwise if necessary do a relocation
0745       // try relocation to refine out_state to correct location after the boundary
0746       ((Impl *)nav)
0747           ->Impl::Relocate(
0748               MovePointAfterBoundary(
0749                   Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i), LaneAt(localpoint.z(), i)),
0750                   Vector3D<Precision>(LaneAt(localdir.x(), i), LaneAt(localdir.y(), i), LaneAt(localdir.z(), i)),
0751                   out_steps[trackid]),
0752               *in_states[trackid], *out_states[trackid]);
0753     }
0754   }
0755
0756   // this kernel is a generic implementation to navigate with chunks of data
0757   // can be used also for the scalar implementation
0758   // same as above but including safety treatment
0759   // TODO: remerge with above to avoid code duplication
0760   template <typename T, unsigned int ChunkSize>
0761   VECCORE_ATT_HOST_DEVICE
0762   static void NavigateAChunkNoReloc(VNavigator const *__restrict__ nav, VPlacedVolume const *__restrict__ pvol,
0763                                     LogicalVolume const *__restrict__ lvol,
0764                                     SOA3D<Precision> const &__restrict__ globalpoints,
0765                                     SOA3D<Precision> const &__restrict__ globaldirs,
0766                                     Precision const *__restrict__ step_limits,
0767                                     NavigationState const *const *__restrict__ in_states,
0768                                     Precision *__restrict__ out_steps, bool const *__restrict__ calcsafeties,
0769                                     Precision *__restrict__ out_safeties, unsigned int from_index)
0770   {
0771     VPlacedVolume const *hitcandidates[ChunkSize] = {}; // initialize all to nullptr
0772
0773     Vector3D<T> localpoint, localdir;
0774     Impl::template DoGlobalToLocalTransformations<T, ChunkSize>(in_states, globalpoints, globaldirs, from_index,
0775                                                                 localpoint, localdir);
0776
0777     // safety part
0778     Impl::template SafetyLooper<T, ChunkSize>(nav, pvol, localpoint, from_index, calcsafeties, out_safeties);
0779
0780     T slimit(vecCore::FromPtr<T>(step_limits + from_index));
0781     // need to calc DistanceToOut first
0782     T step = Impl::template TreatDistanceToMother<T>(pvol, localpoint, localdir, slimit);
0783     vecCore::Store(step, out_steps + from_index);
0784
0785     // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0786     Impl::template DaughterIntersectionsLooper<T, ChunkSize>(nav, lvol, localpoint, localdir, in_states, from_index,
0787                                                               out_steps, hitcandidates);
0788   }
0789
0790   // generic implementation for the vector interface
0791   // this implementation tries to process everything in vector CHUNKS
0792   // at the very least this enables at least the DistanceToOut call to be vectorized
0793   virtual void ComputeStepsAndPropagatedStates(SOA3D<Precision> const &__restrict__ globalpoints,
0794                                                SOA3D<Precision> const &__restrict__ globaldirs,
0795                                                Precision const *__restrict__ step_limit,
0796                                                NavigationState const *const *__restrict__ in_states,
0797                                                NavigationState **__restrict__ out_states,
0798                                                Precision *__restrict__ out_steps) const override
0799   {
0800
0801     // process SIMD part and TAIL part
0802     // something like
0803     using Real_v = vecgeom::VectorBackend::Real_v;
0804     (void)(Real_v *)(nullptr); // Avoid spurrious warning: unused type alias 'Real_v' [-Wunused-local-typedef]
0805     const auto size = globalpoints.size();
0806     auto pvol       = in_states[0]->Top();
0807     auto lvol       = pvol->GetLogicalVolume();
0808     // loop over all tracks in chunks
0809     auto i             = decltype(size){0};
0810     constexpr auto kVS = vecCore::VectorSize<Real_v>();
0811     for (; i < (size - (kVS - 1)); i += kVS) {
0812       NavigateAChunk<Real_v, kVS>(this, pvol, lvol, globalpoints, globaldirs, step_limit, in_states, out_states,
0813                                   out_steps, i);
0814     }
0815
0816     // fall back to scalar interface for tail treatment
0817     for (; i < size; ++i) {
0818       out_steps[i] = ((Impl *)this)
0819                          ->Impl::ComputeStepAndPropagatedState(globalpoints[i], globaldirs[i], step_limit[i],
0820                                                                *in_states[i], *out_states[i]);
0821     }
0822   }
0823
0824   // generic implementation for the vector interface -- using NavStatePool instead of NavigationState**
0825   // this implementation tries to process everything in vector CHUNKS
0826   // at the very least this enables at least the DistanceToOut call to be vectorized
0827   VECCORE_ATT_HOST_DEVICE
0828   virtual void ComputeStepsAndSafetiesAndPropagatedStates(SOA3D<Precision> const &__restrict__ globalpoints,
0829                                                           SOA3D<Precision> const &__restrict__ globaldirs,
0830                                                           Precision const *__restrict__ step_limit,
0831                                                           NavStatePool const &in_states, NavStatePool &out_states,
0832                                                           Precision *__restrict__ out_steps,
0833                                                           bool const *__restrict__ calcsafeties,
0834                                                           Precision *__restrict__ out_safeties) const override
0835   {
0836     // process SIMD part and TAIL part
0837     // something like
0838     const auto size = globalpoints.size();
0839     auto pvol       = in_states[0]->Top();
0840     auto lvol       = pvol->GetLogicalVolume();
0841     // loop over all tracks in chunks
0842     auto i             = decltype(size){0};
0843     constexpr auto kVS = vecCore::VectorSize<Real_v>();
0844     for (; i < (size - (kVS - 1)); i += kVS) {
0845       NavigateAChunk<Real_v, kVS>(this, pvol, lvol, globalpoints, globaldirs, step_limit, in_states, out_states,
0846                                   out_steps, calcsafeties, out_safeties, i);
0847     }
0848
0849     // fall back to scalar interface for tail treatment
0850     for (; i < size; ++i) {
0851       out_steps[i] = ((Impl *)this)
0852                          ->Impl::ComputeStepAndSafetyAndPropagatedState(globalpoints[i], globaldirs[i], step_limit[i],
0853                                                                         *in_states[i], *out_states[i], calcsafeties[i],
0854                                                                         out_safeties[i]);
0855     }
0856   }
0857
0858   // generic implementation for the vector interface
0859   // this implementation tries to process everything in vector CHUNKS
0860   // at the very least this enables at least the DistanceToOut call to be vectorized
0861   VECCORE_ATT_HOST_DEVICE
0862   virtual void ComputeStepsAndSafetiesAndPropagatedStates(
0863       SOA3D<Precision> const &__restrict__ globalpoints, SOA3D<Precision> const &__restrict__ globaldirs,
0864       Precision const *__restrict__ step_limit, NavigationState const *const *__restrict__ in_states,
0865       NavigationState **__restrict__ out_states, Precision *__restrict__ out_steps,
0866       bool const *__restrict__ calcsafeties, Precision *__restrict__ out_safeties) const override
0867   {
0868
0869     // process SIMD part and TAIL part
0870     // something like
0871     using Real_v = vecgeom::VectorBackend::Real_v;
0872     (void)(Real_v *)(nullptr); // Avoid spurrious warning: unused type alias 'Real_v' [-Wunused-local-typedef]
0873     const auto size = globalpoints.size();
0874     auto pvol       = in_states[0]->Top();
0875     auto lvol       = pvol->GetLogicalVolume();
0876     // loop over all tracks in chunks
0877     auto i             = decltype(size){0};
0878     constexpr auto kVS = vecCore::VectorSize<Real_v>();
0879     for (; i < (size - (kVS - 1)); i += kVS) {
0880       NavigateAChunk<Real_v, kVS>(this, pvol, lvol, globalpoints, globaldirs, step_limit, in_states, out_states,
0881                                   out_steps, calcsafeties, out_safeties, i);
0882     }
0883     // fall back to scalar interface for tail treatment
0884     for (; i < size; ++i) {
0885       out_steps[i] = ((Impl *)this)
0886                          ->Impl::ComputeStepAndSafetyAndPropagatedState(globalpoints[i], globaldirs[i], step_limit[i],
0887                                                                         *in_states[i], *out_states[i], calcsafeties[i],
0888                                                                         out_safeties[i]);
0889     }
0890   }
0891
0892 #ifdef TRIVIAL
0893   // another generic implementation for the vector interface
0894   // this implementation just loops over the scalar interface
0895   virtual void ComputeStepsAndPropagatedStates(SOA3D<Precision> const &__restrict__ globalpoints,
0896                                                SOA3D<Precision> const &__restrict__ globaldirs,
0897                                                Precision const *__restrict__ step_limit,
0898                                                NavigationState const *const *__restrict__ in_states,
0899                                                NavigationState **__restrict__ out_states,
0900                                                Precision *__restrict__ out_steps) const override
0901   {
0902     for (unsigned int i = 0; i < globalpoints.size(); ++i) {
0903       out_steps[i] = ((Impl *)this)
0904                          ->Impl::ComputeStepAndPropagatedState(globalpoints[i], globaldirs[i], step_limit[i],
0905                                                                *in_states[i], *out_states[i]);
0906     }
0907   }
0908
0909   VECCORE_ATT_HOST_DEVICE
0910   virtual void ComputeStepsAndSafetiesAndPropagatedStates(
0911       SOA3D<Precision> const &__restrict__ globalpoints, SOA3D<Precision> const &__restrict__ globaldirs,
0912       Precision const *__restrict__ step_limit, NavigationState const *const *__restrict__ in_states,
0913       NavigationState **__restrict__ out_states, Precision *__restrict__ out_steps,
0914       bool const *__restrict__ calcsafeties, Precision *__restrict__ out_safeties) const override
0915   {
0916
0917     for (unsigned int i = 0; i < globalpoints.size(); ++i) {
0918       out_steps[i] = ((Impl *)this)
0919                          ->Impl::ComputeStepAndSafetyAndPropagatedState(globalpoints[i], globaldirs[i], step_limit[i],
0920                                                                         *in_states[i], *out_states[i], calcsafeties[i],
0921                                                                         out_safeties[i]);
0922     }
0923   }
0924 #endif
0925
0926   // a similar interface also returning the safety
0927   // TODO: reduce this evident code duplication with ComputeStepAndPropagatedState
0928   VECCORE_ATT_HOST_DEVICE
0929   virtual Precision ComputeStepAndSafetyAndPropagatedState(Vector3D<Precision> const &globalpoint,
0930                                                            Vector3D<Precision> const &globaldir, Precision step_limit,
0931                                                            NavigationState const &__restrict__ in_state,
0932                                                            NavigationState &__restrict__ out_state, bool calcsafety,
0933                                                            Precision &safety_out) const override
0934   {
0935     // calculate local point/dir from global point/dir
0936     Vector3D<Precision> localpoint;
0937     Vector3D<Precision> localdir;
0938     Impl::DoGlobalToLocalTransformation(in_state, globalpoint, globaldir, localpoint, localdir);
0939
0940     // get safety first ( the benefit here is that we reuse the local points )
0941     using SafetyE_t = typename Impl::SafetyEstimator_t;
0942     safety_out      = 0.;
0943     if (calcsafety) {
0944       // call the appropriate safety Estimator
0945       safety_out = ((SafetyE_t *)fSafetyEstimator)->SafetyE_t::ComputeSafetyForLocalPoint(localpoint, in_state.Top());
0946     }
0947
0948     VPlacedVolume const *hitcandidate = nullptr;
0949     auto pvol                         = in_state.Top();
0950     auto lvol                         = pvol->GetLogicalVolume();
0951     Precision step                    = Impl::TreatDistanceToMother(pvol, localpoint, localdir, step_limit);;
0952     if (lvol->GetDaughters().size() > 0)
0953       // "suck in" algorithm from Impl and treat hit detection in local coordinates for daughters
0954       ((Impl *)this)
0955           ->Impl::CheckDaughterIntersections(lvol, localpoint, localdir, &in_state, &out_state, step, hitcandidate);
0956
0957     // fix state
0958     bool done;
0959     step = Impl::PrepareOutState(in_state, out_state, step, step_limit, hitcandidate, done);
0960     if (done) return step;
0961
0962     // step was physics limited
0963     if (!out_state.IsOnBoundary()) return step;
0964
0965     // otherwise if necessary do a relocation
0966     // try relocation to refine out_state to correct location after the boundary
0967     ((Impl *)this)->Impl::Relocate(MovePointAfterBoundary(localpoint, localdir, step), in_state, out_state);
0968     return step;
0969   }
0970
0971 protected:
0972   // a common relocate method ( to calculate propagated states after the boundary )
0973   VECGEOM_FORCE_INLINE
0974   VECCORE_ATT_HOST_DEVICE
0975   virtual void Relocate(Vector3D<Precision> const &pointafterboundary, NavigationState const &__restrict__ in_state,
0976                         NavigationState &__restrict__ out_state) const override
0977   {
0978     // this means that we are leaving the mother
0979     // alternatively we could use nextvolumeindex like before
0980     if (out_state.Top() == in_state.Top()) {
0981       GlobalLocator::RelocatePointFromPathForceDifferent(pointafterboundary, out_state);
0982 #ifdef CHECK_RELOCATION_ERRORS
0983       assert(in_state.Distance(out_state) != 0 && " error relocating when leaving ");
0984 #endif
0985     } else {
0986       // continue directly further down ( next volume should have been stored in out_state already )
0987       VPlacedVolume const *nextvol = out_state.Top();
0988       out_state.Pop();
0989       GlobalLocator::LocateGlobalPoint(nextvol, nextvol->GetTransformation()->Transform(pointafterboundary), out_state,
0990                                        false);
0991 #ifdef CHECK_RELOCATION_ERRORS
0992       assert(in_state.Distance(out_state) != 0 && " error relocating when entering ");
0993 #endif
0994       return;
0995     }
0996   }
0997
0998 public:
0999   static const char *GetClassName() { return Impl::gClassNameString; }
1000
1001   virtual const char *GetName() const override { return GetClassName(); }
1002 }; // end class VNavigatorHelper
1003 } // namespace VECGEOM_IMPL_NAMESPACE
1004 } // namespace vecgeom
1005
1006 #endif /* NAVIGATION_VNAVIGATOR_H_ */