diff --git a/examples/Example19/electrons.cu b/examples/Example19/electrons.cu index cfd0be44..3abba5fe 100644 --- a/examples/Example19/electrons.cu +++ b/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 ¤tTrack = 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 ¤tTrack = 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()}; diff --git a/examples/Example19/gammas.cu b/examples/Example19/gammas.cu index 2c9cb1f0..eb75061b 100644 --- a/examples/Example19/gammas.cu +++ b/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 ¤tTrack = 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 ¤tTrack = 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{¤tTrack.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);