VecGeom/navigation/SimpleABBoxNavigator.h

0001 /*
0002  * SimpleABBoxNavigator.h
0003  *
0004  *  Created on: Nov 23, 2015
0005  *      Author: swenzel
0006  */
0007
0008 #ifndef NAVIGATION_SIMPLEABBOXNAVIGATOR_H_
0009 #define NAVIGATION_SIMPLEABBOXNAVIGATOR_H_
0010
0011 #include "VNavigator.h"
0012 #include "SimpleABBoxSafetyEstimator.h"
0013 #include "VecGeom/management/ABBoxManager.h"
0014
0015 namespace vecgeom {
0016 inline namespace VECGEOM_IMPL_NAMESPACE {
0017
0018 // A basic implementation of a navigator which is SIMD accelerated by using a flat aligned bounding box list
0019 template <bool MotherIsConvex = false>
0020 class SimpleABBoxNavigator : public VNavigatorHelper<SimpleABBoxNavigator<MotherIsConvex>, MotherIsConvex> {
0021
0022 private:
0023   ABBoxManager &fABBoxManager;
0024   SimpleABBoxNavigator()
0025       : VNavigatorHelper<SimpleABBoxNavigator<MotherIsConvex>, MotherIsConvex>(SimpleABBoxSafetyEstimator::Instance()),
0026         fABBoxManager(ABBoxManager::Instance())
0027   {
0028   }
0029
0030   // convert index to physical daugher
0031   VPlacedVolume const *LookupDaughter(LogicalVolume const *lvol, int id) const
0032   {
0033     assert(id >= 0 && "access with negative index");
0034     assert(size_t(id) < lvol->GetDaughtersp()->size() && "access beyond size of daughterlist ");
0035     return lvol->GetDaughtersp()->operator[](id);
0036   }
0037
0038   // a simple sort class (based on insertionsort)
0039   // template <typename T, typename Cmp>
0040   static void insertionsort(ABBoxManager::BoxIdDistancePair_t *arr, unsigned int N)
0041   {
0042     for (unsigned short i = 1; i < N; ++i) {
0043       ABBoxManager::BoxIdDistancePair_t value = arr[i];
0044       short hole                              = i;
0045
0046       for (; hole > 0 && value.second < arr[hole - 1].second; --hole)
0047         arr[hole] = arr[hole - 1];
0048
0049       arr[hole] = value;
0050     }
0051   }
0052
0053   // vector version
0054   size_t GetHitCandidates_v(LogicalVolume const * /*lvol*/, Vector3D<Precision> const &point,
0055                             Vector3D<Precision> const &dir, ABBoxManager::ABBoxContainer_v const &corners, size_t size,
0056                             ABBoxManager::BoxIdDistancePair_t *hitlist) const
0057   {
0058     size_t vecsize  = size;
0059     size_t hitcount = 0;
0060     Vector3D<float> invdirfloat(1.f / (float)dir.x(), 1.f / (float)dir.y(), 1.f / (float)dir.z());
0061     Vector3D<float> pfloat((float)point.x(), (float)point.y(), (float)point.z());
0062     int sign[3];
0063     sign[0]       = invdirfloat.x() < 0;
0064     sign[1]       = invdirfloat.y() < 0;
0065     sign[2]       = invdirfloat.z() < 0;
0066     using Float_v = ABBoxManager::Float_v;
0067     using Bool_v  = vecCore::Mask_v<Float_v>;
0068     for (size_t box = 0; box < vecsize; ++box) {
0069       Float_v distance = BoxImplementation::IntersectCachedKernel2<Float_v, float>(
0070           &corners[2 * box], pfloat, invdirfloat, sign[0], sign[1], sign[2], 0, InfinityLength<float>());
0071       Bool_v hit = distance < InfinityLength<float>();
0072       if (!vecCore::MaskEmpty(hit)) {
0073         constexpr auto kVS = vecCore::VectorSize<Float_v>();
0074
0075         // VecCore does not have firstOne() function; iterating from zero
0076         // consider putting a firstOne into vecCore or in VecGeom
0077         for (size_t i = 0; i < kVS; ++i) {
0078           if (vecCore::MaskLaneAt(hit, i)) {
0079             assert(hitcount < VECGEOM_MAXDAUGHTERS);
0080             hitlist[hitcount] = (ABBoxManager::BoxIdDistancePair_t(box * kVS + i, vecCore::LaneAt(distance, i)));
0081             hitcount++;
0082           }
0083         }
0084       }
0085     }
0086     return hitcount;
0087   }
0088
0089 public:
0090   // we provide hit detection on the local level and reuse the generic implementations from
0091   // VNavigatorHelper<SimpleABBoxNavigator>
0092
0093   VECGEOM_FORCE_INLINE
0094   virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0095                                           Vector3D<Precision> const &localdir, NavigationState const *in_state,
0096                                           NavigationState * /*out_state*/, Precision &step,
0097                                           VPlacedVolume const *&hitcandidate) const override
0098   {
0099     // The following construct reserves stackspace for objects
0100     // of type IdDistPair_t WITHOUT initializing those objects
0101     using IdDistPair_t = ABBoxManager::BoxIdDistancePair_t;
0102     char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0103     IdDistPair_t *hitlist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0104
0105     if (lvol->GetDaughtersp()->size() == 0) return false;
0106
0107     int size;
0108     ABBoxManager::ABBoxContainer_v bboxes = fABBoxManager.GetABBoxes_v(lvol, size);
0109     auto ncandidates                      = GetHitCandidates_v(lvol, localpoint, localdir, bboxes, size, hitlist);
0110
0111     // sort candidates according to their bounding volume hit distance
0112     insertionsort(hitlist, ncandidates);
0113
0114     for (size_t index = 0; index < ncandidates; ++index) {
0115       auto &hitbox                   = hitlist[index];
0116       VPlacedVolume const *candidate = LookupDaughter(lvol, hitbox.first);
0117
0118       // only consider those hitboxes which are within potential reach of this step
0119       if (!(step < hitbox.second)) {
0120         //      std::cerr << "checking id " << hitbox.first << " at box distance " << hitbox.second << "\n";
0121         //           if( hitbox.second < 0 ){
0122         //            //   std::cerr << "funny2\n";
0123         //              bool checkindaughter = candidate->Contains( localpoint );
0124         //              if( checkindaughter == true ){
0125         //                  std::cerr << "funny\n";
0126         //
0127         //                  // need to relocate
0128         //                  step = 0;
0129         //                  hitcandidate = candidate;
0130         //                  // THE ALTERNATIVE WOULD BE TO PUSH THE CURRENT STATE AND RETURN DIRECTLY
0131         //                  break;
0132         //              }
0133         //          }
0134         Precision ddistance = candidate->DistanceToIn(localpoint, localdir, step);
0135 #ifdef VERBOSE
0136         std::cerr << "distance to " << candidate->GetLabel() << " is " << ddistance << "\n";
0137 #endif
0138         const auto valid = !IsInf(ddistance) && ddistance < step &&
0139                            !((ddistance <= 0.) && in_state && in_state->GetLastExited() == candidate);
0140         hitcandidate = valid ? candidate : hitcandidate;
0141         step         = valid ? ddistance : step;
0142       } else {
0143         break;
0144       }
0145     }
0146     return false;
0147   }
0148
0149   virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0150                                           Vector3D<Precision> const &localdir, VPlacedVolume const *blocked,
0151                                           Precision &step, VPlacedVolume const *&hitcandidate) const override
0152   {
0153     // The following construct reserves stackspace for objects
0154     // of type IdDistPair_t WITHOUT initializing those objects
0155     using IdDistPair_t = ABBoxManager::BoxIdDistancePair_t;
0156     char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0157     IdDistPair_t *hitlist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0158
0159     if (lvol->GetDaughtersp()->size() == 0) return false;
0160
0161     int size;
0162     ABBoxManager::ABBoxContainer_v bboxes = fABBoxManager.GetABBoxes_v(lvol, size);
0163     auto ncandidates                      = GetHitCandidates_v(lvol, localpoint, localdir, bboxes, size, hitlist);
0164
0165     // sort candidates according to their bounding volume hit distance
0166     insertionsort(hitlist, ncandidates);
0167
0168     for (size_t index = 0; index < ncandidates; ++index) {
0169       auto &hitbox                   = hitlist[index];
0170       VPlacedVolume const *candidate = LookupDaughter(lvol, hitbox.first);
0171
0172       // only consider those hitboxes which are within potential reach of this step
0173       if (hitbox.second <= step) { // !(step < hitbox.second)) {
0174 #ifdef VERBOSE
0175         std::cerr << "SimpleABBoxNavigator: checking id " << hitbox.first << " at box distance " << hitbox.second
0176                   << "\n";
0177 #endif
0178         //           if( hitbox.second < 0 ){
0179         //            //   std::cerr << "hit dist < 0. \n";
0180         //              bool checkindaughter = candidate->Contains( localpoint );
0181         //              if( checkindaughter == true ){
0182         //                  std::cerr << " -- unexpected : in daughter vol \n";
0183         //
0184         //                  // need to relocate
0185         //                  step = 0;
0186         //                  hitcandidate = candidate;
0187         //                  // THE ALTERNATIVE WOULD BE TO PUSH THE CURRENT STATE AND RETURN DIRECTLY
0188         //                  break;
0189         //              }
0190         //          }
0191         Precision ddistance = candidate->DistanceToIn(localpoint, localdir, step);
0192         Vector3D<Precision> normal; // To reuse in printing below - else move it into 'if'
0193 #ifdef VERBOSE
0194         std::cerr << " . Distance to " << candidate->GetLabel() << " is " << ddistance;
0195 #endif
0196         if (ddistance <= 0. /* && candidate == blocked */) {
0197           candidate->Normal(localpoint, normal);
0198 #ifdef VERBOSE
0199           std::cerr << "  SimpleABBoxNavigator:  Negative daughter dist = " << ddistance
0200                     << (blocked == candidate ? " Is " : " Not ") << "Blocked. "
0201                     << " normal= " << normal << "  normal.dir= " << normal.Dot(localdir) << "\n";
0202 #endif
0203         }
0204         const auto valid = !IsInf(ddistance) && ddistance < step &&
0205                            !((ddistance <= 0.) && (blocked == candidate || normal.Dot(localdir) > 0.0));
0206
0207         hitcandidate = valid ? candidate : hitcandidate;
0208         step         = valid ? ddistance : step;
0209       } else {
0210         break;
0211       }
0212     }
0213     return false;
0214   }
0215
0216   static VNavigator *Instance()
0217   {
0218     static SimpleABBoxNavigator instance;
0219     return &instance;
0220   }
0221
0222   static constexpr const char *gClassNameString = "SimpleABBoxNavigator";
0223   typedef SimpleABBoxSafetyEstimator SafetyEstimator_t;
0224 }; // end of class
0225 }
0226 } // end namespace
0227
0228 #endif /* NAVIGATION_SIMPLEABBOXNAVIGATOR_H_ */