VecGeom/navigation/NewSimpleNavigator.h

0001 /*
0002  * NewSimpleNavigator.h
0003  *
0004  *  Created on: 17.09.2015
0005  *      Author: swenzel
0006  */
0007
0008 #ifndef NAVIGATION_NEWSIMPLENAVIGATOR_H_
0009 #define NAVIGATION_NEWSIMPLENAVIGATOR_H_
0010
0011 #include "VNavigator.h"
0012 #include "SimpleSafetyEstimator.h"
0013
0014 namespace vecgeom {
0015 inline namespace VECGEOM_IMPL_NAMESPACE {
0016
0017 // A very basic implementation of a navigator ( brute force which scales linearly with the number of daughters )
0018 template <bool MotherIsConvex = false>
0019 class NewSimpleNavigator : public VNavigatorHelper<class NewSimpleNavigator<MotherIsConvex>, MotherIsConvex> {
0020
0021 private:
0022   VECCORE_ATT_DEVICE
0023   NewSimpleNavigator()
0024       : VNavigatorHelper<class NewSimpleNavigator<MotherIsConvex>, MotherIsConvex>(SimpleSafetyEstimator::Instance()) {
0025   } VECCORE_ATT_DEVICE virtual ~NewSimpleNavigator() {}
0026
0027 public:
0028   VECGEOM_FORCE_INLINE
0029   VECCORE_ATT_HOST_DEVICE
0030   virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0031                                           Vector3D<Precision> const &localdir, NavigationState const *in_state,
0032                                           NavigationState * /*out_state*/, Precision &step,
0033                                           VPlacedVolume const *&hitcandidate) const override
0034   {
0035     // iterate over all daughters
0036     auto *daughters = lvol->GetDaughtersp();
0037     auto ndaughters = daughters->size();
0038     for (decltype(ndaughters) d = 0; d < ndaughters; ++d) {
0039       auto daughter = daughters->operator[](d);
0040 //    previous distance becomes step estimate, distance to daughter returned in workspace
0041 // SW: this makes the navigation more robust and it appears that I have to
0042 // put this at the moment since not all shapes respond yet with a negative distance if
0043 // the point is actually inside the daughter
0044 #ifdef CHECKCONTAINS
0045       bool contains = daughter->Contains(localpoint);
0046       if (!contains) {
0047 #endif
0048         Precision ddistance = daughter->DistanceToIn(localpoint, localdir, step);
0049
0050         // if distance is negative; we are inside that daughter and should relocate
0051         // unless distance is minus infinity
0052         const bool valid = (ddistance < step && !IsInf(ddistance)) &&
0053                            !((ddistance <= 0.) && in_state && in_state->GetLastExited() == daughter);
0054         hitcandidate = valid ? daughter : hitcandidate;
0055         step         = valid ? ddistance : step;
0056 #ifdef CHECKCONTAINS
0057       } else {
0058         std::cerr << " INDA "
0059                   << " contained in daughter " << daughter << " - inside = " << daughter->Inside(localpoint)
0060                   << " , distToIn(p,v,s) = " << daughter->DistanceToIn(localpoint, localdir, step) << " \n";
0061
0062         std::cerr << " INDA ";
0063         step         = -1.;
0064         hitcandidate = daughter;
0065         break;
0066       }
0067 #endif
0068     }
0069     return false;
0070   }
0071
0072   VECGEOM_FORCE_INLINE
0073   VECCORE_ATT_HOST_DEVICE
0074   virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0075                                           Vector3D<Precision> const &localdir, VPlacedVolume const *blocked,
0076                                           Precision &step, VPlacedVolume const *&hitcandidate) const override
0077   {
0078     //  New Implementation JA 2021.03.18
0079     static const double kMinExitingCos = 1.e-3;
0080     VPlacedVolume const *excludedVol   = nullptr;
0081     if (blocked) {
0082       Vector3D<Precision> normal;
0083       blocked->Normal(localpoint, normal);
0084       if (normal.Dot(localdir) >= kMinExitingCos) {
0085         excludedVol = blocked;
0086       }
0087     }
0088
0089     // iterate over all daughters
0090     auto *daughters = lvol->GetDaughtersp();
0091     auto ndaughters = daughters->size();
0092     for (decltype(ndaughters) d = 0; d < ndaughters; ++d) {
0093       auto daughter = daughters->operator[](d);
0094 //    previous distance becomes step estimate, distance to daughter returned in workspace
0095 // SW: this makes the navigation more robust and it appears that I have to
0096 // put this at the moment since not all shapes respond yet with a negative distance if
0097 // the point is actually inside the daughter
0098 #ifdef CHECKCONTAINS
0099       bool contains = daughter->Contains(localpoint);
0100       if (!contains) {
0101 #endif
0102         if (daughter != excludedVol) {
0103           Precision ddistance = daughter->DistanceToIn(localpoint, localdir, step);
0104
0105           // if distance is negative; we are inside that daughter and should relocate
0106           // unless distance is minus infinity
0107           const bool valid = (ddistance < step && !IsInf(ddistance));
0108           hitcandidate     = valid ? daughter : hitcandidate;
0109           step             = valid ? ddistance : step;
0110         }
0111 #ifdef CHECKCONTAINS
0112       } else {
0113         std::cerr << " INDA: contained in daughter " << daughter << " - inside = " << daughter->Inside(localpoint)
0114                   << " , distToIn(p,v,s) = " << daughter->DistanceToIn(localpoint, localdir, step) << " \n";
0115         step         = -1.;
0116         hitcandidate = daughter;
0117         break;
0118       }
0119 #endif
0120     }
0121     // assert(false); --- Was not implemented before
0122     return false;
0123   }
0124
0125   // Vector specialization for NewSimpleNavigator
0126   // TODO: unify with scalar use case
0127   template <typename T, unsigned int ChunkSize>
0128   VECCORE_ATT_HOST_DEVICE
0129   static void DaughterIntersectionsLooper(VNavigator const * /*nav*/, LogicalVolume const *lvol,
0130                                           Vector3D<T> const &localpoint, Vector3D<T> const &localdir,
0131                                           NavigationState const* const* /*in_states*/, NavigationState ** /*out_states*/,
0132                                           unsigned int from_index, Precision *out_steps,
0133                                           VPlacedVolume const *hitcandidates[ChunkSize])
0134   {
0135     // dispatch to vector implementation
0136     // iterate over all daughters
0137     // using Real_v = typename Backend::Real_v;
0138     // using Bool_v = Real_v::Mask;
0139     T step(vecCore::FromPtr<T>(out_steps + from_index));
0140     const auto *daughters = lvol->GetDaughtersp();
0141     auto ndaughters       = daughters->size();
0142     for (decltype(ndaughters) d = 0; d < ndaughters; ++d) {
0143       auto daughter = daughters->operator[](d);
0144       // SW: this makes the navigation more robust and it appears that I have to
0145       // put this at the moment since not all shapes respond yet with a negative distance if
0146       // the point is actually inside the daughter
0147       //  #ifdef CHECKCONTAINS
0148       //        Bool_v contains = daughter->Contains(localpoint);
0149       //        if (!contains) {
0150       //  #endif
0151       const T ddistance = daughter->DistanceToIn(localpoint, localdir, step);
0152
0153       // if distance is negative; we are inside that daughter and should relocate
0154       // unless distance is minus infinity
0155       const auto valid = (ddistance < step && !IsInf(ddistance));
0156
0157       // serial treatment of hit candidate until I find a better way
0158       // we might do this using a cast to a double vector with subsequent masked assignment
0159       if (!vecCore::MaskEmpty(valid)) {
0160         for (unsigned int i = 0 /*valid.firstOne()*/; i < ChunkSize; ++i) {
0161           hitcandidates[i] = vecCore::MaskLaneAt(valid, i) ? daughter : hitcandidates[i];
0162         }
0163
0164         vecCore::MaskedAssign(step, valid, ddistance);
0165         //  #ifdef CHECKCONTAINS
0166         //        } else {
0167         //          std::cerr << " INDA ";
0168         //          step = -1.;
0169         //          hitcandidate = daughter;
0170         //          break;
0171         //        }
0172         //  #endif
0173       }
0174     }
0175     vecCore::Store(step, out_steps + from_index);
0176   }
0177
0178   // Vector specialization for NewSimpleNavigator
0179   // TODO: unify with scalar use case
0180   template <typename T, unsigned int ChunkSize>
0181   VECCORE_ATT_HOST_DEVICE
0182   static void DaughterIntersectionsLooper(VNavigator const * /*nav*/, LogicalVolume const *lvol,
0183                                           Vector3D<T> const &localpoint, Vector3D<T> const &localdir,
0184                                           NavStatePool const &/*in_states*/, NavStatePool &/*out_states*/,
0185                                           unsigned int from_index, Precision *out_steps,
0186                                           VPlacedVolume const *hitcandidates[ChunkSize])
0187   {
0188     // dispatch to vector implementation
0189     // iterate over all daughters
0190     // using Real_v = typename Backend::Real_v;
0191     // using Bool_v = Real_v::Mask;
0192     T step(vecCore::FromPtr<T>(out_steps + from_index));
0193     const auto *daughters = lvol->GetDaughtersp();
0194     auto ndaughters       = daughters->size();
0195     for (decltype(ndaughters) d = 0; d < ndaughters; ++d) {
0196       auto daughter = daughters->operator[](d);
0197       // SW: this makes the navigation more robust and it appears that I have to
0198       // put this at the moment since not all shapes respond yet with a negative distance if
0199       // the point is actually inside the daughter
0200       //  #ifdef CHECKCONTAINS
0201       //        Bool_v contains = daughter->Contains(localpoint);
0202       //        if (!contains) {
0203       //  #endif
0204       const T ddistance = daughter->DistanceToIn(localpoint, localdir, step);
0205
0206       // if distance is negative; we are inside that daughter and should relocate
0207       // unless distance is minus infinity
0208       const auto valid = (ddistance < step && !IsInf(ddistance));
0209
0210       // serial treatment of hit candidate until I find a better way
0211       // we might do this using a cast to a double vector with subsequent masked assignment
0212       if (!vecCore::MaskEmpty(valid)) {
0213         for (unsigned int i = 0 /*valid.firstOne()*/; i < ChunkSize; ++i) {
0214           hitcandidates[i] = vecCore::MaskLaneAt(valid, i) ? daughter : hitcandidates[i];
0215         }
0216
0217         vecCore::MaskedAssign(step, valid, ddistance);
0218         //  #ifdef CHECKCONTAINS
0219         //        } else {
0220         //          std::cerr << " INDA ";
0221         //          step = -1.;
0222         //          hitcandidate = daughter;
0223         //          break;
0224         //        }
0225         //  #endif
0226       }
0227     }
0228     vecCore::Store(step, out_steps + from_index);
0229   }
0230
0231   template <typename T, unsigned int ChunkSize>
0232   VECCORE_ATT_HOST_DEVICE
0233   static void DaughterIntersectionsLooper(VNavigator const *nav, LogicalVolume const *lvol,
0234                                           Vector3D<T> const &localpoint, Vector3D<T> const &localdir,
0235                                           NavigationState const *const *in_states, unsigned int from_index,
0236                                           Precision *out_steps, VPlacedVolume const *hitcandidates[ChunkSize])
0237   {
0238     NewSimpleNavigator<MotherIsConvex>::template DaughterIntersectionsLooper<T, ChunkSize>(
0239         nav, lvol, localpoint, localdir, in_states, nullptr, from_index, out_steps, hitcandidates);
0240   }
0241
0242 #ifndef VECCORE_CUDA
0243   static VNavigator *Instance()
0244   {
0245     static NewSimpleNavigator instance;
0246     return &instance;
0247   }
0248 #else
0249   VECCORE_ATT_DEVICE
0250   static VNavigator *Instance();
0251 #endif
0252
0253   static constexpr const char *gClassNameString = "NewSimpleNavigator";
0254   typedef SimpleSafetyEstimator SafetyEstimator_t;
0255 }; // end of class
0256 }
0257 } // end namespace
0258
0259 #endif /* NAVIGATION_NEWSIMPLENAVIGATOR_H_ */