Mark a track's survival more explicit

Adds a lambda function survive() which pushes the track back into the active queue.

Author: bernhardmgruber

examples/Example19/electrons.cu

@@ -47,11 +47,15 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
     const int globalSlot = (*active)[i];
     Track &currentTrack  = electrons[globalSlot];
-    auto volume          = currentTrack.navState.Top();
-    int volumeID         = volume->id();
+    const auto volume    = currentTrack.navState.Top();
+    const int volumeID   = volume->id();
     // the MCC vector is indexed by the logical volume id
-    int lvolID     = volume->GetLogicalVolume()->id();
-    int theMCIndex = MCIndex[lvolID];
+    const int lvolID     = volume->GetLogicalVolume()->id();
+    const int theMCIndex = MCIndex[lvolID];
+
+    auto survive = [&](bool push = true) {
+      if (push) activeQueue->push_back(globalSlot);
+    };
 
     // Signal that this globalSlot doesn't undergo an interaction (yet)
     soaData.nextInteraction[i] = -1;
@@ -252,21 +256,21 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
 
       // Kill the particle if it left the world.
       if (nextState.Top() != nullptr) {
-        activeQueue->push_back(globalSlot);
         BVHNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, nextState);
 
         // Move to the next boundary.
         currentTrack.navState = nextState;
+        survive();
       }
       continue;
     } else if (!propagated || restrictedPhysicalStepLength) {
       // Did not yet reach the interaction point due to error in the magnetic
       // field propagation. Try again next time.
-      activeQueue->push_back(globalSlot);
+      survive();
       continue;
     } else if (winnerProcessIndex < 0) {
       // No discrete process, move on.
-      activeQueue->push_back(globalSlot);
+      survive();
       continue;
     }
 
@@ -277,11 +281,13 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
     // Check if a delta interaction happens instead of the real discrete process.
     if (G4HepEmElectronManager::CheckDelta(&g4HepEmData, theTrack, currentTrack.Uniform())) {
       // A delta interaction happened, move on.
-      activeQueue->push_back(globalSlot);
+      survive();
       continue;
     }
 
     soaData.nextInteraction[i] = winnerProcessIndex;
+
+    survive(false);
   }
 }
 
@@ -307,11 +313,13 @@ __device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaD
                                     GlobalScoring *globalScoring, ScoringPerVolume *scoringPerVolume)
 {
   Track &currentTrack = particles[globalSlot];
-  auto volume         = currentTrack.navState.Top();
+  const auto volume   = currentTrack.navState.Top();
   // the MCC vector is indexed by the logical volume id
   const int lvolID     = volume->GetLogicalVolume()->id();
   const int theMCIndex = MCIndex[lvolID];
 
+  auto survive = [&] { activeQueue->push_back(globalSlot); };
+
   const double energy   = currentTrack.energy;
   const double theElCut = g4HepEmData.fTheMatCutData->fMatCutData[theMCIndex].fSecElProdCutE;
 
@@ -337,8 +345,7 @@ __device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaD
 
     currentTrack.energy = energy - deltaEkin;
     currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
-    // The current track continues to live.
-    activeQueue->push_back(globalSlot);
+    survive();
   } else if constexpr (ProcessIndex == 1) {
     // Invoke model for Bremsstrahlung: either SB- or Rel-Brem.
     double logEnergy = std::log(energy);
@@ -362,8 +369,7 @@ __device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaD
 
     currentTrack.energy = energy - deltaEkin;
     currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
-    // The current track continues to live.
-    activeQueue->push_back(globalSlot);
+    survive();
   } else if constexpr (ProcessIndex == 2) {
     // Invoke annihilation (in-flight) for e+
     double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};

examples/Example19/gammas.cu

@@ -32,11 +32,15 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
   for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
     const int slot      = (*active)[i];
     Track &currentTrack = gammas[slot];
-    auto volume         = currentTrack.navState.Top();
-    int volumeID        = volume->id();
+    const auto volume   = currentTrack.navState.Top();
+    const int volumeID  = volume->id();
     // the MCC vector is indexed by the logical volume id
-    int lvolID     = volume->GetLogicalVolume()->id();
-    int theMCIndex = MCIndex[lvolID];
+    const int lvolID     = volume->GetLogicalVolume()->id();
+    const int theMCIndex = MCIndex[lvolID];
+
+    auto survive = [&](bool push = true) {
+      if (push) activeQueue->push_back(slot);
+    };
 
     // Signal that this slot doesn't undergo an interaction (yet)
     soaData.nextInteraction[i] = -1;
@@ -95,16 +99,16 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
 
       // Kill the particle if it left the world.
       if (nextState.Top() != nullptr) {
-        activeQueue->push_back(slot);
         BVHNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, nextState);
 
         // Move to the next boundary.
         currentTrack.navState = nextState;
+        survive();
       }
       continue;
     } else if (winnerProcessIndex < 0) {
       // No discrete process, move on.
-      activeQueue->push_back(slot);
+      survive();
       continue;
     }
 
@@ -114,6 +118,7 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
 
     soaData.nextInteraction[i] = winnerProcessIndex;
     soaData.gamma_PEmxSec[i] = gammaTrack.GetPEmxSec();
+    survive(false);
   }
 }
 
@@ -123,20 +128,22 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
                                  ScoringPerVolume *scoringPerVolume)
 {
   Track &currentTrack = particles[globalSlot];
-  auto volume         = currentTrack.navState.Top();
+  const auto volume   = currentTrack.navState.Top();
   const int volumeID  = volume->id();
   // the MCC vector is indexed by the logical volume id
   const int lvolID     = volume->GetLogicalVolume()->id();
   const int theMCIndex = MCIndex[lvolID];
   const auto energy    = currentTrack.energy;
 
+  auto survive = [&] { activeQueue->push_back(globalSlot); };
+
   RanluxppDouble newRNG{currentTrack.rngState.Branch()};
   G4HepEmRandomEngine rnge{&currentTrack.rngState};
 
   if constexpr (ProcessIndex == 0) {
     // Invoke gamma conversion to e-/e+ pairs, if the energy is above the threshold.
     if (energy < 2 * copcore::units::kElectronMassC2) {
-      activeQueue->push_back(globalSlot);
+      survive();
       return;
     }
 
@@ -171,7 +178,7 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
     // Invoke Compton scattering of gamma.
     constexpr double LowEnergyThreshold = 100 * copcore::units::eV;
     if (energy < LowEnergyThreshold) {
-      activeQueue->push_back(globalSlot);
+      survive();
       return;
     }
     const double origDirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
@@ -200,9 +207,7 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
     if (newEnergyGamma > LowEnergyThreshold) {
       currentTrack.energy = newEnergyGamma;
       currentTrack.dir    = newDirGamma;
-
-      // The current track continues to live.
-      activeQueue->push_back(globalSlot);
+      survive();
     } else {
       atomicAdd(&globalScoring->energyDeposit, newEnergyGamma);
       atomicAdd(&scoringPerVolume->energyDeposit[volumeID], newEnergyGamma);