VecGeom/navigation/VoxelSafetyEstimator.h

0001 // This file is part of VecGeom and is distributed under the
0002 // conditions in the file LICENSE.txt in the top directory.
0003 // For the full list of authors see CONTRIBUTORS.txt and `git log`.
0004
0005 /// \file navigation/VoxelSafetyEstimator.h
0006 /// \author Sandro Wenzel (CERN)
0007
0008 #ifndef NAVIGATION_VOXELSAFETYESTIMATOR_H_
0009 #define NAVIGATION_VOXELSAFETYESTIMATOR_H_
0010
0011 #include "VecGeom/navigation/VSafetyEstimator.h"
0012 #include "VecGeom/management/FlatVoxelManager.h"
0013 #include <cmath>
0014
0015 namespace vecgeom {
0016 inline namespace VECGEOM_IMPL_NAMESPACE {
0017
0018 //! a safety estimator using a voxel lookup table of precomputed safety candidates
0019 //! (or in the extreme case of precomputed safety values)
0020 class VoxelSafetyEstimator : public VSafetyEstimatorHelper<VoxelSafetyEstimator> {
0021
0022 private:
0023   const FlatVoxelManager &fAccStructureManager;
0024
0025   VoxelSafetyEstimator()
0026       : VSafetyEstimatorHelper<VoxelSafetyEstimator>(), fAccStructureManager{FlatVoxelManager::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   static constexpr const char *gClassNameString = "VoxelSafetyEstimator";
0040
0041   VECGEOM_FORCE_INLINE
0042   VECCORE_ATT_HOST_DEVICE
0043   Precision TreatSafetyToIn(Vector3D<Precision> const &localpoint, LogicalVolume const *lvol, Precision outsafety) const
0044   {
0045     throw std::runtime_error("unimplemented function called");
0046   }
0047
0048   // #define CROSSCHECK
0049
0050 #ifdef CROSSCHECK
0051   void printProcedure(int const *safetycandidates, int length, LogicalVolume const *lvol,
0052                       Vector3D<Precision> const &localpoint) const
0053   {
0054     int size{0};
0055     auto abboxcorners = ABBoxManager::Instance().GetABBoxes(lvol, size);
0056     Vector3D<float> lp(localpoint.x(), localpoint.y(), localpoint.z()); // in float
0057     const bool needmother        = (safetycandidates[0] == -1);
0058     const Precision safetymother = needmother ? lvol->GetUnplacedVolume()->SafetyToOut(localpoint) : 0.;
0059     std::cerr << "safetymother " << safetymother;
0060     Precision safetyToDaughters{1E20};
0061     // now talk to daughter candidates - with bounding box speedup
0062     // TODO: this could actually happen with SIMD speedup
0063     float bestSafetySqr{1E20};
0064     for (int i = 0; i < length; ++i) {
0065       if (safetycandidates[i] == -1) continue;
0066       const auto candidateid = safetycandidates[i];
0067       const auto &lower      = abboxcorners[2 * candidateid];
0068       const auto &upper      = abboxcorners[2 * candidateid + 1];
0069       const float ssqr       = ABBoxImplementation::ABBoxSafetySqr(lower, upper, lp);
0070       std::cerr << i << " boxsqr " << ssqr << "\t" << std::sqrt(ssqr) << "\n";
0071       if (ssqr > 0) {
0072         bestSafetySqr = std::min(bestSafetySqr, ssqr);
0073       } else { // inside
0074         const auto daughter = LookupDaughter(lvol, safetycandidates[i]);
0075         auto sin            = daughter->SafetyToIn(localpoint);
0076         std::cerr << i << " SI " << sin << "\n";
0077         safetyToDaughters = std::min(safetyToDaughters, daughter->SafetyToIn(localpoint));
0078       }
0079     }
0080     safetyToDaughters = std::min(safetyToDaughters, (Precision)std::sqrt(bestSafetySqr));
0081     std::cerr << "final safetyd " << safetyToDaughters << "\n";
0082     const float returnvalue = needmother ? std::min(safetyToDaughters, safetymother) : safetyToDaughters;
0083     std::cerr << "returning " << returnvalue << "\n";
0084   }
0085 #endif
0086
0087   VECGEOM_FORCE_INLINE
0088   VECCORE_ATT_HOST_DEVICE
0089   virtual Precision ComputeSafetyForLocalPoint(Vector3D<Precision> const &localpoint,
0090                                                VPlacedVolume const *pvol) const override
0091   {
0092     static int counter = 0;
0093     counter++;
0094     const auto lvol      = pvol->GetLogicalVolume();
0095     const auto structure = fAccStructureManager.GetStructure(lvol);
0096
0097     int size{0};
0098     auto abboxcorners = ABBoxManager::Instance().GetABBoxes(lvol, size);
0099
0100     if (structure != nullptr) {
0101       const Vector3D<float> lp(localpoint.x(), localpoint.y(), localpoint.z()); // in float
0102       // fetch safety properties for localpoint;
0103       const int *safetycandidates{nullptr};
0104
0105       int length{0};
0106       auto voxelHashMap = structure->fVoxelToCandidate;
0107       if (voxelHashMap == nullptr) {
0108         std::cerr << "ERROR> VoxelSafetyEstimator::ComputeSafetyForLocalPoint call# " << counter
0109                   << " no structure VoxelToCandidate structure found for volume "
0110                   << " phys: " << pvol->GetName() << "  physvol id = " << pvol->id() << " log : " << lvol->GetName()
0111                   << "  log vol id = " << lvol->id() << std::endl;
0112         return 0.0;
0113       }
0114       safetycandidates = voxelHashMap->getProperties(lp, length);
0115       if (length > 0) {
0116         const bool needmother        = true; //(safetycandidates[0] == -1);
0117         const Precision safetymother = needmother ? lvol->GetUnplacedVolume()->SafetyToOut(localpoint) : 0.;
0118         Precision safetyToDaughters{1E20};
0119         // now talk to daughter candidates - with bounding box speedup
0120         // TODO: this could actually happen with SIMD speedup
0121         float bestSafetySqr{1E20};
0122         for (int i = 0; i < length; ++i) {
0123           if (safetycandidates[i] == -1) continue;
0124           const auto candidateid = safetycandidates[i];
0125           const auto &lower      = abboxcorners[2 * candidateid];
0126           const auto &upper      = abboxcorners[2 * candidateid + 1];
0127           const float ssqr       = ABBoxImplementation::ABBoxSafetySqr(lower, upper, lp);
0128           if (ssqr > 0) {
0129             bestSafetySqr = std::min(bestSafetySqr, ssqr);
0130           } else { // inside
0131             const auto daughter = LookupDaughter(lvol, safetycandidates[i]);
0132             auto s              = daughter->SafetyToIn(localpoint);
0133             if (s <= 0.) {
0134               // can break here
0135               s = 0.;
0136             }
0137             safetyToDaughters = std::min(safetyToDaughters, s);
0138           }
0139         }
0140         safetyToDaughters       = std::min(safetyToDaughters, (Precision)std::sqrt(bestSafetySqr));
0141         const float returnvalue = needmother ? std::min(safetyToDaughters, safetymother) : safetyToDaughters;
0142
0143 #ifdef CROSSCHECK
0144         if (returnvalue > lvol->GetUnplacedVolume()->SafetyToOut(localpoint) + 1E-4) {
0145           std::cerr << returnvalue << " " << lvol->GetUnplacedVolume()->SafetyToOut(localpoint)
0146                     << "Problem in voxel safety for point " << lp << " voxelkey "
0147                     << structure->fVoxelToCandidate->getKey(lp.x(), lp.y(), lp.z()) << " ";
0148           std::cerr << "{ ";
0149           for (int i = 0; i < length; ++i) {
0150             std::cout << safetycandidates[i] << " , ";
0151           }
0152           std::cerr << " }\n";
0153           printProcedure(safetycandidates, length, lvol, localpoint);
0154         }
0155 #endif
0156
0157         if (returnvalue < -1.0e-10) {
0158           std::cerr << "VoxelSafetyEstimator: returning negative value.  saf-to-daughters= " << safetyToDaughters
0159                     << " saf-to-mother = " << safetymother << "\n";
0160         }
0161
0162         return returnvalue;
0163       } else {
0164         // no information for this voxel present
0165         std::cerr << "Warning> ComputeSafetyForLocalPoint call# " << counter
0166                   << " no information for this voxel present " << localpoint << " at key "
0167                   << structure->fVoxelToCandidate->getKey(lp.x(), lp.y(), lp.z()) << " \n ";
0168         return 0.;
0169       }
0170     }
0171     std::cerr << " No voxel information found; Cannot use voxel safety\n";
0172     return 0.;
0173   }
0174
0175   VECGEOM_FORCE_INLINE
0176   VECCORE_ATT_HOST_DEVICE
0177   virtual Precision ComputeSafetyToDaughtersForLocalPoint(Vector3D<Precision> const &localpoint,
0178                                                           LogicalVolume const *lvol) const override
0179   {
0180     return TreatSafetyToIn(localpoint, lvol, kInfLength);
0181   }
0182
0183   //
0184   // These interfaces for baskets are only dummy-implemented at the moment
0185   //
0186   VECCORE_ATT_HOST_DEVICE
0187   Real_v ComputeSafetyForLocalPoint(Vector3D<Real_v> const & /*localpoint*/, VPlacedVolume const * /*pvol*/,
0188                                     Bool_v /*m*/) const override
0189   {
0190     throw std::runtime_error("unimplemented function called");
0191     return Real_v(0.);
0192   }
0193
0194   // interfaces to treat vectors/collections of points (uses the approach with intermediate storage and passing down the
0195   // loops to shapes)
0196   void ComputeVectorSafety(SOA3D<Precision> const & /*globalpoints*/, NavStatePool &states,
0197                            SOA3D<Precision> & /*workspace*/, Precision * /*safeties*/) const override
0198   {
0199     throw std::runtime_error("unimplemented function called");
0200   }
0201
0202   // interfaces to treat vectors/collections of points (uses the approach without intermediate storage; requires access
0203   // to new SIMD interface)
0204   void ComputeVectorSafety(SOA3D<Precision> const & /*globalpoints*/, NavStatePool & /*states*/,
0205                            Precision * /*safeties*/) const override
0206   {
0207     throw std::runtime_error("unimplemented function called");
0208   }
0209
0210   void ComputeSafetyForLocalPoints(SOA3D<Precision> const & /*localpoints*/, VPlacedVolume const * /*pvol*/,
0211                                    Precision * /*safeties*/) const override
0212   {
0213     throw std::runtime_error("unimplemented function called");
0214   }
0215
0216   static VSafetyEstimator *Instance()
0217   {
0218     static VoxelSafetyEstimator instance;
0219     return &instance;
0220   }
0221
0222 }; // end class
0223 } // namespace VECGEOM_IMPL_NAMESPACE
0224 } // namespace vecgeom
0225
0226 #endif /* NAVIGATION_VOXELSAFETYESTIMATOR_H_ */