VecGeom/navigation/SimpleABBoxSafetyEstimator.h

0001 /*
0002  * SimpleABBoxSafetyEstimator.h
0003  *
0004  *  Created on: Aug 28, 2015
0005  *      Author: swenzel
0006  */
0007
0008 #ifndef NAVIGATION_SIMPLEABBOXSAFETYESTIMATOR_H_
0009 #define NAVIGATION_SIMPLEABBOXSAFETYESTIMATOR_H_
0010
0011 #include "VecGeom/navigation/VSafetyEstimator.h"
0012 #include "VecGeom/management/ABBoxManager.h"
0013
0014 namespace vecgeom {
0015 inline namespace VECGEOM_IMPL_NAMESPACE {
0016
0017 //! a safety estimator using a (vectorized) search through bounding boxes to exclude certain daughter volumes
0018 //! to talk to
0019 class SimpleABBoxSafetyEstimator : public VSafetyEstimatorHelper<SimpleABBoxSafetyEstimator> {
0020
0021 private:
0022   // we keep a reference to the ABBoxManager ( avoids calling Instance() on this guy all the time )
0023   ABBoxManager &fABBoxManager;
0024
0025   SimpleABBoxSafetyEstimator()
0026       : VSafetyEstimatorHelper<SimpleABBoxSafetyEstimator>(), 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 public:
0039   // helper function calculating some candidate volumes
0040   VECCORE_ATT_HOST_DEVICE
0041   static size_t GetSafetyCandidates_v(Vector3D<Precision> const &point, ABBoxManager::ABBoxContainer_v const &corners,
0042                                       size_t size, ABBoxManager::BoxIdDistancePair_t *boxsafetypairs,
0043                                       Precision upper_squared_limit)
0044   {
0045     size_t count = 0;
0046     Vector3D<float> pointfloat((float)point.x(), (float)point.y(), (float)point.z());
0047     size_t vecsize = size;
0048     for (size_t box = 0; box < vecsize; ++box) {
0049       ABBoxManager::Float_v safetytoboxsqr =
0050           ABBoxImplementation::ABBoxSafetySqr(corners[2 * box], corners[2 * box + 1], pointfloat);
0051
0052       auto hit           = safetytoboxsqr < ABBoxManager::Real_t(upper_squared_limit);
0053       constexpr auto kVS = vecCore::VectorSize<ABBoxManager::Float_v>();
0054       if (!vecCore::MaskEmpty(hit)) {
0055         for (size_t i = 0; i < kVS; ++i) {
0056           if (vecCore::MaskLaneAt(hit, i)) {
0057             boxsafetypairs[count] =
0058                 ABBoxManager::BoxIdDistancePair_t(box * kVS + i, vecCore::LaneAt(safetytoboxsqr, i));
0059             count++;
0060           }
0061         }
0062       }
0063     }
0064     return count;
0065   }
0066
0067   VECGEOM_FORCE_INLINE
0068   VECCORE_ATT_HOST_DEVICE
0069   Precision TreatSafetyToIn(Vector3D<Precision> const &localpoint, LogicalVolume const *lvol, Precision outsafety) const
0070   {
0071     // a stack based workspace array
0072     // The following construct reserves stackspace for objects
0073     // of type IdDistPair_t WITHOUT initializing those objects
0074     using IdDistPair_t = ABBoxManager::BoxIdDistancePair_t;
0075     char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0076     IdDistPair_t *boxsafetylist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0077
0078     double safety    = outsafety; // we use the outsafety estimate as starting point
0079     double safetysqr = safety * safety;
0080
0081     // safety to bounding boxes
0082     if (safety > 0. && lvol->GetDaughtersp()->size() > 0) {
0083       int size;
0084
0085       ABBoxManager::ABBoxContainer_v bboxes = fABBoxManager.GetABBoxes_v(lvol, size);
0086       // calculate squared bounding box safeties in vectorized way
0087       auto ncandidates = GetSafetyCandidates_v(localpoint, bboxes, size, boxsafetylist, safetysqr);
0088       // not sorting the candidate list ( which one could do )
0089       for (unsigned int candidate = 0; candidate < ncandidates; ++candidate) {
0090         auto boxsafetypair = boxsafetylist[candidate];
0091         if (boxsafetypair.second < safetysqr) {
0092           VPlacedVolume const *cand = LookupDaughter(lvol, boxsafetypair.first);
0093           if (boxsafetypair.first > lvol->GetDaughtersp()->size()) break;
0094           auto candidatesafety = cand->SafetyToIn(localpoint);
0095 #ifdef VERBOSE
0096           if (candidatesafety * candidatesafety > boxsafetypair.second && boxsafetypair.second > 0)
0097             std::cerr << "real safety smaller than boxsafety \n";
0098 #endif
0099           if (candidatesafety < safety) {
0100             safety    = candidatesafety;
0101             safetysqr = safety * safety;
0102           }
0103         }
0104       }
0105     }
0106     return safety;
0107   }
0108
0109 public:
0110   static constexpr const char *gClassNameString = "SimpleABBoxSafetyEstimator";
0111
0112   VECGEOM_FORCE_INLINE
0113   VECCORE_ATT_HOST_DEVICE
0114   virtual Precision ComputeSafetyForLocalPoint(Vector3D<Precision> const &localpoint,
0115                                                VPlacedVolume const *pvol) const override
0116   {
0117
0118     // safety to mother
0119     double safety = pvol->SafetyToOut(localpoint);
0120     return TreatSafetyToIn(localpoint, pvol->GetLogicalVolume(), safety);
0121   }
0122
0123   // estimate just the safety to daughters for a local point with respect to a logical volume
0124   VECCORE_ATT_HOST_DEVICE
0125   virtual Precision ComputeSafetyToDaughtersForLocalPoint(Vector3D<Precision> const &localpoint,
0126                                                           LogicalVolume const *lvol) const override
0127   {
0128     return TreatSafetyToIn(localpoint, lvol, kInfLength);
0129   }
0130
0131   VECGEOM_FORCE_INLINE
0132   VECCORE_ATT_HOST_DEVICE
0133   virtual Real_v ComputeSafetyForLocalPoint(Vector3D<Real_v> const &localpoint, VPlacedVolume const *pvol,
0134                                             Bool_v m) const override
0135   {
0136     using vecCore::AssignLane;
0137     using vecCore::LaneAt;
0138     Real_v safety(0.);
0139     if (!vecCore::MaskEmpty(m)) {
0140       // SIMD safety to mother
0141       safety                    = pvol->SafetyToOut(localpoint);
0142       LogicalVolume const *lvol = pvol->GetLogicalVolume();
0143       // now loop over the voxelized treatment of safety to in
0144       for (unsigned int i = 0; i < vecCore::VectorSize<Real_v>(); ++i) {
0145         if (vecCore::MaskLaneAt(m, i)) {
0146           vecCore::AssignLane(safety, i,
0147                               TreatSafetyToIn(Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i),
0148                                                                   LaneAt(localpoint.z(), i)),
0149                                               lvol, vecCore::LaneAt(safety, i)));
0150         } else {
0151           vecCore::AssignLane(safety, i, 0.);
0152         }
0153       }
0154     }
0155     return safety;
0156   }
0157
0158   // vector interface
0159   VECGEOM_FORCE_INLINE
0160   virtual void ComputeSafetyForLocalPoints(SOA3D<Precision> const &localpoints, VPlacedVolume const *pvol,
0161                                            Precision *safeties) const override
0162   {
0163     // a stack based workspace array
0164     // The following construct reserves stackspace for objects
0165     // of type IdDistPair_t WITHOUT initializing those objects
0166     using IdDistPair_t = ABBoxManager::BoxIdDistancePair_t;
0167     char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0168     IdDistPair_t *boxsafetylist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0169
0170     // safety to mother -- using vector interface
0171     pvol->SafetyToOut(localpoints, safeties);
0172
0173     // safety to bounding boxes
0174     LogicalVolume const *lvol = pvol->GetLogicalVolume();
0175     if (!(lvol->GetDaughtersp()->size() > 0)) return;
0176
0177     // get bounding boxes (they are the same for all tracks)
0178     int numberofboxes;
0179     auto bboxes = fABBoxManager.GetABBoxes_v(lvol, numberofboxes);
0180
0181     // now loop over particles
0182     for (int i = 0, ntracks = localpoints.size(); i < ntracks; ++i) {
0183       double safety = safeties[i];
0184       if (safeties[i] > 0.) {
0185         double safetysqr = safeties[i] * safeties[i];
0186         auto lpoint      = localpoints[i];
0187         // vectorized search through bounding boxes -- quickly excluding many candidates
0188         auto ncandidates = GetSafetyCandidates_v(lpoint, bboxes, numberofboxes, boxsafetylist, safetysqr);
0189         // loop over remaining candidates
0190         for (unsigned int candidate = 0; candidate < ncandidates; ++candidate) {
0191           auto boxsafetypair = boxsafetylist[candidate];
0192           if (boxsafetypair.second < safetysqr) {
0193             VPlacedVolume const *cand = LookupDaughter(lvol, boxsafetypair.first);
0194             if (boxsafetypair.first > lvol->GetDaughtersp()->size()) break;
0195             auto candidatesafety = cand->SafetyToIn(lpoint);
0196 #ifdef VERBOSE
0197             if (candidatesafety * candidatesafety > boxsafetypair.second && boxsafetypair.second > 0)
0198               std::cerr << "real safety smaller than boxsafety \n";
0199 #endif
0200             if (candidatesafety < safety) {
0201               safety    = candidatesafety;
0202               safetysqr = safety * safety;
0203             }
0204           }
0205         }
0206       }
0207       // write back result
0208       safeties[i] = safety;
0209     }
0210   }
0211
0212   static VSafetyEstimator *Instance()
0213   {
0214     static SimpleABBoxSafetyEstimator instance;
0215     return &instance;
0216   }
0217
0218 }; // end class
0219 } // namespace VECGEOM_IMPL_NAMESPACE
0220 } // namespace vecgeom
0221
0222 #endif /* NAVIGATION_SIMPLEABBOXSAFETYESTIMATOR_H_ */