Skip to content

Mark a track's survival more explicit #228

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Oct 11, 2022
Merged

Mark a track's survival more explicit #228

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 19 additions & 13 deletions examples/Example19/electrons.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}

Expand All @@ -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);
}
}

Expand All @@ -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;

Expand All @@ -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);
Expand All @@ -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()};
Expand Down
29 changes: 17 additions & 12 deletions examples/Example19/gammas.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}

Expand All @@ -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);
}
}

Expand All @@ -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;
}

Expand Down Expand Up @@ -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()};
Expand Down Expand Up @@ -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);
Expand Down