Provide secondary vector in G4Step

include/AdePT/core/AsyncAdePTTransport.icc

@@ -288,7 +288,12 @@ void AsyncAdePTTransport<IntegrationLayer>::ProcessGPUSteps(int threadId, int ev
 
   while ((range = async_adept_impl::GetGPUHitsFromBuffer(threadId, eventId, *fGPUstate, dataOnBuffer)).first !=
          nullptr) {
-    for (auto it = range.first; it != range.second; ++it) {
+
+    // Loop over returned GPU steps. The steps are ordered like this:
+    // [parent1][secondary1_1][secondary1_2][parent2][parent3][secondary3_1]
+    // In processing the steps of the parent tracks, the secondaries are also consumed and
+    // therefore the counter must advance by 1 + it->fNumSecondaries
+    for (auto it = range.first; it != range.second;) {
       // important sanity check: thread should only process its own hits and only from the current event
       if (it->threadId != threadId)
         std::cerr << "\033[1;31mError, threadId doesn't match it->threadId " << it->threadId << " threadId " << threadId
@@ -301,7 +306,10 @@ void AsyncAdePTTransport<IntegrationLayer>::ProcessGPUSteps(int threadId, int ev
                   << " ptype " << static_cast<short>(it->fParticleType) << " stepLimit / creator process "
                   << it->fStepLimProcessId << "\033[0m" << std::endl;
       }
-      integrationInstance.ProcessGPUStep(*it, fReturnAllSteps, fReturnFirstAndLastStep);
+      auto blockSize = 1 + it->fNumSecondaries;
+      std::span<const GPUHit> gpuStepWithSecondaries(it, blockSize);
+      integrationInstance.ProcessGPUStep(gpuStepWithSecondaries, fReturnAllSteps, fReturnFirstAndLastStep);
+      it += 1 + it->fNumSecondaries;
     }
     async_adept_impl::CloseGPUBuffer(threadId, *fGPUstate, range.first, dataOnBuffer);
   }

include/AdePT/core/PerEventScoringImpl.cuh

@@ -724,28 +724,25 @@ __device__ void RecordHit(AsyncAdePT::PerEventScoring * /*scoring*/, uint64_t aT
   // allocate hit slots: one for the parent and then one for each secondary
   auto slotStartIndex = AsyncAdePT::gHitScoringBuffer_dev.ReserveHitSlots(threadID, 1u + nSecondaries);
 
-  // NOTE: to be consistent with the previous implementation and to ensure that a new secondary arrives before the last
-  // step of the parent (otherwise, the hostTrackmapper might delete the parent, which is then not accessible in the
-  // processing of the secondary step anymore), the secondaries are processed before the parent. Next step will be to
-  // switch that order, such that the secondaries can be associated to the parent directly
+  // The ProcessGPUSteps on the Host expects the step of the parent track first, and then all secondaries
+  // that were generated in that step.
+  GPUHit &parentStep = AsyncAdePT::gHitScoringBuffer_dev.GetSlot(threadID, slotStartIndex);
+  // Fill the required data for the parent step
+  FillHit(parentStep, aTrackID, aParentID, stepLimProcessId, aParticleType, aStepLength, aTotalEnergyDeposit,
+          aTrackWeight, aPreState, aPrePosition, aPreMomentumDirection, aPreEKin, aPostState, aPostPosition,
+          aPostMomentumDirection, aPostEKin, aGlobalTime, aLocalTime, aPreGlobalTime, eventID, threadID, isLastStep,
+          stepCounter, nSecondaries);
 
   // Fill the steps for the secondaries
   for (unsigned int i = 0; i < nSecondaries; ++i) {
-    // The index should be slotStartInidex + 1u + i when the parent step is processed first
-    GPUHit &secondaryStep = AsyncAdePT::gHitScoringBuffer_dev.GetSlot(threadID, slotStartIndex + i);
+    // The index is the startIndex + 1 (for the parent) + i for the current secondary
+    GPUHit &secondaryStep = AsyncAdePT::gHitScoringBuffer_dev.GetSlot(threadID, slotStartIndex + 1u + i);
     FillHit(secondaryStep, secondaryData[i].trackId, aTrackID, stepLimProcessId, secondaryData[i].particleType,
             /*steplength*/ 0., /*energydeposit*/ 0., aTrackWeight, aPostState, aPostPosition, secondaryData[i].dir,
             secondaryData[i].eKin, aPostState, aPostPosition, secondaryData[i].dir, secondaryData[i].eKin, aGlobalTime,
-            /*localTime*/ 0., aGlobalTime, eventID, threadID, /*isLastStep*/ false, /*stepCounter*/ 0);
+            /*localTime*/ 0., aGlobalTime, eventID, threadID, /*isLastStep*/ false, /*stepCounter*/ 0,
+            /*nSecondaries*/ 0);
   }
-
-  // The index should simply be slotStartIndex when the parent is processed before the secondaries
-  GPUHit &parentStep = AsyncAdePT::gHitScoringBuffer_dev.GetSlot(threadID, slotStartIndex + nSecondaries);
-  // Fill the required data for the parent step
-  FillHit(parentStep, aTrackID, aParentID, stepLimProcessId, aParticleType, aStepLength, aTotalEnergyDeposit,
-          aTrackWeight, aPreState, aPrePosition, aPreMomentumDirection, aPreEKin, aPostState, aPostPosition,
-          aPostMomentumDirection, aPostEKin, aGlobalTime, aLocalTime, aPreGlobalTime, eventID, threadID, isLastStep,
-          stepCounter);
 }
 
 /// @brief Account for the number of produced secondaries

include/AdePT/core/ScoringCommons.hh

@@ -34,10 +34,10 @@ struct GPUHit {
   short fStepLimProcessId{-1};
   int fEventId{0};
   short threadId{-1};
-  // bool fFirstStepInVolume{false};
-  bool fLastStepOfTrack{false};
   unsigned short fStepCounter{0};
+  bool fLastStepOfTrack{false};
   char fParticleType{0}; // Particle type ID
+  unsigned char fNumSecondaries{0};
 };
 
 /// @brief Minimal data struct that is needed along with the parent track to provide the initial track information that
@@ -86,7 +86,8 @@ __device__ __forceinline__ void FillHit(
     vecgeom::Vector3D<double> const &aPrePosition, vecgeom::Vector3D<double> const &aPreMomentumDirection,
     double aPreEKin, vecgeom::NavigationState const &aPostState, vecgeom::Vector3D<double> const &aPostPosition,
     vecgeom::Vector3D<double> const &aPostMomentumDirection, double aPostEKin, double aGlobalTime, double aLocalTime,
-    double aPreGlobalTime, unsigned int eventID, short threadID, bool isLastStep, unsigned short stepCounter)
+    double aPreGlobalTime, unsigned int eventID, short threadID, bool isLastStep, unsigned short stepCounter,
+    unsigned char aNumSecondaries)
 {
   aGPUHit.fEventId = eventID;
   aGPUHit.threadId = threadID;
@@ -104,6 +105,7 @@ __device__ __forceinline__ void FillHit(
   aGPUHit.fGlobalTime         = aGlobalTime;
   aGPUHit.fLocalTime          = aLocalTime;
   aGPUHit.fPreGlobalTime      = aPreGlobalTime;
+  aGPUHit.fNumSecondaries     = aNumSecondaries;
   // Pre step point
   aGPUHit.fPreStepPoint.fNavigationState = aPreState;
   Copy3DVector(aPrePosition, aGPUHit.fPreStepPoint.fPosition);

include/AdePT/integration/AdePTGeant4Integration.hh

@@ -18,6 +18,7 @@
 #include <G4Event.hh>
 
 #include <unordered_map>
+#include <span>
 
 struct G4HepEmState;
 
@@ -63,7 +64,7 @@ public:
                              std::vector<G4LogicalVolume const *> &vecgeomLvToG4Map);
 
   /// @brief Reconstructs GPU hits on host and calls the user-defined sensitive detector code
-  void ProcessGPUStep(GPUHit const &hit, bool const callUserSteppingAction = false,
+  void ProcessGPUStep(std::span<const GPUHit> gpuSteps, bool const callUserSteppingAction = false,
                       bool const callUserTrackingaction = false);
 
   /// @brief Takes a range of tracks coming from the device and gives them back to Geant4
@@ -100,9 +101,19 @@ private:
 
   G4TouchableHandle MakeTouchableFromNavState(vecgeom::NavigationState const &navState) const;
 
-  void FillG4Step(GPUHit const *aGPUHit, G4Step *aG4Step, G4TouchableHandle &aPreG4TouchableHandle,
-                  G4TouchableHandle &aPostG4TouchableHandle, G4StepStatus aPreStepStatus, G4StepStatus aPostStepStatus,
-                  bool callUserTrackingAction, bool callUserSteppingAction) const;
+  /// @brief Construct the temporary secondary track that is attached to the secondary vector of the parent step
+  G4Track *ConstructSecondaryTrackInPlace(GPUHit const *secHit) const;
+
+  void InitSecondaryHostTrackDataFromParent(GPUHit const *secHit, HostTrackData &secTData, int g4ParentID,
+                                            G4TouchableHandle &preTouchable) const;
+
+  void FillG4Track(GPUHit const *aGPUHit, G4Track *aG4Track, const HostTrackData &hostTData,
+                   G4TouchableHandle &aPreG4TouchableHandle, G4TouchableHandle &aPostG4TouchableHandle) const;
+
+  void FillG4Step(GPUHit const *aGPUHit, G4Step *aG4Step, const HostTrackData &hostTData,
+                  G4TouchableHandle &aPreG4TouchableHandle, G4TouchableHandle &aPostG4TouchableHandle,
+                  G4StepStatus aPreStepStatus, G4StepStatus aPostStepStatus, bool callUserTrackingAction,
+                  bool callUserSteppingAction) const;
 
   void ReturnTrack(adeptint::TrackData const &track, unsigned int trackIndex, int debugLevel,
                    bool callUserActions = false) const;

include/AdePT/kernels/electrons.cuh

@@ -769,7 +769,7 @@ static __device__ __forceinline__ void TransportElectrons(ParticleManager &parti
       slotManager.MarkSlotForFreeing(slot);
     }
 
-    assert(nSecondaries <= 2);
+    assert(nSecondaries <= 3);
 
     // Record the step. Edep includes the continuous energy loss and edep from secondaries which were cut
     if ((energyDeposit > 0 && auxData.fSensIndex >= 0) || returnAllSteps ||