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<Precision> &fABBoxManager;
0024   SimpleABBoxNavigator()
0025       : VNavigatorHelper<SimpleABBoxNavigator<MotherIsConvex>, MotherIsConvex>(SimpleABBoxSafetyEstimator::Instance()),
0026         fABBoxManager(ABBoxManager<Precision>::Instance())
0027   {
0028   }
0029
0030   // convert index to physical daugher
0031   VPlacedVolume const *LookupDaughter(LogicalVolume const *lvol, int id) const
0032   {
0033     VECGEOM_VALIDATE(id >= 0, << "access with negative index");
0034     VECGEOM_VALIDATE(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<Precision>::BoxIdDistancePair_t *arr, unsigned int N)
0041   {
0042     for (unsigned short i = 1; i < N; ++i) {
0043       ABBoxManager<Precision>::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<Precision>::ABBoxContainer_v const &corners,
0056                             size_t size, ABBoxManager<Precision>::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<Precision>::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             VECGEOM_ASSERT(hitcount < VECGEOM_MAXDAUGHTERS);
0080             hitlist[hitcount] =
0081                 (ABBoxManager<Precision>::BoxIdDistancePair_t(box * kVS + i, vecCore::LaneAt(distance, i)));
0082             hitcount++;
0083           }
0084         }
0085       }
0086     }
0087     return hitcount;
0088   }
0089
0090 public:
0091   // we provide hit detection on the local level and reuse the generic implementations from
0092   // VNavigatorHelper<SimpleABBoxNavigator>
0093
0094   VECGEOM_FORCE_INLINE
0095   virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0096                                           Vector3D<Precision> const &localdir, NavigationState const *in_state,
0097                                           NavigationState * /*out_state*/, Precision &step,
0098                                           VPlacedVolume const *&hitcandidate) const override
0099   {
0100     // The following construct reserves stackspace for objects
0101     // of type IdDistPair_t WITHOUT initializing those objects
0102     using IdDistPair_t = ABBoxManager<Precision>::BoxIdDistancePair_t;
0103     char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0104     IdDistPair_t *hitlist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0105
0106     if (lvol->GetDaughtersp()->size() == 0) return false;
0107
0108     int size;
0109     ABBoxManager<Precision>::ABBoxContainer_v bboxes = fABBoxManager.GetABBoxes_v(lvol, size);
0110     auto ncandidates = GetHitCandidates_v(lvol, localpoint, localdir, bboxes, size, hitlist);
0111
0112     // sort candidates according to their bounding volume hit distance
0113     insertionsort(hitlist, ncandidates);
0114
0115     for (size_t index = 0; index < ncandidates; ++index) {
0116       auto &hitbox                   = hitlist[index];
0117       VPlacedVolume const *candidate = LookupDaughter(lvol, hitbox.first);
0118       if (in_state && in_state->GetLastExited() == candidate) continue;
0119
0120       // only consider those hitboxes which are within potential reach of this step
0121       if (!(step < hitbox.second)) {
0122         //      std::cerr << "checking id " << hitbox.first << " at box distance " << hitbox.second << "\n";
0123         //           if( hitbox.second < 0 ){
0124         //            //   std::cerr << "funny2\n";
0125         //              bool checkindaughter = candidate->Contains( localpoint );
0126         //              if( checkindaughter == true ){
0127         //                  std::cerr << "funny\n";
0128         //
0129         //                  // need to relocate
0130         //                  step = 0;
0131         //                  hitcandidate = candidate;
0132         //                  // THE ALTERNATIVE WOULD BE TO PUSH THE CURRENT STATE AND RETURN DIRECTLY
0133         //                  break;
0134         //              }
0135         //          }
0136         Precision ddistance = candidate->DistanceToIn(localpoint, localdir, step);
0137 #ifdef VERBOSE
0138         std::cerr << "distance to " << candidate->GetLabel() << " is " << ddistance << "\n";
0139 #endif
0140         const auto valid = ddistance < step;
0141         hitcandidate     = valid ? candidate : hitcandidate;
0142         step             = valid ? ddistance : step;
0143       } else {
0144         break;
0145       }
0146     }
0147     return false;
0148   }
0149
0150   virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0151                                           Vector3D<Precision> const &localdir, VPlacedVolume const *blocked,
0152                                           Precision &step, VPlacedVolume const *&hitcandidate) const override
0153   {
0154     // The following construct reserves stackspace for objects
0155     // of type IdDistPair_t WITHOUT initializing those objects
0156     using IdDistPair_t = ABBoxManager<Precision>::BoxIdDistancePair_t;
0157     char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0158     IdDistPair_t *hitlist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0159
0160     if (lvol->GetDaughtersp()->size() == 0) return false;
0161
0162     int size;
0163     ABBoxManager<Precision>::ABBoxContainer_v bboxes = fABBoxManager.GetABBoxes_v(lvol, size);
0164     auto ncandidates = GetHitCandidates_v(lvol, localpoint, localdir, bboxes, size, hitlist);
0165
0166     // sort candidates according to their bounding volume hit distance
0167     insertionsort(hitlist, ncandidates);
0168
0169     for (size_t index = 0; index < ncandidates; ++index) {
0170       auto &hitbox                   = hitlist[index];
0171       VPlacedVolume const *candidate = LookupDaughter(lvol, hitbox.first);
0172
0173       // only consider those hitboxes which are within potential reach of this step
0174       if (hitbox.second <= step) { // !(step < hitbox.second)) {
0175         Precision ddistance = candidate->DistanceToIn(localpoint, localdir, step);
0176         Vector3D<Precision> normal; // To reuse in printing below - else move it into 'if'
0177         if (ddistance <= 0.) candidate->Normal(localpoint, normal);
0178         const auto valid = !IsInf(ddistance) && ddistance < step &&
0179                            !((ddistance <= 0.) && (blocked == candidate || normal.Dot(localdir) > 0.0));
0180
0181         hitcandidate = valid ? candidate : hitcandidate;
0182         step         = valid ? ddistance : step;
0183       } else {
0184         break;
0185       }
0186     }
0187     return false;
0188   }
0189
0190   static VNavigator *Instance()
0191   {
0192     static SimpleABBoxNavigator instance;
0193     return &instance;
0194   }
0195
0196   static constexpr const char *gClassNameString = "SimpleABBoxNavigator";
0197   typedef SimpleABBoxSafetyEstimator SafetyEstimator_t;
0198 }; // end of class
0199 } // namespace VECGEOM_IMPL_NAMESPACE
0200 } // namespace vecgeom
0201
0202 #endif /* NAVIGATION_SIMPLEABBOXNAVIGATOR_H_ */