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<Precision> &fABBoxManager;
0024
0025   SimpleABBoxSafetyEstimator()
0026       : VSafetyEstimatorHelper<SimpleABBoxSafetyEstimator>(), 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 public:
0039   // helper function calculating some candidate volumes
0040   VECCORE_ATT_HOST_DEVICE
0041   static size_t GetSafetyCandidates_v(Vector3D<Precision> const &point,
0042                                       ABBoxManager<Precision>::ABBoxContainer_v const &corners, size_t size,
0043                                       ABBoxManager<Precision>::BoxIdDistancePair_t *boxsafetypairs,
0044                                       Precision upper_squared_limit)
0045   {
0046     size_t count = 0;
0047     Vector3D<float> pointfloat((float)point.x(), (float)point.y(), (float)point.z());
0048     size_t vecsize = size;
0049     for (size_t box = 0; box < vecsize; ++box) {
0050       ABBoxManager<Precision>::Float_v safetytoboxsqr =
0051           ABBoxImplementation::ABBoxSafetySqr(corners[2 * box], corners[2 * box + 1], pointfloat);
0052
0053       auto hit           = safetytoboxsqr < ABBoxManager<Precision>::Real_s(upper_squared_limit);
0054       constexpr auto kVS = vecCore::VectorSize<ABBoxManager<Precision>::Float_v>();
0055       if (!vecCore::MaskEmpty(hit)) {
0056         for (size_t i = 0; i < kVS; ++i) {
0057           if (vecCore::MaskLaneAt(hit, i)) {
0058             boxsafetypairs[count] =
0059                 ABBoxManager<Precision>::BoxIdDistancePair_t(box * kVS + i, vecCore::LaneAt(safetytoboxsqr, i));
0060             count++;
0061           }
0062         }
0063       }
0064     }
0065     return count;
0066   }
0067
0068   VECGEOM_FORCE_INLINE
0069   VECCORE_ATT_HOST_DEVICE
0070   Precision TreatSafetyToIn(Vector3D<Precision> const &localpoint, LogicalVolume const *lvol, Precision outsafety) const
0071   {
0072     // a stack based workspace array
0073     // The following construct reserves stackspace for objects
0074     // of type IdDistPair_t WITHOUT initializing those objects
0075     using IdDistPair_t = ABBoxManager<Precision>::BoxIdDistancePair_t;
0076     char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0077     IdDistPair_t *boxsafetylist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0078
0079     double safety    = outsafety; // we use the outsafety estimate as starting point
0080     double safetysqr = safety * safety;
0081
0082     // safety to bounding boxes
0083     if (safety > 0. && lvol->GetDaughtersp()->size() > 0) {
0084       int size;
0085
0086       ABBoxManager<Precision>::ABBoxContainer_v bboxes = fABBoxManager.GetABBoxes_v(lvol, size);
0087       // calculate squared bounding box safeties in vectorized way
0088       auto ncandidates = GetSafetyCandidates_v(localpoint, bboxes, size, boxsafetylist, safetysqr);
0089       // not sorting the candidate list ( which one could do )
0090       for (unsigned int candidate = 0; candidate < ncandidates; ++candidate) {
0091         auto boxsafetypair = boxsafetylist[candidate];
0092         if (boxsafetypair.second < safetysqr) {
0093           VPlacedVolume const *cand = LookupDaughter(lvol, boxsafetypair.first);
0094           if (boxsafetypair.first > lvol->GetDaughtersp()->size()) break;
0095           auto candidatesafety = cand->SafetyToIn(localpoint);
0096 #ifdef VERBOSE
0097           if (candidatesafety * candidatesafety > boxsafetypair.second && boxsafetypair.second > 0)
0098             std::cerr << "real safety smaller than boxsafety \n";
0099 #endif
0100           if (candidatesafety < safety) {
0101             safety    = candidatesafety;
0102             safetysqr = safety * safety;
0103           }
0104         }
0105       }
0106     }
0107     return safety;
0108   }
0109
0110 public:
0111   static constexpr const char *gClassNameString = "SimpleABBoxSafetyEstimator";
0112
0113   VECGEOM_FORCE_INLINE
0114   VECCORE_ATT_HOST_DEVICE
0115   virtual Precision ComputeSafetyForLocalPoint(Vector3D<Precision> const &localpoint,
0116                                                VPlacedVolume const *pvol) const override
0117   {
0118
0119     // safety to mother
0120     double safety = pvol->SafetyToOut(localpoint);
0121     return TreatSafetyToIn(localpoint, pvol->GetLogicalVolume(), safety);
0122   }
0123
0124   // estimate just the safety to daughters for a local point with respect to a logical volume
0125   VECCORE_ATT_HOST_DEVICE
0126   virtual Precision ComputeSafetyToDaughtersForLocalPoint(Vector3D<Precision> const &localpoint,
0127                                                           LogicalVolume const *lvol) const override
0128   {
0129     return TreatSafetyToIn(localpoint, lvol, kInfLength);
0130   }
0131
0132   static VSafetyEstimator *Instance()
0133   {
0134     static SimpleABBoxSafetyEstimator instance;
0135     return &instance;
0136   }
0137
0138 }; // end class
0139 } // namespace VECGEOM_IMPL_NAMESPACE
0140 } // namespace vecgeom
0141
0142 #endif /* NAVIGATION_SIMPLEABBOXSAFETYESTIMATOR_H_ */