VecGeom/navigation/NavStateIndex.h

0001 /// \file NavStateIndex
0002 /// \author Andrei Gheata (andrei.gheata@cern.ch)
0003 /// \date 12.03.2014
0004
0005 #ifndef VECGEOM_NAVIGATION_NAVSTATEINDEX_H_
0006 #define VECGEOM_NAVIGATION_NAVSTATEINDEX_H_
0007
0008 #include "VecGeom/base/Config.h"
0009 #include "VecGeom/base/Global.h"
0010 #include "VecGeom/base/Transformation3D.h"
0011 #include "VecGeom/volumes/PlacedVolume.h"
0012 #include "VecGeom/management/GeoManager.h"
0013
0014 #ifdef VECGEOM_ENABLE_CUDA
0015 #include "VecGeom/management/CudaManager.h"
0016 #endif
0017
0018 #include <iostream>
0019 #include <string>
0020
0021 class TGeoBranchArray;
0022
0023 namespace vecgeom {
0024 inline namespace VECGEOM_IMPL_NAMESPACE {
0025
0026 /**
0027  * A class describing a current geometry state based on a single index
0028  * likely there will be such an object for each particle/track currently treated.
0029  */
0030 class NavStateIndex {
0031 public:
0032   using Value_t = unsigned long;
0033
0034 private:
0035   NavIndex_t fNavInd     = 0;     ///< Navigation state index
0036   NavIndex_t fLastExited = 0;     ///< Navigation state index of the last exited state
0037   bool fOnBoundary       = false; ///< flag indicating whether track is on boundary of the "Top()" placed volume
0038
0039 public:
0040   VECCORE_ATT_HOST_DEVICE
0041   NavStateIndex(NavIndex_t nav_ind = 0) { fNavInd = nav_ind; }
0042
0043   VECGEOM_FORCE_INLINE
0044   VECCORE_ATT_HOST_DEVICE
0045   static unsigned char GetMaxLevel()
0046   {
0047 #ifdef VECCORE_CUDA_DEVICE_COMPILATION
0048     return vecgeom::globaldevicegeomdata::gMaxDepth;
0049 #else
0050     return (unsigned char)GeoManager::Instance().getMaxDepth();
0051 #endif
0052   }
0053
0054   // Static accessors
0055   VECCORE_ATT_HOST_DEVICE
0056   static NavStateIndex *MakeInstance(int)
0057   {
0058     // MaxLevel is 'zero' based (i.e. maxlevel==0 requires one value)
0059     return new NavStateIndex();
0060   }
0061
0062   VECCORE_ATT_HOST_DEVICE
0063   static NavStateIndex *MakeCopy(NavStateIndex const &other) { return new NavStateIndex(other); }
0064
0065   VECCORE_ATT_HOST_DEVICE
0066   static NavStateIndex *MakeInstanceAt(int, void *addr) { return new (addr) NavStateIndex(); }
0067
0068   VECCORE_ATT_HOST_DEVICE
0069   static NavStateIndex *MakeCopy(NavStateIndex const &other, void *addr) { return new (addr) NavStateIndex(other); }
0070
0071   VECCORE_ATT_HOST_DEVICE
0072   static void ReleaseInstance(NavStateIndex *state)
0073   {
0074     // MaxLevel is 'zero' based (i.e. maxlevel==0 requires one value)
0075     delete state;
0076   }
0077
0078   // returns the size in bytes of a NavStateIndex object with internal
0079   // path depth maxlevel
0080   VECCORE_ATT_HOST_DEVICE
0081   static size_t SizeOfInstance(int)
0082   {
0083     // MaxLevel is 'zero' based (i.e. maxlevel==0 requires one value)
0084     return sizeof(NavStateIndex);
0085   }
0086
0087   // returns the size in bytes of a NavStateIndex object with internal
0088   // path depth maxlevel -- including space needed for padding to next aligned object
0089   // of same kind
0090   VECCORE_ATT_HOST_DEVICE
0091   static size_t SizeOfInstanceAlignAware(int)
0092   {
0093     // MaxLevel is 'zero' based (i.e. maxlevel==0 requires one value)
0094     return sizeof(NavStateIndex);
0095   }
0096
0097   VECCORE_ATT_HOST_DEVICE
0098   VECGEOM_FORCE_INLINE
0099   NavIndex_t GetNavIndex() const { return fNavInd; }
0100
0101   VECCORE_ATT_HOST_DEVICE
0102   int GetObjectSize() const { return (int)sizeof(NavStateIndex); }
0103
0104   VECCORE_ATT_HOST_DEVICE
0105   static size_t SizeOf(size_t) { return sizeof(NavStateIndex); }
0106
0107   VECCORE_ATT_HOST_DEVICE
0108   void CopyTo(NavStateIndex *other) const { *other = *this; }
0109
0110   // copies a fixed and predetermined number of bytes
0111   // might be useful for specialized navigators which know the depth + SizeOf in advance
0112   // N is number of bytes to be copied and can be obtained by a prior call to constexpr NavStateIndex::SizeOf( ... );
0113   template <size_t N>
0114   void CopyToFixedSize(NavStateIndex *other) const
0115   {
0116     *other = *this;
0117   }
0118
0119   VECCORE_ATT_HOST_DEVICE
0120   VECGEOM_FORCE_INLINE
0121   static NavIndex_t const *NavIndAddr(NavIndex_t nav_ind)
0122   {
0123 #ifdef VECCORE_CUDA_DEVICE_COMPILATION
0124     // checking here for NVCC_DEVICE since the global variable globaldevicegeomgata::gCompact...
0125     // is marked __device__ and can only be compiled within device compiler passes
0126     assert(vecgeom::globaldevicegeomdata::gNavIndex != nullptr);
0127     return &vecgeom::globaldevicegeomdata::gNavIndex[nav_ind];
0128 #else
0129 #ifndef VECCORE_CUDA
0130     assert(vecgeom::GeoManager::gNavIndex != nullptr);
0131     return &vecgeom::GeoManager::gNavIndex[nav_ind];
0132 #else
0133     // this is the case when we compile with nvcc for host side
0134     // (failed previously due to undefined symbol vecgeom::cuda::GeoManager::gCompactPlacedVolBuffer)
0135     assert(false && "reached unimplement code");
0136     return nullptr;
0137 #endif
0138 #endif
0139   }
0140
0141   VECCORE_ATT_HOST_DEVICE
0142   VECGEOM_FORCE_INLINE
0143   static NavIndex_t NavInd(NavIndex_t nav_ind) { return *NavIndAddr(nav_ind); }
0144
0145   VECCORE_ATT_HOST_DEVICE
0146   VECGEOM_FORCE_INLINE
0147   static VPlacedVolume const *ToPlacedVolume(size_t index)
0148   {
0149 #ifdef VECCORE_CUDA_DEVICE_COMPILATION
0150     // checking here for NVCC_DEVICE since the global variable globaldevicegeomgata::gCompact...
0151     // is marked __device__ and can only be compiled within device compiler passes
0152     assert(vecgeom::globaldevicegeomdata::gCompactPlacedVolBuffer != nullptr);
0153     return &vecgeom::globaldevicegeomdata::gCompactPlacedVolBuffer[index];
0154 #else
0155 #ifndef VECCORE_CUDA
0156     assert(vecgeom::GeoManager::gCompactPlacedVolBuffer == nullptr ||
0157            vecgeom::GeoManager::gCompactPlacedVolBuffer[index].id() == index);
0158     return &vecgeom::GeoManager::gCompactPlacedVolBuffer[index];
0159 #else
0160     // this is the case when we compile with nvcc for host side
0161     // (failed previously due to undefined symbol vecgeom::cuda::GeoManager::gCompactPlacedVolBuffer)
0162     assert(false && "reached unimplement code");
0163     (void)index; // avoid unused parameter warning.
0164     return nullptr;
0165 #endif
0166 #endif
0167   }
0168
0169   VECCORE_ATT_HOST_DEVICE
0170   VECGEOM_FORCE_INLINE
0171   static unsigned short GetNdaughtersImpl(NavIndex_t nav_ind)
0172   {
0173     constexpr unsigned int kOffsetNd = 2 * sizeof(NavIndex_t) + 2;
0174     auto content_nd                  = (unsigned short *)((unsigned char *)(NavIndAddr(nav_ind)) + kOffsetNd);
0175     return *content_nd;
0176   }
0177
0178   VECGEOM_FORCE_INLINE
0179   VECCORE_ATT_HOST_DEVICE
0180   static unsigned char GetLevelImpl(NavIndex_t nav_ind)
0181   {
0182     constexpr unsigned int kOffsetLevel = 2 * sizeof(NavIndex_t);
0183     auto content_level                  = (unsigned char *)(NavIndAddr(nav_ind)) + kOffsetLevel;
0184     return *content_level;
0185   }
0186
0187   VECCORE_ATT_HOST_DEVICE
0188   VECGEOM_FORCE_INLINE
0189   static NavIndex_t GetNavIndexImpl(NavIndex_t nav_ind, int level)
0190   {
0191     int up            = GetLevelImpl(nav_ind) - level;
0192     NavIndex_t mother = nav_ind;
0193     while (mother && up--)
0194       mother = NavInd(mother);
0195     return mother;
0196   }
0197
0198   VECGEOM_FORCE_INLINE
0199   VECCORE_ATT_HOST_DEVICE
0200   static NavIndex_t PopImpl(NavIndex_t nav_ind) { return (nav_ind > 0) ? NavInd(nav_ind) : 0; }
0201
0202   VECGEOM_FORCE_INLINE
0203   VECCORE_ATT_HOST_DEVICE
0204   static NavIndex_t PushImpl(NavIndex_t nav_ind, VPlacedVolume const *v)
0205   {
0206     return (nav_ind > 0) ? NavInd(nav_ind + 3 + v->GetChildId()) : 1;
0207   }
0208
0209   VECGEOM_FORCE_INLINE
0210   VECCORE_ATT_HOST_DEVICE
0211   static VPlacedVolume const *TopImpl(NavIndex_t nav_ind)
0212   {
0213     return (nav_ind > 0) ? ToPlacedVolume(NavInd(nav_ind + 1)) : nullptr;
0214   }
0215
0216   VECCORE_ATT_HOST_DEVICE
0217   static void TopMatrixImpl(NavIndex_t nav_ind, Transformation3D &trans)
0218   {
0219     constexpr unsigned int kOffsetHasm = 2 * sizeof(NavIndex_t) + 1;
0220
0221     unsigned char hasm;
0222     while (true) {
0223       if (nav_ind == 0) return;
0224       hasm            = *((unsigned char *)(NavIndAddr(nav_ind)) + kOffsetHasm);
0225       bool has_matrix = (hasm & 0x04) > 0;
0226       if (has_matrix) break;
0227       auto t = *TopImpl(nav_ind)->GetTransformation();
0228       t.MultiplyFromRight(trans);
0229       trans   = t;
0230       nav_ind = NavInd(nav_ind);
0231     }
0232
0233     if ((hasm & 0x03) == 0) return;
0234     bool has_trans = (hasm & 0x02) > 0;
0235     bool has_rot   = (hasm & 0x01) > 0;
0236     auto nd        = GetNdaughtersImpl(nav_ind);
0237
0238     // Potentially skip one NavIndex_t to ensure alignment of transformation data
0239     auto transformationDataIndex     = nav_ind + 3 + nd + ((nd + 1) & 1);
0240     const bool padTransformationData = (transformationDataIndex * sizeof(NavIndex_t)) % sizeof(::Precision) != 0;
0241     transformationDataIndex += unsigned{padTransformationData};
0242
0243     const auto address = reinterpret_cast<const Precision *>(NavIndAddr(transformationDataIndex));
0244     assert(reinterpret_cast<uintptr_t>(address) % sizeof(Precision) == 0);
0245
0246     Transformation3D t;
0247     t.Set(address, address + 3, has_trans, has_rot);
0248     t.MultiplyFromRight(trans);
0249     trans = t;
0250   }
0251
0252   VECCORE_ATT_HOST_DEVICE
0253   static Vector3D<Precision> GlobalToLocalImpl(NavIndex_t nav_ind, Vector3D<Precision> const &globalpoint)
0254   {
0255     Transformation3D trans;
0256     TopMatrixImpl(nav_ind, trans);
0257     Vector3D<Precision> local = trans.Transform(globalpoint);
0258     return local;
0259   }
0260
0261   // Intrerface methods
0262   VECGEOM_FORCE_INLINE
0263   VECCORE_ATT_HOST_DEVICE
0264   VPlacedVolume const *GetLastExited() const { return TopImpl(fLastExited); }
0265
0266   VECGEOM_FORCE_INLINE
0267   VECCORE_ATT_HOST_DEVICE
0268   void SetLastExited() { fLastExited = fNavInd; }
0269
0270   VECCORE_ATT_HOST_DEVICE
0271   VECGEOM_FORCE_INLINE
0272   unsigned short GetNdaughters() const { return GetNdaughtersImpl(fNavInd); }
0273
0274   VECGEOM_FORCE_INLINE
0275   VECCORE_ATT_HOST_DEVICE
0276   void Push(VPlacedVolume const *v)
0277   {
0278     // fLastExited = fNavInd;
0279     fNavInd = PushImpl(fNavInd, v);
0280   }
0281
0282   VECGEOM_FORCE_INLINE
0283   VECCORE_ATT_HOST_DEVICE
0284   void Pop()
0285   {
0286     // fLastExited = fNavInd;
0287     fNavInd = PopImpl(fNavInd);
0288   }
0289
0290   VECGEOM_FORCE_INLINE
0291   VECCORE_ATT_HOST_DEVICE
0292   VPlacedVolume const *Top() const { return TopImpl(fNavInd); }
0293
0294   /**
0295    * returns the number of FILLED LEVELS such that
0296    * state.GetNode( state.GetLevel() ) == state.Top()
0297    */
0298   VECGEOM_FORCE_INLINE
0299   VECCORE_ATT_HOST_DEVICE
0300   unsigned char GetLevel() const { return GetLevelImpl(fNavInd); }
0301
0302   /** Compatibility getter for NavigationState interface */
0303   VECGEOM_FORCE_INLINE
0304   VECCORE_ATT_HOST_DEVICE
0305   unsigned char GetCurrentLevel() const { return GetLevel() + 1; }
0306
0307   /**
0308    * Returns the navigation index for a level smaller/equal than the current level.
0309    */
0310   VECCORE_ATT_HOST_DEVICE
0311   VECGEOM_FORCE_INLINE
0312   NavIndex_t GetNavIndex(int level) const { return GetNavIndexImpl(fNavInd, level); }
0313
0314   /**
0315    * Returns the placed volume at a evel smaller/equal than the current level.
0316    */
0317   VECGEOM_FORCE_INLINE
0318   VECCORE_ATT_HOST_DEVICE
0319   VPlacedVolume const *At(int level) const
0320   {
0321     auto parent = GetNavIndexImpl(fNavInd, level);
0322     return (parent > 0) ? ToPlacedVolume(NavInd(parent + 1)) : nullptr;
0323   }
0324
0325   /**
0326    * Returns the index of a placed volume at a evel smaller/equal than the current level.
0327    */
0328   VECGEOM_FORCE_INLINE
0329   VECCORE_ATT_HOST_DEVICE
0330   size_t ValueAt(int level) const
0331   {
0332     auto parent = GetNavIndexImpl(fNavInd, level);
0333     return (parent > 0) ? (size_t)NavInd(parent + 1) : 0;
0334   }
0335
0336   VECCORE_ATT_HOST_DEVICE
0337   void TopMatrix(Transformation3D &trans) const { TopMatrixImpl(fNavInd, trans); }
0338
0339   VECCORE_ATT_HOST_DEVICE
0340   void TopMatrix(int tolevel, Transformation3D &trans) const
0341   {
0342     TopMatrixImpl(GetNavIndexImpl(fNavInd, tolevel), trans);
0343   }
0344
0345   // returning a "delta" transformation that can transform
0346   // coordinates given in reference frame of this->Top() to the reference frame of other->Top()
0347   // simply with otherlocalcoordinate = delta.Transform( thislocalcoordinate )
0348   VECCORE_ATT_HOST_DEVICE
0349   void DeltaTransformation(NavStateIndex const &other, Transformation3D &delta) const;
0350
0351   // VECGEOM_FORCE_INLINE
0352   VECCORE_ATT_HOST_DEVICE
0353   Vector3D<Precision> GlobalToLocal(Vector3D<Precision> const &localpoint) const
0354   {
0355     return GlobalToLocalImpl(fNavInd, localpoint);
0356   }
0357
0358   VECCORE_ATT_HOST_DEVICE
0359   Vector3D<Precision> GlobalToLocal(Vector3D<Precision> const &localpoint, int tolevel) const
0360   {
0361     return GlobalToLocalImpl(GetNavIndexImpl(fNavInd, tolevel), localpoint);
0362   }
0363
0364   /**
0365    * calculates if other navigation state takes a different branch in geometry path or is on same branch
0366    * ( two states are on same branch if one can connect the states just by going upwards or downwards ( or do nothing ))
0367    */
0368   VECCORE_ATT_HOST_DEVICE
0369   int Distance(NavStateIndex const &) const;
0370
0371   // returns a string representation of a (relative) sequence of operations/moves
0372   // that transforms this navigation state into the other navigation state
0373   // example:
0374   // state1 = /0/1/1/
0375   // state2 = /0/2/2/3
0376   // results in string
0377   // "/up/horiz/1/down/2/down/3" with 4 operations "up", "horiz", "down", "down"
0378   // the sequence of moves is the following
0379   // up: /0/1/1 --> /0/1/
0380   // horiz/1 : 0/1 --> /0/2 ( == /0/(1+1) )   "we are hopping from daughter 1 to 2 (which corresponds to a step of 1)"
0381   // down/2 : /0/2 --> /0/2/2   "going further down 2nd daughter"
0382   // down/3 : /0/2/2/3 --> /0/2/2/3  "going further down 2nd daughter"
0383   std::string RelativePath(NavStateIndex const & /*other*/) const;
0384
0385   // functions useful to "serialize" navigationstate
0386   // the Vector-of-Indices basically describes the path on the tree taken from top to bottom
0387   // an index corresponds to a daughter
0388
0389   void GetPathAsListOfIndices(std::list<uint> &indices) const;
0390   void ResetPathFromListOfIndices(VPlacedVolume const *world, std::list<uint> const &indices);
0391
0392   // replaces the volume pointers from CPU volumes in fPath
0393   // to the equivalent pointers on the GPU
0394   // uses the CudaManager to do so
0395   void ConvertToGPUPointers() {}
0396
0397   // replaces the pointers from GPU volumes in fPath
0398   // to the equivalent pointers on the CPU
0399   // uses the CudaManager to do so
0400   void ConvertToCPUPointers() {}
0401
0402   // clear all information
0403   VECGEOM_FORCE_INLINE
0404   VECCORE_ATT_HOST_DEVICE
0405   void Clear()
0406   {
0407     fNavInd     = 0;
0408     fLastExited = 0;
0409     fOnBoundary = false;
0410   }
0411
0412   VECCORE_ATT_HOST_DEVICE
0413   void Print() const;
0414
0415   VECCORE_ATT_HOST_DEVICE
0416   void Dump() const { Print(); }
0417
0418   VECGEOM_FORCE_INLINE
0419   VECCORE_ATT_HOST_DEVICE
0420   bool HasSamePathAsOther(NavStateIndex const &other) const { return (fNavInd == other.fNavInd); }
0421
0422   void printValueSequence(std::ostream & = std::cerr) const;
0423
0424   // calculates a checksum along the path
0425   // can be used (as a quick criterion) to see whether 2 states are same
0426   unsigned long getCheckSum() const { return (unsigned long)fNavInd; }
0427
0428   /**
0429     function returning whether the point (current navigation state) is outside the detector setup
0430   */
0431   VECGEOM_FORCE_INLINE
0432   VECCORE_ATT_HOST_DEVICE
0433   bool IsOutside() const { return (fNavInd == 0); }
0434
0435   VECGEOM_FORCE_INLINE
0436   VECCORE_ATT_HOST_DEVICE
0437   bool IsOnBoundary() const { return fOnBoundary; }
0438
0439   VECGEOM_FORCE_INLINE
0440   VECCORE_ATT_HOST_DEVICE
0441   void SetBoundaryState(bool b) { fOnBoundary = b; }
0442 };
0443
0444 /**
0445  * encodes the geometry path as a concatenated string of ( Value_t ) present in fPath
0446  */
0447 inline void NavStateIndex::printValueSequence(std::ostream &stream) const
0448 {
0449   auto level = GetLevel();
0450   for (int i = 0; i < level + 1; ++i) {
0451     auto pvol = At(i);
0452     if (pvol) stream << "/" << ValueAt(i) << "(" << pvol->GetLabel() << ")";
0453   }
0454 }
0455
0456 } // namespace VECGEOM_IMPL_NAMESPACE
0457 } // namespace vecgeom
0458
0459 #endif // VECGEOM_NAVIGATION_NAVSTATEINDEX_H_