VecGeom/navigation/SimpleSafetyEstimator.h

0001 /*
0002  * VSimpleSafetyEstimator.h
0003  *
0004  *  Created on: 28.08.2015
0005  *      Author: swenzel
0006  */
0007
0008 #ifndef NAVIGATION_SIMPLESAFETYESTIMATOR_H_
0009 #define NAVIGATION_SIMPLESAFETYESTIMATOR_H_
0010
0011 #include "VecGeom/navigation/VSafetyEstimator.h"
0012 //#include "VecGeom/base/Array.h"
0013
0014 namespace vecgeom {
0015 inline namespace VECGEOM_IMPL_NAMESPACE {
0016
0017 /// Keep in dest the minimum between dest,temp for each corresponding element
0018 static void VectMin(unsigned int nelem, Precision *__restrict__ dest, const Precision *__restrict__ temp) {
0019   using vecCore::FromPtr;
0020   using Real_v = vecgeom::VectorBackend::Real_v;
0021
0022   // loop over all elements
0023   unsigned int i = 0;
0024   constexpr unsigned int nlanes = vecCore::VectorSize<Real_v>();
0025   const auto ilast = nelem - (nlanes-1);
0026   Real_v result;
0027   for ( ; i < ilast; i += nlanes) {
0028     result = vecCore::math::Min(FromPtr<Real_v>(dest + i), FromPtr<Real_v>(temp + i));
0029     vecCore::Store(result, &dest[i]);
0030   }
0031
0032   // fall back to scalar interface for tail treatment
0033   for (; i < nelem; ++i) {
0034     if ( dest[i] > temp[i] ) dest[i] = temp[i];
0035   }
0036 }
0037
0038 //! a simple safety estimator based on a brute force (O(N)) approach
0039 class SimpleSafetyEstimator : public VSafetyEstimatorHelper<SimpleSafetyEstimator> {
0040
0041 public:
0042   static constexpr const char *gClassNameString = "SimpleSafetyEstimator";
0043
0044   // estimate just the safety to daughters for a local point with respect to a logical volume
0045   // TODO: use this function in other interfaces to avoid code duplication
0046   VECCORE_ATT_HOST_DEVICE
0047   virtual Precision ComputeSafetyToDaughtersForLocalPoint(Vector3D<Precision> const &localpoint,
0048                                                           LogicalVolume const *lvol) const override
0049   {
0050     // safety to daughters
0051     double safety(kInfLength);
0052     auto daughters       = lvol->GetDaughtersp();
0053     auto numberdaughters = daughters->size();
0054     for (decltype(numberdaughters) d = 0; d < numberdaughters; ++d) {
0055       VPlacedVolume const *daughter = daughters->operator[](d);
0056       double tmp                    = daughter->SafetyToIn(localpoint);
0057       safety                        = Min(safety, tmp);
0058     }
0059     return safety;
0060   }
0061
0062   VECGEOM_FORCE_INLINE
0063   VECCORE_ATT_HOST_DEVICE
0064   virtual Precision ComputeSafetyForLocalPoint(Vector3D<Precision> const &localpoint,
0065                                                VPlacedVolume const *pvol) const override
0066   {
0067     // safety to mother
0068     double safety = pvol->SafetyToOut(localpoint);
0069
0070     // safety to daughters
0071     auto daughters       = pvol->GetLogicalVolume()->GetDaughtersp();
0072     auto numberdaughters = daughters->size();
0073     for (decltype(numberdaughters) d = 0; d < numberdaughters; ++d) {
0074       VPlacedVolume const *daughter = daughters->operator[](d);
0075       double tmp                    = daughter->SafetyToIn(localpoint);
0076       safety                        = Min(safety, tmp);
0077     }
0078     return safety;
0079   }
0080
0081   VECGEOM_FORCE_INLINE
0082   VECCORE_ATT_HOST_DEVICE
0083   virtual Real_v ComputeSafetyForLocalPoint(Vector3D<Real_v> const &localpoint, VPlacedVolume const *pvol,
0084                                             Bool_v m) const override
0085   {
0086     // safety to mother
0087     Real_v safety(0.);
0088     if (!vecCore::MaskEmpty(m)) {
0089       safety = pvol->SafetyToOut(localpoint);
0090
0091       // safety to daughters
0092       auto daughters       = pvol->GetLogicalVolume()->GetDaughtersp();
0093       auto numberdaughters = daughters->size();
0094       for (decltype(numberdaughters) d = 0; d < numberdaughters; ++d) {
0095         VPlacedVolume const *daughter = daughters->operator[](d);
0096         auto tmp                      = daughter->SafetyToIn(localpoint);
0097         safety                        = Min(safety, tmp);
0098       }
0099     }
0100     return safety;
0101   }
0102
0103
0104   VECGEOM_FORCE_INLINE
0105   virtual void ComputeSafetyForLocalPoints(SOA3D<Precision> const &localpoints, VPlacedVolume const *pvol,
0106                                            Precision *safeties) const override
0107   {
0108     auto npoints = localpoints.size();
0109     // a stack based workspace array
0110     // The following construct reserves stackspace WITHOUT initializing it
0111     char stackspace[npoints * sizeof(Precision)];
0112     Precision *tmpSafeties = reinterpret_cast<Precision *>(&stackspace);
0113
0114     // // Array-based - transparently ensures proper alignment
0115     // // NavigBenchmark performance is actually worse than above though
0116     // Array<Precision> stackspace(npoints);
0117     // Precision *tmpSafeties = reinterpret_cast<Precision *>(stackspace.begin());
0118
0119     pvol->SafetyToOut(localpoints, safeties);
0120
0121     // safety to daughters; brute force but each function (possibly) vectorized
0122     Vector<Daughter> const *daughters = pvol->GetLogicalVolume()->GetDaughtersp();
0123     auto numberdaughters              = daughters->size();
0124     for (decltype(numberdaughters) d = 0; d < numberdaughters; ++d) {
0125       VPlacedVolume const *daughter = daughters->operator[](d);
0126       daughter->SafetyToIn(localpoints, tmpSafeties);
0127
0128       //.. keep track of closest daughter
0129       VectMin(npoints, safeties, tmpSafeties);
0130     }
0131     // here the safeties array automatically contains the right values
0132   }
0133
0134   // bring in interface from higher up
0135   using VSafetyEstimatorHelper<SimpleSafetyEstimator>::ComputeVectorSafety;
0136
0137   // implementation (using VC that does not need a workspace) by doing the whole algorithm in little vector chunks
0138   virtual void ComputeVectorSafety(SOA3D<Precision> const &globalpoints, NavStatePool &states,
0139                                    Precision *safeties) const override
0140   {
0141     VPlacedVolume const *pvol = states[0]->Top();
0142     auto npoints              = globalpoints.size();
0143
0144     constexpr auto kVS = vecCore::VectorSize<Real_v>();
0145     for (decltype(npoints) i = 0; i < npoints; i += kVS) {
0146       // do the transformation to local
0147       Vector3D<Real_v> local;
0148       for (size_t j = 0; j < kVS; ++j) {
0149         Transformation3D m;
0150         states[i + j]->TopMatrix(m);
0151         auto v = m.Transform(globalpoints[i + j]);
0152
0153         using vecCore::AssignLane;
0154         AssignLane(local.x(), j, v.x());
0155         AssignLane(local.y(), j, v.y());
0156         AssignLane(local.z(), j, v.z());
0157       }
0158
0159       auto safety          = pvol->SafetyToOut(local);
0160       auto daughters       = pvol->GetLogicalVolume()->GetDaughtersp();
0161       auto numberdaughters = daughters->size();
0162       for (decltype(numberdaughters) d = 0; d < numberdaughters; ++d) {
0163         VPlacedVolume const *daughter = daughters->operator[](d);
0164         safety                        = Min(safety, daughter->SafetyToIn(local));
0165       }
0166       vecCore::Store(safety, safeties + i);
0167     }
0168   }
0169
0170 #ifndef VECCORE_CUDA
0171   static VSafetyEstimator *Instance()
0172   {
0173     static SimpleSafetyEstimator instance;
0174     return &instance;
0175   }
0176 #else
0177   VECCORE_ATT_DEVICE
0178   static VSafetyEstimator *Instance();
0179 #endif
0180
0181 }; // end class
0182 }
0183 } // end namespaces
0184
0185 #endif /* NAVIGATION_SIMPLESAFETYESTIMATOR_H_ */