example19: Use interaction queues instead of SOAData

This revives an old concept from the toy examples 1, 2, and 4 (deleted
in commit 9230c4a) to have one queue per discrete interaction. This
appears to be equally fast as the quite complex scheduling kernel based
on SOAData.

Author: hahnjo

base/inc/AdePT/MParray.h

@@ -101,16 +101,14 @@ class MParray : protected copcore::VariableSizeObjectInterface<MParray, int> {
   const_reference operator[](size_t index) const { return fData[index]; }
 
   /** @brief Dispatch next free element, nullptr if none left */
-  __host__ __device__
-  __forceinline__
-  bool push_back(value_type val)
+  __host__ __device__ __forceinline__ int push_back(value_type val)
   {
     // Operation may fail if the max size is exceeded. Has to be checked by the user.
     int index = fNbooked.fetch_add(1);
-    if (index >= fCapacity) return false;
+    if (index >= fCapacity) return -1;
     fData[index] = val;
     fNused++;
-    return true;
+    return index;
   }
 
   /** @brief Check if container is fully distributed */

examples/Example19/electrons.cu

@@ -31,8 +31,8 @@
 template <bool IsElectron>
 static __device__ __forceinline__ void TransportElectrons(Track *electrons, const adept::MParray *active,
                                                           Secondaries &secondaries, adept::MParray *activeQueue,
-                                                          GlobalScoring *globalScoring,
-                                                          ScoringPerVolume *scoringPerVolume, SOAData soaData)
+                                                          Interactions interactions, GlobalScoring *globalScoring,
+                                                          ScoringPerVolume *scoringPerVolume)
 {
 #ifdef VECGEOM_FLOAT_PRECISION
   const Precision kPush = 10 * vecgeom::kTolerance;
@@ -57,9 +57,6 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
       if (push) activeQueue->push_back(globalSlot);
     };
 
-    // Signal that this globalSlot doesn't undergo an interaction (yet)
-    soaData.nextInteraction[i] = -1;
-
     // Init a track with the needed data to call into G4HepEm.
     G4HepEmElectronTrack elTrack;
     G4HepEmTrack *theTrack = elTrack.GetTrack();
@@ -285,151 +282,150 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
       continue;
     }
 
-    soaData.nextInteraction[i] = winnerProcessIndex;
+    interactions.queues[winnerProcessIndex]->push_back(globalSlot);
 
     survive(false);
   }
 }
 
 // Instantiate kernels for electrons and positrons.
 __global__ void TransportElectrons(Track *electrons, const adept::MParray *active, Secondaries secondaries,
-                                   adept::MParray *activeQueue, GlobalScoring *globalScoring,
-                                   ScoringPerVolume *scoringPerVolume, SOAData soaData)
+                                   adept::MParray *activeQueue, Interactions interactions, GlobalScoring *globalScoring,
+                                   ScoringPerVolume *scoringPerVolume)
 {
-  TransportElectrons</*IsElectron*/ true>(electrons, active, secondaries, activeQueue, globalScoring, scoringPerVolume,
-                                          soaData);
+  TransportElectrons</*IsElectron*/ true>(electrons, active, secondaries, activeQueue, interactions, globalScoring,
+                                          scoringPerVolume);
 }
 __global__ void TransportPositrons(Track *positrons, const adept::MParray *active, Secondaries secondaries,
-                                   adept::MParray *activeQueue, GlobalScoring *globalScoring,
-                                   ScoringPerVolume *scoringPerVolume, SOAData soaData)
+                                   adept::MParray *activeQueue, Interactions interactions, GlobalScoring *globalScoring,
+                                   ScoringPerVolume *scoringPerVolume)
 {
-  TransportElectrons</*IsElectron*/ false>(positrons, active, secondaries, activeQueue, globalScoring, scoringPerVolume,
-                                           soaData);
+  TransportElectrons</*IsElectron*/ false>(positrons, active, secondaries, activeQueue, interactions, globalScoring,
+                                           scoringPerVolume);
 }
 
 template <bool IsElectron, int ProcessIndex>
-__device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaData*/, int const /*soaSlot*/,
-                                    Track *particles, Secondaries secondaries, adept::MParray *activeQueue,
-                                    GlobalScoring *globalScoring, ScoringPerVolume *scoringPerVolume)
+__device__ void ElectronInteraction(Track *particles, const adept::MParray *active, Secondaries secondaries,
+                                    adept::MParray *activeQueue, GlobalScoring *globalScoring,
+                                    ScoringPerVolume *scoringPerVolume)
 {
-  Track &currentTrack = particles[globalSlot];
-  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;
-
-  RanluxppDouble newRNG{currentTrack.rngState.Branch()};
-  G4HepEmRandomEngine rnge{&currentTrack.rngState};
-
-  if constexpr (ProcessIndex == 0) {
-    // Invoke ionization (for e-/e+):
-    double deltaEkin = (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, energy, &rnge)
-                                    : G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, energy, &rnge);
-
-    double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
-    double dirSecondary[3];
-    G4HepEmElectronInteractionIoni::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
-
-    Track &secondary = secondaries.electrons.NextTrack();
-    atomicAdd(&globalScoring->numElectrons, 1);
-
-    secondary.InitAsSecondary(/*parent=*/currentTrack);
-    secondary.rngState = newRNG;
-    secondary.energy   = deltaEkin;
-    secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
-
-    currentTrack.energy = energy - deltaEkin;
-    currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
-    survive();
-  } else if constexpr (ProcessIndex == 1) {
-    // Invoke model for Bremsstrahlung: either SB- or Rel-Brem.
-    double logEnergy = std::log(energy);
-    double deltaEkin = energy < g4HepEmPars.fElectronBremModelLim
-                           ? G4HepEmElectronInteractionBrem::SampleETransferSB(&g4HepEmData, energy, logEnergy,
-                                                                               theMCIndex, &rnge, IsElectron)
-                           : G4HepEmElectronInteractionBrem::SampleETransferRB(&g4HepEmData, energy, logEnergy,
-                                                                               theMCIndex, &rnge, IsElectron);
-
-    double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
-    double dirSecondary[3];
-    G4HepEmElectronInteractionBrem::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
-
-    Track &gamma = secondaries.gammas.NextTrack();
-    atomicAdd(&globalScoring->numGammas, 1);
-
-    gamma.InitAsSecondary(/*parent=*/currentTrack);
-    gamma.rngState = newRNG;
-    gamma.energy   = deltaEkin;
-    gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
-
-    currentTrack.energy = energy - deltaEkin;
-    currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
-    survive();
-  } else if constexpr (ProcessIndex == 2) {
-    // Invoke annihilation (in-flight) for e+
-    double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
-    double theGamma1Ekin, theGamma2Ekin;
-    double theGamma1Dir[3], theGamma2Dir[3];
-    G4HepEmPositronInteractionAnnihilation::SampleEnergyAndDirectionsInFlight(
-        energy, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge);
-
-    Track &gamma1 = secondaries.gammas.NextTrack();
-    Track &gamma2 = secondaries.gammas.NextTrack();
-    atomicAdd(&globalScoring->numGammas, 2);
-
-    gamma1.InitAsSecondary(/*parent=*/currentTrack);
-    gamma1.rngState = newRNG;
-    gamma1.energy   = theGamma1Ekin;
-    gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]);
-
-    gamma2.InitAsSecondary(/*parent=*/currentTrack);
-    // Reuse the RNG state of the dying track.
-    gamma2.rngState = currentTrack.rngState;
-    gamma2.energy   = theGamma2Ekin;
-    gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]);
-
-    // The current track is killed by not enqueuing into the next activeQueue.
+  int activeSize = active->size();
+  for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
+    const int globalSlot = (*active)[i];
+    Track &currentTrack  = particles[globalSlot];
+    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;
+
+    RanluxppDouble newRNG{currentTrack.rngState.Branch()};
+    G4HepEmRandomEngine rnge{&currentTrack.rngState};
+
+    if constexpr (ProcessIndex == 0) {
+      // Invoke ionization (for e-/e+):
+      double deltaEkin = (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, energy, &rnge)
+                                      : G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, energy, &rnge);
+
+      double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
+      double dirSecondary[3];
+      G4HepEmElectronInteractionIoni::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
+
+      Track &secondary = secondaries.electrons.NextTrack();
+      atomicAdd(&globalScoring->numElectrons, 1);
+
+      secondary.InitAsSecondary(/*parent=*/currentTrack);
+      secondary.rngState = newRNG;
+      secondary.energy   = deltaEkin;
+      secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
+
+      currentTrack.energy = energy - deltaEkin;
+      currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
+      survive();
+    } else if constexpr (ProcessIndex == 1) {
+      // Invoke model for Bremsstrahlung: either SB- or Rel-Brem.
+      double logEnergy = std::log(energy);
+      double deltaEkin = energy < g4HepEmPars.fElectronBremModelLim
+                             ? G4HepEmElectronInteractionBrem::SampleETransferSB(&g4HepEmData, energy, logEnergy,
+                                                                                 theMCIndex, &rnge, IsElectron)
+                             : G4HepEmElectronInteractionBrem::SampleETransferRB(&g4HepEmData, energy, logEnergy,
+                                                                                 theMCIndex, &rnge, IsElectron);
+
+      double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
+      double dirSecondary[3];
+      G4HepEmElectronInteractionBrem::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
+
+      Track &gamma = secondaries.gammas.NextTrack();
+      atomicAdd(&globalScoring->numGammas, 1);
+
+      gamma.InitAsSecondary(/*parent=*/currentTrack);
+      gamma.rngState = newRNG;
+      gamma.energy   = deltaEkin;
+      gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
+
+      currentTrack.energy = energy - deltaEkin;
+      currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
+      survive();
+    } else if constexpr (ProcessIndex == 2) {
+      // Invoke annihilation (in-flight) for e+
+      double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
+      double theGamma1Ekin, theGamma2Ekin;
+      double theGamma1Dir[3], theGamma2Dir[3];
+      G4HepEmPositronInteractionAnnihilation::SampleEnergyAndDirectionsInFlight(
+          energy, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge);
+
+      Track &gamma1 = secondaries.gammas.NextTrack();
+      Track &gamma2 = secondaries.gammas.NextTrack();
+      atomicAdd(&globalScoring->numGammas, 2);
+
+      gamma1.InitAsSecondary(/*parent=*/currentTrack);
+      gamma1.rngState = newRNG;
+      gamma1.energy   = theGamma1Ekin;
+      gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]);
+
+      gamma2.InitAsSecondary(/*parent=*/currentTrack);
+      // Reuse the RNG state of the dying track.
+      gamma2.rngState = currentTrack.rngState;
+      gamma2.energy   = theGamma2Ekin;
+      gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]);
+
+      // The current track is killed by not enqueuing into the next activeQueue.
+    }
   }
 }
 
 __global__ void IonizationEl(Track *particles, const adept::MParray *active, Secondaries secondaries,
                              adept::MParray *activeQueue, GlobalScoring *globalScoring,
-                             ScoringPerVolume *scoringPerVolume, SOAData const soaData)
+                             ScoringPerVolume *scoringPerVolume)
 {
-  InteractionLoop<0>(&ElectronInteraction<true, 0>, active, soaData, particles, secondaries, activeQueue, globalScoring,
-                     scoringPerVolume);
+  ElectronInteraction<true, 0>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
 }
 __global__ void BremsstrahlungEl(Track *particles, const adept::MParray *active, Secondaries secondaries,
                                  adept::MParray *activeQueue, GlobalScoring *globalScoring,
-                                 ScoringPerVolume *scoringPerVolume, SOAData const soaData)
+                                 ScoringPerVolume *scoringPerVolume)
 {
-  InteractionLoop<1>(&ElectronInteraction<true, 1>, active, soaData, particles, secondaries, activeQueue, globalScoring,
-                     scoringPerVolume);
+  ElectronInteraction<true, 1>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
 }
 
 __global__ void IonizationPos(Track *particles, const adept::MParray *active, Secondaries secondaries,
                               adept::MParray *activeQueue, GlobalScoring *globalScoring,
-                              ScoringPerVolume *scoringPerVolume, SOAData const soaData)
+                              ScoringPerVolume *scoringPerVolume)
 {
-  InteractionLoop<0>(&ElectronInteraction<false, 0>, active, soaData, particles, secondaries, activeQueue,
-                     globalScoring, scoringPerVolume);
+  ElectronInteraction<false, 0>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
 }
 __global__ void BremsstrahlungPos(Track *particles, const adept::MParray *active, Secondaries secondaries,
                                   adept::MParray *activeQueue, GlobalScoring *globalScoring,
-                                  ScoringPerVolume *scoringPerVolume, SOAData const soaData)
+                                  ScoringPerVolume *scoringPerVolume)
 {
-  InteractionLoop<1>(&ElectronInteraction<false, 1>, active, soaData, particles, secondaries, activeQueue,
-                     globalScoring, scoringPerVolume);
+  ElectronInteraction<false, 1>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
 }
 __global__ void AnnihilationPos(Track *particles, const adept::MParray *active, Secondaries secondaries,
                                 adept::MParray *activeQueue, GlobalScoring *globalScoring,
-                                ScoringPerVolume *scoringPerVolume, SOAData const soaData)
+                                ScoringPerVolume *scoringPerVolume)
 {
-  InteractionLoop<2>(&ElectronInteraction<false, 2>, active, soaData, particles, secondaries, activeQueue,
-                     globalScoring, scoringPerVolume);
+  ElectronInteraction<false, 2>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
 }