[integration] Refactor hit scoring and AdePTGeant4Integration.

- Refactor the processing of hits.
  Instead of processing hits by passing a pointer/reference to a
  HostScoring instance, a loop over iterators to hits is used.
  In this way, hit scoring is decoupled from the specific implementation
  of HostScoring, and all classes with the same interface as the original
  GPUHit can be used for scoring.
  This facilitates hit scoring for the AsyncTransport implementation.
- Move several Geant4 objects into the .cpp to make the integration headers
  simpler.
- Place temporary scoring objects into a struct to go around G4's pool
  allocators. This prevents a destruction order fiasco (where the pool
  is gone but the object isn't), and keeps the scoring objects closer
  in memory.
  A few objects need to leak, unfortunately, since they are allocated in
  G4 pools, and the handles don't support them being on a stack.
- Improve const correctness in a few places.
- Add information about threadID and eventID to the scoring interface.
  This information is required for AsyncAdePT to score correctly, but is
  unused in the thread-local transport for now.

include/AdePT/core/AdePTScoringTemplate.cuh

@@ -18,25 +18,30 @@ namespace adept_scoring
   void FreeGPU(Scoring *scoring, Scoring *scoring_dev){}
 
   template <typename Scoring>
-  __device__ void RecordHit(Scoring *scoring_dev, int aParentID, char aParticleType, double aStepLength, double aTotalEnergyDeposit,
-                          vecgeom::NavigationState const *aPreState, vecgeom::Vector3D<Precision> *aPrePosition,
-                          vecgeom::Vector3D<Precision> *aPreMomentumDirection,
-                          vecgeom::Vector3D<Precision> *aPrePolarization, double aPreEKin, double aPreCharge,
-                          vecgeom::NavigationState const *aPostState, vecgeom::Vector3D<Precision> *aPostPosition,
-                          vecgeom::Vector3D<Precision> *aPostMomentumDirection,
-                          vecgeom::Vector3D<Precision> *aPostPolarization, double aPostEKin, double aPostCharge){}
+  __device__ void RecordHit(Scoring *scoring_dev, int aParentID, char aParticleType, double aStepLength,
+                            double aTotalEnergyDeposit, vecgeom::NavigationState const *aPreState,
+                            vecgeom::Vector3D<Precision> const *aPrePosition,
+                            vecgeom::Vector3D<Precision> const *aPreMomentumDirection,
+                            vecgeom::Vector3D<Precision> const *aPrePolarization, double aPreEKin, double aPreCharge,
+                            vecgeom::NavigationState const *aPostState,
+                            vecgeom::Vector3D<Precision> const *aPostPosition,
+                            vecgeom::Vector3D<Precision> const *aPostMomentumDirection,
+                            vecgeom::Vector3D<Precision> const *aPostPolarization, double aPostEKin, double aPostCharge,
+                            unsigned int eventId, short threadId);
 
-template <typename Scoring>
-__device__ void AccountProduced(Scoring *scoring_dev, int num_ele, int num_pos, int num_gam);
+  template <typename Scoring>
+  __device__ void AccountProduced(Scoring *scoring_dev, int num_ele, int num_pos, int num_gam);
 
-template <typename Scoring>
-__device__ __forceinline__ void EndOfIterationGPU(Scoring *scoring_dev);
+  template <typename Scoring>
+  __device__ __forceinline__ void EndOfIterationGPU(Scoring *scoring_dev);
 
-template <typename Scoring, typename IntegrationLayer>
-inline void EndOfIteration(Scoring &scoring, Scoring *scoring_dev, cudaStream_t &stream, IntegrationLayer &integration);
+  template <typename Scoring, typename IntegrationLayer>
+  inline void EndOfIteration(Scoring &scoring, Scoring *scoring_dev, cudaStream_t &stream,
+                             IntegrationLayer &integration);
 
-template <typename Scoring, typename IntegrationLayer>
-inline void EndOfTransport(Scoring &scoring, Scoring *scoring_dev, cudaStream_t &stream, IntegrationLayer &integration);
+  template <typename Scoring, typename IntegrationLayer>
+  inline void EndOfTransport(Scoring &scoring, Scoring *scoring_dev, cudaStream_t &stream,
+                             IntegrationLayer &integration);
 }
 
 #endif

include/AdePT/core/AdePTTransport.icc

@@ -249,7 +249,7 @@ void AdePTTransport<IntegrationLayer>::Shower(int event)
     std::cout << "[" << tid << "] fromDevice: " << nelec << " elec, " << nposi << " posi, " << ngamma << " gamma\n";
   }
 
-  fIntegrationLayer.ReturnTracks(&(fBuffer.fromDevice), fDebugLevel);
+  fIntegrationLayer.ReturnTracks(fBuffer.fromDevice.begin(), fBuffer.fromDevice.end(), fDebugLevel);
 
   fBuffer.Clear();
 }

include/AdePT/core/HostScoringImpl.cuh

@@ -46,8 +46,8 @@ unsigned int aHitIndex = GetNextFreeHitIndex(hostScoring_dev);
 }
 
 /// @brief Utility function to copy a 3D vector, used for filling the Step Points
-__device__ __forceinline__ void Copy3DVector(vecgeom::Vector3D<Precision> *source,
-                                            vecgeom::Vector3D<Precision> *destination)
+__device__ __forceinline__ void Copy3DVector(vecgeom::Vector3D<Precision> const *source,
+                                             vecgeom::Vector3D<Precision> *destination)
 {
   destination->x() = source->x();
   destination->y() = source->y();
@@ -151,14 +151,13 @@ namespace adept_scoring
 
   /// @brief Record a hit
   template <>
-  __device__ void RecordHit(HostScoring *hostScoring_dev, int aParentID, char aParticleType, double aStepLength,
-                          double aTotalEnergyDeposit, vecgeom::NavigationState const *aPreState,
-                          vecgeom::Vector3D<Precision> *aPrePosition,
-                          vecgeom::Vector3D<Precision> *aPreMomentumDirection,
-                          vecgeom::Vector3D<Precision> *aPrePolarization, double aPreEKin, double aPreCharge,
-                          vecgeom::NavigationState const *aPostState, vecgeom::Vector3D<Precision> *aPostPosition,
-                          vecgeom::Vector3D<Precision> *aPostMomentumDirection,
-                          vecgeom::Vector3D<Precision> *aPostPolarization, double aPostEKin, double aPostCharge)
+  __device__ void RecordHit(
+      HostScoring *hostScoring_dev, int aParentID, char aParticleType, double aStepLength, double aTotalEnergyDeposit,
+      vecgeom::NavigationState const *aPreState, vecgeom::Vector3D<Precision> const *aPrePosition,
+      vecgeom::Vector3D<Precision> const *aPreMomentumDirection, vecgeom::Vector3D<Precision> const *aPrePolarization,
+      double aPreEKin, double aPreCharge, vecgeom::NavigationState const *aPostState,
+      vecgeom::Vector3D<Precision> const *aPostPosition, vecgeom::Vector3D<Precision> const *aPostMomentumDirection,
+      vecgeom::Vector3D<Precision> const *aPostPolarization, double aPostEKin, double aPostCharge, unsigned int, short)
   {
     // Acquire a hit slot
     GPUHit *aGPUHit = GetNextFreeHit(hostScoring_dev);
@@ -217,7 +216,9 @@ namespace adept_scoring
       // Synchronize the stream used to copy back the hits
       COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream));
       // Process the hits on CPU
-      integration.ProcessGPUHits(hostScoring, hostScoring.fStats);
+      for (const auto &hit : hostScoring) {
+        integration.ProcessGPUHit(hit);
+      }
     }
   }
 
@@ -231,7 +232,9 @@ namespace adept_scoring
     CopyGlobalCountersToHost(hostScoring, stream);
     COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream));
     // Process the last hits on CPU
-    integration.ProcessGPUHits(hostScoring, hostScoring.fStats);
+    for (const auto &hit : hostScoring) {
+      integration.ProcessGPUHit(hit);
+    }
   }
 }
 

include/AdePT/core/HostScoringStruct.cuh

@@ -22,7 +22,6 @@ struct GPUStepPoint {
 // call the user-defined Geant4 sensitive detector code
 struct GPUHit {
   int fParentID{0}; // Track ID
-  char fParticleType{0}; // Particle type ID
   // Data needed to reconstruct G4 Step
   double fStepLength{0};
   double fTotalEnergyDeposit{0};
@@ -32,6 +31,9 @@ struct GPUHit {
   // Data needed to reconstruct pre-post step points
   GPUStepPoint fPreStepPoint;
   GPUStepPoint fPostStepPoint;
+  unsigned int fEventId{0};
+  short threadId{-1};
+  char fParticleType{0}; // Particle type ID
 };
 
 /// @brief Stores information used for comparison with Geant4 (Number of steps, Number of produced particles, etc)
@@ -65,6 +67,29 @@ struct HostScoring {
     unsigned int fBufferStart; ///< Index of first used hit slot in the buffer
   };
 
+  struct Iterator {
+    std::size_t counter;
+    std::size_t const modulus;
+    GPUHit *const storage;
+
+    GPUHit &operator*() { return storage[counter % modulus]; }
+    Iterator &operator++()
+    {
+      counter++;
+      return *this;
+    }
+    Iterator operator++(int)
+    {
+      Iterator result = *this;
+      counter++;
+      return result;
+    }
+    bool operator!=(Iterator const &other)
+    {
+      return counter != other.counter || modulus != other.modulus || storage != other.storage;
+    }
+  };
+
   HostScoring(unsigned int aBufferCapacity = 1024 * 1024, float aFlushLimit = 0.8)
       : fBufferCapacity(aBufferCapacity), fFlushLimit(aFlushLimit)
   {
@@ -81,6 +106,12 @@ struct HostScoring {
     free(fGlobalCounters_host);
   }
 
+  Iterator begin() const { return Iterator{fStats.fBufferStart, fBufferCapacity, fGPUHitsBuffer_host}; }
+  Iterator end() const
+  {
+    return Iterator{fStats.fBufferStart + fStats.fUsedSlots, fBufferCapacity, fGPUHitsBuffer_host};
+  }
+
   /// @brief Print scoring info
   void Print();
 

include/AdePT/integration/AdePTGeant4Integration.hh

@@ -9,23 +9,22 @@
 #ifndef ADEPTGEANT4_INTEGRATION_H
 #define ADEPTGEANT4_INTEGRATION_H
 
-#include <unordered_map>
-
-#include <G4HepEmState.hh>
-
 #include <AdePT/core/CommonStruct.h>
 #include <AdePT/core/HostScoringStruct.cuh>
 
-#include <G4VPhysicalVolume.hh>
-#include <G4LogicalVolume.hh>
-#include <G4VPhysicalVolume.hh>
-#include <G4NavigationHistory.hh>
-#include <G4Step.hh>
-#include <G4Event.hh>
+#include <G4HepEmState.hh>
+
 #include <G4EventManager.hh>
+#include <G4Event.hh>
 
-#include <VecGeom/volumes/PlacedVolume.h>
-#include <VecGeom/volumes/LogicalVolume.h>
+#include <unordered_map>
+
+namespace AdePTGeant4Integration_detail {
+struct ScoringObjects;
+struct Deleter {
+  void operator()(ScoringObjects *ptr);
+};
+} // namespace AdePTGeant4Integration_detail
 
 class AdePTGeant4Integration {
 public:
@@ -42,43 +41,47 @@ public:
 
   /// @brief Fills the auxiliary data needed for AdePT
   static void InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmState *hepEmState, bool trackInAllRegions,
-                             std::vector<std::string> *gpuRegionNames);
+                             std::vector<std::string> const *gpuRegionNames);
 
   /// @brief Initializes the mapping of VecGeom to G4 volumes for sensitive volumes and their parents
   void InitScoringData(adeptint::VolAuxData *volAuxData);
 
   /// @brief Reconstructs GPU hits on host and calls the user-defined sensitive detector code
-  void ProcessGPUHits(HostScoring &aScoring, HostScoring::Stats &aStats);
-
-  /// @brief Takes a buffer of tracks coming from the device and gives them back to Geant4
-  void ReturnTracks(std::vector<adeptint::TrackData> *tracksFromDevice, int debugLevel);
+  void ProcessGPUHit(GPUHit const &hit);
+
+  /// @brief Takes a range of tracks coming from the device and gives them back to Geant4
+  template <typename Iterator>
+  void ReturnTracks(Iterator begin, Iterator end, int debugLevel) const
+  {
+    if (debugLevel > 1) {
+      G4cout << "Returning " << end - begin << " tracks from device" << G4endl;
+    }
+    for (Iterator it = begin; it != end; ++it) {
+      ReturnTrack(*it, it - begin, debugLevel);
+    }
+  }
 
   /// @brief Returns the Z value of the user-defined uniform magnetic field
   /// @details This function can only be called when the user-defined field is a G4UniformMagField
-  double GetUniformFieldZ();
+  double GetUniformFieldZ() const;
 
-  int GetEventID() { return G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID(); }
+  int GetEventID() const { return G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID(); }
 
-  int GetThreadID() { return G4Threading::G4GetThreadId(); }
+  int GetThreadID() const { return G4Threading::G4GetThreadId(); }
 
 private:
   /// @brief Reconstruct G4TouchableHistory from a VecGeom Navigation index
-  void FillG4NavigationHistory(vecgeom::NavigationState aNavState, G4NavigationHistory *aG4NavigationHistory);
-
-  void FillG4Step(GPUHit *aGPUHit, G4Step *aG4Step, G4TouchableHandle &aPreG4TouchableHandle,
-                  G4TouchableHandle &aPostG4TouchableHandle);
-
-  std::unordered_map<size_t, const G4VPhysicalVolume *> fglobal_vecgeom_to_g4_map; ///< Maps Vecgeom PV IDs to G4 PV IDs
-
-  bool fScoringObjectsInitialized{false};
-  G4NavigationHistory *fPreG4NavigationHistory{nullptr};
-  G4NavigationHistory *fPostG4NavigationHistory{nullptr};
-  G4Step *fG4Step{nullptr};
-  G4TouchableHandle fPreG4TouchableHistoryHandle;
-  G4TouchableHandle fPostG4TouchableHistoryHandle;
-  G4Track *fElectronTrack{nullptr};
-  G4Track *fPositronTrack{nullptr};
-  G4Track *fGammaTrack{nullptr};
+  void FillG4NavigationHistory(vecgeom::NavigationState aNavState, G4NavigationHistory *aG4NavigationHistory) const;
+
+  void FillG4Step(GPUHit const *aGPUHit, G4Step *aG4Step, G4TouchableHandle &aPreG4TouchableHandle,
+                  G4TouchableHandle &aPostG4TouchableHandle) const;
+
+  void ReturnTrack(adeptint::TrackData const &track, unsigned int trackIndex, int debugLevel) const;
+
+  std::unordered_map<size_t,
+                     const G4VPhysicalVolume *> fglobal_vecgeom_to_g4_map; ///< Maps Vecgeom PV IDs to G4 PV IDs
+  std::unique_ptr<AdePTGeant4Integration_detail::ScoringObjects, AdePTGeant4Integration_detail::Deleter>
+      fScoringObjects{nullptr};
 };
 
 #endif

include/AdePT/kernels/electrons.cuh

@@ -279,7 +279,8 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
                                &dir,                     // Post-step point momentum direction
                                nullptr,                  // Post-step point polarization
                                eKin,                     // Post-step point kinetic energy
-                               IsElectron ? -1 : 1);     // Post-step point charge
+                               IsElectron ? -1 : 1,      // Post-step point charge
+                               0, -1);                   // eventID and threadID (not needed here)
 
     // Save the `number-of-interaction-left` in our track.
     for (int ip = 0; ip < 3; ++ip) {

include/AdePT/kernels/gammas.cuh

@@ -259,7 +259,8 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
                                    &dir,                  // Post-step point momentum direction
                                    nullptr,               // Post-step point polarization
                                    newEnergyGamma,        // Post-step point kinetic energy
-                                   0);                    // Post-step point charge
+                                   0,                     // Post-step point charge
+                                   0, -1);                // event and thread ID
       }
 
       // Check the new gamma energy and deposit if below threshold.
@@ -285,7 +286,8 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
                                    &dir,                  // Post-step point momentum direction
                                    nullptr,               // Post-step point polarization
                                    newEnergyGamma,        // Post-step point kinetic energy
-                                   0);                    // Post-step point charge
+                                   0,                     // Post-step point charge
+                                   0, -1);                // event and thread ID
         // The current track is killed by not enqueuing into the next activeQueue.
       }
       break;
@@ -333,7 +335,8 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
                                  &dir,                  // Post-step point momentum direction
                                  nullptr,               // Post-step point polarization
                                  0,                     // Post-step point kinetic energy
-                                 0);                    // Post-step point charge
+                                 0,                     // Post-step point charge
+                                 0, -1);                // event and thread ID
       // The current track is killed by not enqueuing into the next activeQueue.
       break;
     }