From 1331728bbd64eb85107181bc6813edcfbb22a5ba Mon Sep 17 00:00:00 2001 From: JuanGonzalezCaminero Date: Tue, 24 Sep 2024 16:33:22 +0200 Subject: [PATCH] Some cleanup --- CMakeLists.txt | 8 +- include/AdePT/core/AdePTTransport.cuh | 64 +- include/AdePT/core/AdePTTransportStruct.cuh | 8 +- include/AdePT/core/Track.cuh | 25 +- .../AdePT/kernels/electrons_experimental.cuh | 933 ------------------ include/AdePT/kernels/gammas_experimental.cuh | 411 -------- 6 files changed, 53 insertions(+), 1396 deletions(-) delete mode 100644 include/AdePT/kernels/electrons_experimental.cuh delete mode 100644 include/AdePT/kernels/gammas_experimental.cuh diff --git a/CMakeLists.txt b/CMakeLists.txt index 18418eed..a0dc8cea 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,6 +48,7 @@ set(CMAKE_INCLUDE_DIRECTORIES_PROJECT_BEFORE ON) # Options option(ADEPT_USE_SURF "Enable surface model navigation on GPU" OFF) option(ADEPT_USE_SURF_SINGLE "Use surface model in single precision" OFF) +option(USE_SPLIT_KERNELS "Run split version of the transport kernels" OFF) option(DEBUG_SINGLE_THREAD "Run transport kernels in single thread mode" OFF) option(WITH_FLUCT "Switch on the energy loss fluctuations" OFF) @@ -83,6 +84,11 @@ endif() if(NOT TARGET VecGeom::vgdml) message(FATAL_ERROR "AdePT requires VecGeom compiled with GDML support") endif() +# Run split kernels +if (USE_SPLIT_KERNELS) + add_compile_definitions(USE_SPLIT_KERNELS) + message(STATUS "${Green}AdePT will run with split kernels${ColorReset}") +endif() # Debugging in single-thread mode if (DEBUG_SINGLE_THREAD) add_compile_definitions("$<$:DEBUG_SINGLE_THREAD>") @@ -153,8 +159,6 @@ if(NOT WITH_FLUCT) add_compile_definitions(NOFLUCTUATION) endif() -# string(APPEND CMAKE_CUDA_FLAGS " -Xptxas=-v") - #----------------------------------------------------------------------------# # Build Targets #----------------------------------------------------------------------------# diff --git a/include/AdePT/core/AdePTTransport.cuh b/include/AdePT/core/AdePTTransport.cuh index 05672bb4..4b0a4c99 100644 --- a/include/AdePT/core/AdePTTransport.cuh +++ b/include/AdePT/core/AdePTTransport.cuh @@ -9,18 +9,14 @@ #include #include #include -// #include -// #include - - -#include -// #include -// #include - -#include -// #include -// #include +#ifndef USE_SPLIT_KERNELS +#include +#include +#else +#include +#include +#endif #include #ifdef VECGEOM_ENABLE_CUDA @@ -272,12 +268,14 @@ GPUstate *InitializeGPU(adeptint::TrackBuffer &buffer, int capacity, int maxbatc COPCORE_CUDA_CHECK(cudaEventCreate(&gpuState.particles[i].event)); } +#ifdef USE_SPLIT_KERNELS // Init HepEM tracks // e+ / e- COPCORE_CUDA_CHECK(cudaMalloc(&gpuState.hepEMBuffers_d.electronsHepEm, capacity * sizeof(G4HepEmElectronTrack))); COPCORE_CUDA_CHECK(cudaMalloc(&gpuState.hepEMBuffers_d.positronsHepEm, capacity * sizeof(G4HepEmElectronTrack))); // Gammas COPCORE_CUDA_CHECK(cudaMalloc(&gpuState.hepEMBuffers_d.gammasHepEm, capacity * sizeof(G4HepEmGammaTrack))); +#endif InitLeakedQueues<<<1, 1, 0, gpuState.stream>>>(gpuState.allmgr_d, kQueueSize); COPCORE_CUDA_CHECK(cudaDeviceSynchronize()); @@ -305,9 +303,11 @@ void FreeGPU(GPUstate &gpuState, G4HepEmState *g4hepem_state) COPCORE_CUDA_CHECK(cudaFreeHost(gpuState.stats)); COPCORE_CUDA_CHECK(cudaFree(gpuState.toDevice_dev)); - COPCORE_CUDA_CHECK(cudaFree(gpuState.hepEMBuffers_d.electronsHepEm)); +#ifdef USE_SPLIT_KERNELS + COPCORE_CUDA_CHECK(cudaFree(gpuState.hepEMBuffers_d.electronsHepEm)); COPCORE_CUDA_CHECK(cudaFree(gpuState.hepEMBuffers_d.positronsHepEm)); COPCORE_CUDA_CHECK(cudaFree(gpuState.hepEMBuffers_d.gammasHepEm)); +#endif COPCORE_CUDA_CHECK(cudaStreamDestroy(gpuState.stream)); @@ -402,13 +402,11 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & transportBlocks = (numElectrons + TransportThreads - 1) / TransportThreads; transportBlocks = std::min(transportBlocks, MaxBlocks); #endif - // TransportElectrons<<>>( - // electrons.trackmgr, secondaries, electrons.leakedTracks, scoring_dev, - // VolAuxArray::GetInstance().fAuxData_dev); - // TransportElectrons<<>>( - // electrons.trackmgr, gpuState.hepEMBuffers_d.electronsHepEm, secondaries, electrons.leakedTracks, scoring_dev, - // VolAuxArray::GetInstance().fAuxData_dev); - +#ifndef USE_SPLIT_KERNELS + TransportElectrons<<>>( + electrons.trackmgr, secondaries, electrons.leakedTracks, scoring_dev, + VolAuxArray::GetInstance().fAuxData_dev); +#else ElectronPhysics1<<>>( electrons.trackmgr, gpuState.hepEMBuffers_d.electronsHepEm, VolAuxArray::GetInstance().fAuxData_dev ); @@ -421,7 +419,7 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & ElectronInteractions<<>>( electrons.trackmgr, gpuState.hepEMBuffers_d.electronsHepEm, secondaries, electrons.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); - +#endif COPCORE_CUDA_CHECK(cudaEventRecord(electrons.event, electrons.stream)); COPCORE_CUDA_CHECK(cudaStreamWaitEvent(gpuState.stream, electrons.event, 0)); } @@ -433,13 +431,11 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & transportBlocks = (numPositrons + TransportThreads - 1) / TransportThreads; transportBlocks = std::min(transportBlocks, MaxBlocks); #endif - // TransportPositrons<<>>( - // positrons.trackmgr, secondaries, positrons.leakedTracks, scoring_dev, - // VolAuxArray::GetInstance().fAuxData_dev); - // TransportElectrons<<>>( - // positrons.trackmgr, gpuState.hepEMBuffers_d.positronsHepEm, secondaries, positrons.leakedTracks, scoring_dev, - // VolAuxArray::GetInstance().fAuxData_dev); - +#ifndef USE_SPLIT_KERNELS + TransportPositrons<<>>( + positrons.trackmgr, secondaries, positrons.leakedTracks, scoring_dev, + VolAuxArray::GetInstance().fAuxData_dev); +#else ElectronPhysics1<<>>( positrons.trackmgr, gpuState.hepEMBuffers_d.positronsHepEm, VolAuxArray::GetInstance().fAuxData_dev ); @@ -452,7 +448,7 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & ElectronInteractions<<>>( positrons.trackmgr, gpuState.hepEMBuffers_d.positronsHepEm, secondaries, positrons.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); - +#endif COPCORE_CUDA_CHECK(cudaEventRecord(positrons.event, positrons.stream)); COPCORE_CUDA_CHECK(cudaStreamWaitEvent(gpuState.stream, positrons.event, 0)); } @@ -464,6 +460,10 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & transportBlocks = (numGammas + TransportThreads - 1) / TransportThreads; transportBlocks = std::min(transportBlocks, MaxBlocks); #endif +#ifndef USE_SPLIT_KERNELS + TransportGammas<<>>( + gammas.trackmgr, secondaries, gammas.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); +#else GammaPhysics1<<>>( gammas.trackmgr, gpuState.hepEMBuffers_d.gammasHepEm, VolAuxArray::GetInstance().fAuxData_dev ); @@ -476,13 +476,7 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & GammaPhysics2<<>>( gammas.trackmgr, gpuState.hepEMBuffers_d.gammasHepEm, secondaries, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev ); - - // TransportGammas<<>>( - // gammas.trackmgr, secondaries, gammas.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); - - // TransportGammas<<>>( - // gammas.trackmgr, gpuState.hepEMBuffers_d.gammasHepEm, secondaries, gammas.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); - +#endif COPCORE_CUDA_CHECK(cudaEventRecord(gammas.event, gammas.stream)); COPCORE_CUDA_CHECK(cudaStreamWaitEvent(gpuState.stream, gammas.event, 0)); } diff --git a/include/AdePT/core/AdePTTransportStruct.cuh b/include/AdePT/core/AdePTTransportStruct.cuh index 52eb75e7..ff80e9f7 100644 --- a/include/AdePT/core/AdePTTransportStruct.cuh +++ b/include/AdePT/core/AdePTTransportStruct.cuh @@ -14,8 +14,10 @@ #include #include +#ifdef USE_SPLIT_KERNELS #include #include +#endif #ifdef __CUDA_ARCH__ // Define inline implementations of the RNG methods for the device. @@ -68,11 +70,13 @@ struct AllTrackManagers { MParrayTracks *leakedTracks[ParticleType::NumParticleTypes]; }; +#ifdef USE_SPLIT_KERNELS struct HepEmBuffers { G4HepEmElectronTrack *electronsHepEm; G4HepEmElectronTrack *positronsHepEm; G4HepEmGammaTrack *gammasHepEm; }; +#endif // A data structure to transfer statistics after each iteration. struct Stats { @@ -87,7 +91,9 @@ struct GPUstate { ParticleType particles[ParticleType::NumParticleTypes]; AllTrackManagers allmgr_h; ///< Host pointers for track managers AllTrackManagers allmgr_d; ///< Device pointers for track managers - HepEmBuffers hepEMBuffers_d; +#ifdef USE_SPLIT_KERNELS + HepEmBuffers hepEMBuffers_d; +#endif // Create a stream to synchronize kernels of all particle types. cudaStream_t stream; ///< all-particle sync stream TrackData *toDevice_dev{nullptr}; ///< toDevice buffer of tracks diff --git a/include/AdePT/core/Track.cuh b/include/AdePT/core/Track.cuh index 88734086..f41d6f0b 100644 --- a/include/AdePT/core/Track.cuh +++ b/include/AdePT/core/Track.cuh @@ -19,11 +19,9 @@ struct Track { using Precision = vecgeom::Precision; int parentID{0}; // Stores the track id of the initial particle given to AdePT - long hitsurfID{0}; RanluxppDouble rngState; double eKin; - double preStepEKin; double numIALeft[3]; double initialRange; double dynamicRangeFactor; @@ -34,30 +32,29 @@ struct Track { double properTime{0}; vecgeom::Vector3D pos; - vecgeom::Vector3D preStepPos; vecgeom::Vector3D dir; - vecgeom::Vector3D preStepDir; vecgeom::NavigationState navState; + +#ifdef USE_SPLIT_KERNELS + RanluxppDouble newRNG; + + // Variables used to store track info needed for scoring + double preStepEKin; + vecgeom::Vector3D preStepPos; + vecgeom::Vector3D preStepDir; vecgeom::NavigationState nextState; vecgeom::NavigationState preStepNavState; // Variables used to store navigation results + long hitsurfID{0}; bool propagated{false}; double geometryStepLength{0}; - // Todo: check whether it's needed to keep safety here or we can use the one stored in the HepEM track double safety{0}; - RanluxppDouble newRNG; - - // Variables used to store physics results from G4HepEM - // double geometricalStepLengthFromPhysics{0}; - // int winnerProcessIndex{0}; - // double physicalStepLength{0}; - // double preStepMFPs[3]; - // double PEmxSec{0}; - // G4HepEmMSCTrackData mscData; + // Variables used to store results from G4HepEM bool restrictedPhysicalStepLength{false}; bool stopped{false}; +#endif __host__ __device__ double Uniform() { return rngState.Rndm(); } diff --git a/include/AdePT/kernels/electrons_experimental.cuh b/include/AdePT/kernels/electrons_experimental.cuh deleted file mode 100644 index bb28b6c8..00000000 --- a/include/AdePT/kernels/electrons_experimental.cuh +++ /dev/null @@ -1,933 +0,0 @@ -// SPDX-FileCopyrightText: 2022 CERN -// SPDX-License-Identifier: Apache-2.0 - -#include -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -// Pull in implementation. -#include -#include -#include -#include -#include -#include -#include -#include - -using VolAuxData = adeptint::VolAuxData; - -// Compute velocity based on the kinetic energy of the particle -__device__ double GetVelocity(double eKin) -{ - // Taken from G4DynamicParticle::ComputeBeta - double T = eKin / copcore::units::kElectronMassC2; - double beta = sqrt(T * (T + 2.)) / (T + 1.0); - return copcore::units::kCLight * beta; -} - -template -__global__ void ElectronPhysics1(adept::TrackManager *electrons, G4HepEmElectronTrack *hepEMTracks, VolAuxData const *auxDataArray) -{ - constexpr int Charge = IsElectron ? -1 : 1; - constexpr double restMass = copcore::units::kElectronMassC2; - fieldPropagatorConstBz fieldPropagatorBz(BzFieldValue); - - int activeSize = electrons->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*electrons->fActiveTracks)[i]; - Track ¤tTrack = (*electrons)[slot]; - currentTrack.preStepEKin = currentTrack.eKin; - currentTrack.preStepPos = currentTrack.pos; - currentTrack.preStepDir = currentTrack.dir; - currentTrack.preStepNavState = currentTrack.navState; - // the MCC vector is indexed by the logical volume id -#ifndef ADEPT_USE_SURF - const int lvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -#else - const int lvolID = currentTrack.navState.GetLogicalId(); -#endif - VolAuxData const &auxData = auxDataArray[lvolID]; - - // Init a track with the needed data to call into G4HepEm. - G4HepEmElectronTrack &elTrack = hepEMTracks[slot]; - elTrack.ReSet(); - G4HepEmTrack *theTrack = elTrack.GetTrack(); - theTrack->SetEKin(currentTrack.eKin); - theTrack->SetMCIndex(auxData.fMCIndex); - theTrack->SetOnBoundary(currentTrack.navState.IsOnBoundary()); - theTrack->SetCharge(Charge); - G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData(); - mscData->fIsFirstStep = currentTrack.initialRange < 0; - mscData->fInitialRange = currentTrack.initialRange; - mscData->fDynamicRangeFactor = currentTrack.dynamicRangeFactor; - mscData->fTlimitMin = currentTrack.tlimitMin; - - // Prepare a branched RNG state while threads are synchronized. Even if not - // used, this provides a fresh round of random numbers and reduces thread - // divergence because the RNG state doesn't need to be advanced later. - currentTrack.newRNG = RanluxppDouble(currentTrack.rngState.BranchNoAdvance()); - - // Compute safety, needed for MSC step limit. - currentTrack.safety = 0; - if (!currentTrack.navState.IsOnBoundary()) { - currentTrack.safety = AdePTNavigator::ComputeSafety(currentTrack.pos, currentTrack.navState); - } - theTrack->SetSafety(currentTrack.safety); - - G4HepEmRandomEngine rnge(¤tTrack.rngState); - - // Sample the `number-of-interaction-left` and put it into the track. - for (int ip = 0; ip < 3; ++ip) { - double numIALeft = currentTrack.numIALeft[ip]; - if (numIALeft <= 0) { - numIALeft = -std::log(currentTrack.Uniform()); - // currentTrack.numIALeft[ip] = numIALeft; - } - theTrack->SetNumIALeft(numIALeft, ip); - } - - G4HepEmElectronManager::HowFarToDiscreteInteraction(&g4HepEmData, &g4HepEmPars, &elTrack); - - currentTrack.restrictedPhysicalStepLength = false; - if (BzFieldValue != 0) { - const double momentumMag = sqrt(currentTrack.eKin * (currentTrack.eKin + 2.0 * restMass)); - // Distance along the track direction to reach the maximum allowed error - const double safeLength = fieldPropagatorBz.ComputeSafeLength(momentumMag, Charge, currentTrack.dir); - - constexpr int MaxSafeLength = 10; - double limit = MaxSafeLength * safeLength; - limit = currentTrack.safety > limit ? currentTrack.safety : limit; - - double physicalStepLength = elTrack.GetPStepLength(); - if (physicalStepLength > limit) { - physicalStepLength = limit; - currentTrack.restrictedPhysicalStepLength = true; - elTrack.SetPStepLength(physicalStepLength); - // Note: We are limiting the true step length, which is converted to - // a shorter geometry step length in HowFarToMSC. In that sense, the - // limit is an over-approximation, but that is fine for our purpose. - } - } - - G4HepEmElectronManager::HowFarToMSC(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge); - - // Remember MSC values for the next step(s). - currentTrack.initialRange = mscData->fInitialRange; - currentTrack.dynamicRangeFactor = mscData->fDynamicRangeFactor; - currentTrack.tlimitMin = mscData->fTlimitMin; - } -} - -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// - -template -static __global__ void ElectronTransport1(adept::TrackManager *electrons, G4HepEmElectronTrack *hepEMTracks) -{ -#ifdef VECGEOM_FLOAT_PRECISION - const Precision kPush = 10 * vecgeom::kTolerance; -#else - const Precision kPush = 0.; -#endif - constexpr int Charge = IsElectron ? -1 : 1; - constexpr double restMass = copcore::units::kElectronMassC2; - fieldPropagatorConstBz fieldPropagatorBz(BzFieldValue); - - int activeSize = electrons->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*electrons->fActiveTracks)[i]; - Track ¤tTrack = (*electrons)[slot]; - - // Retrieve HepEM track - G4HepEmElectronTrack &elTrack = hepEMTracks[slot]; - G4HepEmTrack *theTrack = elTrack.GetTrack(); - - // Check if there's a volume boundary in between. - currentTrack.propagated = true; - currentTrack.hitsurfID = -1; - // double geometryStepLength; - // vecgeom::NavigationState nextState; - if (BzFieldValue != 0) { - currentTrack.geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume( - currentTrack.eKin, restMass, Charge, theTrack->GetGStepLength(), currentTrack.pos, currentTrack.dir, - currentTrack.navState, currentTrack.nextState, currentTrack.hitsurfID, currentTrack.propagated, currentTrack.safety); - } else { -#ifdef ADEPT_USE_SURF - currentTrack.geometryStepLength = - AdePTNavigator::ComputeStepAndNextVolume(currentTrack.pos, currentTrack.dir, theTrack->GetGStepLength(), - currentTrack.navState, currentTrack.nextState, currentTrack.hitsurfID, kPush); -#else - currentTrack.geometryStepLength = - AdePTNavigator::ComputeStepAndNextVolume(currentTrack.pos, currentTrack.dir, theTrack->GetGStepLength(), - currentTrack.navState, currentTrack.nextState, kPush); -#endif - currentTrack.pos += currentTrack.geometryStepLength * currentTrack.dir; - } - - // Set boundary state in navState so the next step and secondaries get the - // correct information (navState = nextState only if relocated - // in case of a boundary; see below) - currentTrack.navState.SetBoundaryState(currentTrack.nextState.IsOnBoundary()); - - // Propagate information from geometrical step to MSC. - theTrack->SetDirection(currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()); - theTrack->SetGStepLength(currentTrack.geometryStepLength); - theTrack->SetOnBoundary(currentTrack.nextState.IsOnBoundary()); - } -} - -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// - -// May work well to separate MSC1 (Continuous effects) and MSC2 (Checks + safety) -template -static __global__ void MSC1(adept::TrackManager *electrons, G4HepEmElectronTrack *hepEMTracks) -{ - constexpr double restMass = copcore::units::kElectronMassC2; - int activeSize = electrons->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*electrons->fActiveTracks)[i]; - Track ¤tTrack = (*electrons)[slot]; - - // Init a track with the needed data to call into G4HepEm. - G4HepEmElectronTrack &elTrack = hepEMTracks[slot]; - G4HepEmTrack *theTrack = elTrack.GetTrack(); - G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData(); - - - G4HepEmRandomEngine rnge(¤tTrack.rngState); - - // Apply continuous effects. - currentTrack.stopped = G4HepEmElectronManager::PerformContinuous(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge); - - // Collect the direction change and displacement by MSC. - const double *direction = theTrack->GetDirection(); - currentTrack.dir.Set(direction[0], direction[1], direction[2]); - if (!currentTrack.nextState.IsOnBoundary()) { - const double *mscDisplacement = mscData->GetDisplacement(); - vecgeom::Vector3D displacement(mscDisplacement[0], mscDisplacement[1], mscDisplacement[2]); - const double dLength2 = displacement.Length2(); - constexpr double kGeomMinLength = 5 * copcore::units::nm; // 0.05 [nm] - constexpr double kGeomMinLength2 = kGeomMinLength * kGeomMinLength; // (0.05 [nm])^2 - if (dLength2 > kGeomMinLength2) { - const double dispR = std::sqrt(dLength2); - // Estimate safety by subtracting the geometrical step length. - currentTrack.safety -= currentTrack.geometryStepLength; - constexpr double sFact = 0.99; - double reducedSafety = sFact * currentTrack.safety; - - // Apply displacement, depending on how close we are to a boundary. - // 1a. Far away from geometry boundary: - if (reducedSafety > 0.0 && dispR <= reducedSafety) { - currentTrack.pos += displacement; - } else { - // Recompute safety. - currentTrack.safety = AdePTNavigator::ComputeSafety(currentTrack.pos, currentTrack.navState); - reducedSafety = sFact * currentTrack.safety; - - // 1b. Far away from geometry boundary: - if (reducedSafety > 0.0 && dispR <= reducedSafety) { - currentTrack.pos += displacement; - // 2. Push to boundary: - } else if (reducedSafety > kGeomMinLength) { - currentTrack.pos += displacement * (reducedSafety / dispR); - } - // 3. Very small safety: do nothing. - } - } - } - - // Collect the charged step length (might be changed by MSC). Collect the changes in energy and deposit. - currentTrack.eKin = theTrack->GetEKin(); - // double energyDeposit = theTrack->GetEnergyDeposit(); - - // Update the flight times of the particle - // By calculating the velocity here, we assume that all the energy deposit is done at the PreStepPoint, and - // the velocity depends on the remaining energy - double deltaTime = elTrack.GetPStepLength() / GetVelocity(currentTrack.eKin); - currentTrack.globalTime += deltaTime; - currentTrack.localTime += deltaTime; - currentTrack.properTime += deltaTime * (restMass / currentTrack.eKin); - } -} - -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// - -template -static __global__ void ElectronRelocation(adept::TrackManager *electrons) -{ - int activeSize = electrons->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*electrons->fActiveTracks)[i]; - Track ¤tTrack = (*electrons)[slot]; - - if (currentTrack.nextState.IsOnBoundary()) { - // if the particle hit a boundary, and is neither stopped or outside, relocate to have the correct next state - // before RecordHit is called - if (!currentTrack.stopped && !currentTrack.nextState.IsOutside()) { -#ifdef ADEPT_USE_SURF - AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.hitsurfID, currentTrack.nextState); -#else - AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.nextState); -#endif - } - } - } -} - -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// - -template -static __global__ void ElectronInteractions(adept::TrackManager *electrons, G4HepEmElectronTrack *hepEMTracks, - Secondaries secondaries, MParrayTracks *leakedQueue, - Scoring *userScoring, VolAuxData const *auxDataArray) -{ - constexpr Precision kPushOutRegion = 10 * vecgeom::kTolerance; - constexpr int Pdg = IsElectron ? 11 : -11; - fieldPropagatorConstBz fieldPropagatorBz(BzFieldValue); - - int activeSize = electrons->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*electrons->fActiveTracks)[i]; - Track ¤tTrack = (*electrons)[slot]; - adeptint::TrackData trackdata; - // the MCC vector is indexed by the logical volume id -#ifndef ADEPT_USE_SURF - const int lvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -#else - const int lvolID = currentTrack.navState.GetLogicalId(); -#endif - VolAuxData const &auxData = auxDataArray[lvolID]; - - auto survive = [&](bool leak = false) { - currentTrack.CopyTo(trackdata, Pdg); - if (leak) - leakedQueue->push_back(trackdata); - else - electrons->fNextTracks->push_back(slot); - }; - - // Init a track with the needed data to call into G4HepEm. - G4HepEmElectronTrack &elTrack = hepEMTracks[slot]; - G4HepEmTrack *theTrack = elTrack.GetTrack(); - G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData(); - - // // Prepare a branched RNG state while threads are synchronized. Even if not - // // used, this provides a fresh round of random numbers and reduces thread - // // divergence because the RNG state doesn't need to be advanced later. - // RanluxppDouble newRNG(currentTrack.rngState.BranchNoAdvance()); - - G4HepEmRandomEngine rnge(¤tTrack.rngState); - - if (auxData.fSensIndex >= 0) - adept_scoring::RecordHit(userScoring, currentTrack.parentID, - IsElectron ? 0 : 1, // Particle type - elTrack.GetPStepLength(), // Step length - theTrack->GetEnergyDeposit(), // Total Edep - ¤tTrack.preStepNavState, // Pre-step point navstate - ¤tTrack.preStepPos, // Pre-step point position - ¤tTrack.preStepDir, // Pre-step point momentum direction - nullptr, // Pre-step point polarization - currentTrack.preStepEKin, // Pre-step point kinetic energy - IsElectron ? -1 : 1, // Pre-step point charge - ¤tTrack.navState, // Post-step point navstate - ¤tTrack.pos, // Post-step point position - ¤tTrack.dir, // Post-step point momentum direction - nullptr, // Post-step point polarization - currentTrack.eKin, // Post-step point kinetic energy - IsElectron ? -1 : 1); // Post-step point charge - - // Save the `number-of-interaction-left` in our track. - for (int ip = 0; ip < 3; ++ip) { - double numIALeft = theTrack->GetNumIALeft(ip); - currentTrack.numIALeft[ip] = numIALeft; - } - - if (currentTrack.stopped) { - if (!IsElectron) { - // Annihilate the stopped positron into two gammas heading to opposite - // directions (isotropic). - Track &gamma1 = secondaries.gammas->NextTrack(); - Track &gamma2 = secondaries.gammas->NextTrack(); - - adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 0, /*numPositrons*/ 0, /*numGammas*/ 2); - - const double cost = 2 * currentTrack.Uniform() - 1; - const double sint = sqrt(1 - cost * cost); - const double phi = k2Pi * currentTrack.Uniform(); - double sinPhi, cosPhi; - sincos(phi, &sinPhi, &cosPhi); - - gamma1.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - currentTrack.newRNG.Advance(); - gamma1.parentID = currentTrack.parentID; - gamma1.rngState = currentTrack.newRNG; - gamma1.eKin = copcore::units::kElectronMassC2; - gamma1.dir.Set(sint * cosPhi, sint * sinPhi, cost); - - gamma2.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - // Reuse the RNG state of the dying track. - gamma2.parentID = currentTrack.parentID; - gamma2.rngState = currentTrack.rngState; - gamma2.eKin = copcore::units::kElectronMassC2; - gamma2.dir = -gamma1.dir; - } - // Particles are killed by not enqueuing them into the new activeQueue. - continue; - } - - if (currentTrack.nextState.IsOnBoundary()) { - // For now, just count that we hit something. - - // Kill the particle if it left the world. - if (!currentTrack.nextState.IsOutside()) { - - // Move to the next boundary. - currentTrack.navState = currentTrack.nextState; - // Check if the next volume belongs to the GPU region and push it to the appropriate queue -#ifndef ADEPT_USE_SURF - const int nextlvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -#else - const int nextlvolID = currentTrack.navState.GetLogicalId(); -#endif - VolAuxData const &nextauxData = auxDataArray[nextlvolID]; - if (nextauxData.fGPUregion > 0) - survive(); - else { - // To be safe, just push a bit the track exiting the GPU region to make sure - // Geant4 does not relocate it again inside the same region - currentTrack.pos += kPushOutRegion * currentTrack.dir; - survive(/*leak*/ true); - } - } - continue; - } else if (!currentTrack.propagated || currentTrack.restrictedPhysicalStepLength) { - // Did not yet reach the interaction point due to error in the magnetic - // field propagation. Try again next time. - survive(); - continue; - } else if (theTrack->GetWinnerProcessIndex() < 0) { - // No discrete process, move on. - survive(); - continue; - } - - // Reset number of interaction left for the winner discrete process. - // (Will be resampled in the next iteration.) - currentTrack.numIALeft[theTrack->GetWinnerProcessIndex()] = -1.0; - - // 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. - survive(); - continue; - } - - // Perform the discrete interaction, make sure the branched RNG state is - // ready to be used. - currentTrack.newRNG.Advance(); - // Also advance the current RNG state to provide a fresh round of random - // numbers after MSC used up a fair share for sampling the displacement. - currentTrack.rngState.Advance(); - - const double theElCut = g4HepEmData.fTheMatCutData->fMatCutData[auxData.fMCIndex].fSecElProdCutE; - - switch (theTrack->GetWinnerProcessIndex()) { - case 0: { - // Invoke ionization (for e-/e+): - double deltaEkin = - (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, currentTrack.eKin, &rnge) - : G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, currentTrack.eKin, &rnge); - - double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()}; - double dirSecondary[3]; - G4HepEmElectronInteractionIoni::SampleDirections(currentTrack.eKin, deltaEkin, dirSecondary, dirPrimary, &rnge); - - Track &secondary = secondaries.electrons->NextTrack(); - - adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 1, /*numPositrons*/ 0, /*numGammas*/ 0); - - secondary.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - secondary.parentID = currentTrack.parentID; - secondary.rngState = currentTrack.newRNG; - secondary.eKin = deltaEkin; - secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); - - currentTrack.eKin -= deltaEkin; - currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]); - survive(); - break; - } - case 1: { - // Invoke model for Bremsstrahlung: either SB- or Rel-Brem. - double logEnergy = std::log(currentTrack.eKin); - double deltaEkin = currentTrack.eKin < g4HepEmPars.fElectronBremModelLim - ? G4HepEmElectronInteractionBrem::SampleETransferSB( - &g4HepEmData, currentTrack.eKin, logEnergy, auxData.fMCIndex, &rnge, IsElectron) - : G4HepEmElectronInteractionBrem::SampleETransferRB( - &g4HepEmData, currentTrack.eKin, logEnergy, auxData.fMCIndex, &rnge, IsElectron); - - double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()}; - double dirSecondary[3]; - G4HepEmElectronInteractionBrem::SampleDirections(currentTrack.eKin, deltaEkin, dirSecondary, dirPrimary, &rnge); - - Track &gamma = secondaries.gammas->NextTrack(); - adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 0, /*numPositrons*/ 0, /*numGammas*/ 1); - - gamma.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - gamma.parentID = currentTrack.parentID; - gamma.rngState = currentTrack.newRNG; - gamma.eKin = deltaEkin; - gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); - - currentTrack.eKin -= deltaEkin; - currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]); - survive(); - break; - } - case 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( - currentTrack.eKin, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge); - - Track &gamma1 = secondaries.gammas->NextTrack(); - Track &gamma2 = secondaries.gammas->NextTrack(); - adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 0, /*numPositrons*/ 0, /*numGammas*/ 2); - - gamma1.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - gamma1.parentID = currentTrack.parentID; - gamma1.rngState = currentTrack.newRNG; - gamma1.eKin = theGamma1Ekin; - gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]); - - gamma2.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - // Reuse the RNG state of the dying track. - gamma2.parentID = currentTrack.parentID; - gamma2.rngState = currentTrack.rngState; - gamma2.eKin = theGamma2Ekin; - gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]); - - // The current track is killed by not enqueuing into the next activeQueue. - break; - } - } - } -} - -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// -////////////////////////////////////////////////////////////////////////////////// - -// template -// static __global__ void TransportElectrons(adept::TrackManager *electrons, G4HepEmElectronTrack *hepEMTracks, -// Secondaries secondaries, MParrayTracks *leakedQueue, -// Scoring *userScoring, VolAuxData const *auxDataArray) -// { -// #ifdef VECGEOM_FLOAT_PRECISION -// const Precision kPush = 10 * vecgeom::kTolerance; -// #else -// const Precision kPush = 0.; -// #endif -// constexpr Precision kPushOutRegion = 10 * vecgeom::kTolerance; -// constexpr int Charge = IsElectron ? -1 : 1; -// constexpr double restMass = copcore::units::kElectronMassC2; -// constexpr int Pdg = IsElectron ? 11 : -11; -// fieldPropagatorConstBz fieldPropagatorBz(BzFieldValue); - -// int activeSize = electrons->fActiveTracks->size(); -// for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { -// const int slot = (*electrons->fActiveTracks)[i]; -// Track ¤tTrack = (*electrons)[slot]; -// adeptint::TrackData trackdata; -// // the MCC vector is indexed by the logical volume id -// #ifndef ADEPT_USE_SURF -// const int lvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -// #else -// const int lvolID = currentTrack.navState.GetLogicalId(); -// #endif -// VolAuxData const &auxData = auxDataArray[lvolID]; - -// auto survive = [&](bool leak = false) { -// currentTrack.CopyTo(trackdata, Pdg); -// if (leak) -// leakedQueue->push_back(trackdata); -// else -// electrons->fNextTracks->push_back(slot); -// }; - -// // Init a track with the needed data to call into G4HepEm. -// G4HepEmElectronTrack &elTrack = hepEMTracks[slot]; -// G4HepEmTrack *theTrack = elTrack.GetTrack(); -// G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData(); - -// // // Prepare a branched RNG state while threads are synchronized. Even if not -// // // used, this provides a fresh round of random numbers and reduces thread -// // // divergence because the RNG state doesn't need to be advanced later. -// // RanluxppDouble newRNG(currentTrack.rngState.BranchNoAdvance()); - -// G4HepEmRandomEngine rnge(¤tTrack.rngState); - - // Check if there's a volume boundary in between. -// currentTrack.propagated = true; -// currentTrack.hitsurfID = -1; -// // double geometryStepLength; -// // vecgeom::NavigationState nextState; -// if (BzFieldValue != 0) { -// currentTrack.geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume( -// currentTrack.eKin, restMass, Charge, theTrack->GetGStepLength(), currentTrack.pos, currentTrack.dir, -// currentTrack.navState, currentTrack.nextState, currentTrack.hitsurfID, currentTrack.propagated, currentTrack.safety); -// } else { -// #ifdef ADEPT_USE_SURF -// currentTrack.geometryStepLength = -// AdePTNavigator::ComputeStepAndNextVolume(currentTrack.pos, currentTrack.dir, theTrack->GetGStepLength(), -// currentTrack.navState, currentTrack.nextState, currentTrack.hitsurfID, kPush); -// #else -// currentTrack.geometryStepLength = -// AdePTNavigator::ComputeStepAndNextVolume(currentTrack.pos, currentTrack.dir, theTrack->GetGStepLength(), -// currentTrack.navState, currentTrack.nextState, kPush); -// #endif -// currentTrack.pos += currentTrack.geometryStepLength * currentTrack.dir; -// } - -// // Set boundary state in navState so the next step and secondaries get the -// // correct information (navState = currentTrack.nextState only if relocated -// // in case of a boundary; see below) -// currentTrack.navState.SetBoundaryState(currentTrack.nextState.IsOnBoundary()); - -// // Propagate information from geometrical step to MSC. -// theTrack->SetDirection(currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()); -// theTrack->SetGStepLength(currentTrack.geometryStepLength); -// theTrack->SetOnBoundary(currentTrack.nextState.IsOnBoundary()); - - // Apply continuous effects. - // currentTrack.stopped = G4HepEmElectronManager::PerformContinuous(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge); - - // // Collect the direction change and displacement by MSC. - // const double *direction = theTrack->GetDirection(); - // currentTrack.dir.Set(direction[0], direction[1], direction[2]); - // if (!currentTrack.nextState.IsOnBoundary()) { - // const double *mscDisplacement = mscData->GetDisplacement(); - // vecgeom::Vector3D displacement(mscDisplacement[0], mscDisplacement[1], mscDisplacement[2]); - // const double dLength2 = displacement.Length2(); - // constexpr double kGeomMinLength = 5 * copcore::units::nm; // 0.05 [nm] - // constexpr double kGeomMinLength2 = kGeomMinLength * kGeomMinLength; // (0.05 [nm])^2 - // if (dLength2 > kGeomMinLength2) { - // const double dispR = std::sqrt(dLength2); - // // Estimate safety by subtracting the geometrical step length. - // currentTrack.safety -= currentTrack.geometryStepLength; - // constexpr double sFact = 0.99; - // double reducedSafety = sFact * currentTrack.safety; - - // // Apply displacement, depending on how close we are to a boundary. - // // 1a. Far away from geometry boundary: - // if (reducedSafety > 0.0 && dispR <= reducedSafety) { - // currentTrack.pos += displacement; - // } else { - // // Recompute safety. - // currentTrack.safety = AdePTNavigator::ComputeSafety(currentTrack.pos, currentTrack.navState); - // reducedSafety = sFact * currentTrack.safety; - - // // 1b. Far away from geometry boundary: - // if (reducedSafety > 0.0 && dispR <= reducedSafety) { - // currentTrack.pos += displacement; - // // 2. Push to boundary: - // } else if (reducedSafety > kGeomMinLength) { - // currentTrack.pos += displacement * (reducedSafety / dispR); - // } - // // 3. Very small safety: do nothing. - // } - // } - // } - - // // Collect the charged step length (might be changed by MSC). Collect the changes in energy and deposit. - // currentTrack.eKin = theTrack->GetEKin(); - // // double energyDeposit = theTrack->GetEnergyDeposit(); - - // // Update the flight times of the particle - // // By calculating the velocity here, we assume that all the energy deposit is done at the PreStepPoint, and - // // the velocity depends on the remaining energy - // double deltaTime = elTrack.GetPStepLength() / GetVelocity(currentTrack.eKin); - // currentTrack.globalTime += deltaTime; - // currentTrack.localTime += deltaTime; - // currentTrack.properTime += deltaTime * (restMass / currentTrack.eKin); - -// if (currentTrack.nextState.IsOnBoundary()) { -// // if the particle hit a boundary, and is neither stopped or outside, relocate to have the correct next state -// // before RecordHit is called -// if (!currentTrack.stopped && !currentTrack.nextState.IsOutside()) { -// #ifdef ADEPT_USE_SURF -// AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.hitsurfID, currentTrack.nextState); -// #else -// AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.nextState); -// #endif -// } -// } - -// if (auxData.fSensIndex >= 0) -// adept_scoring::RecordHit(userScoring, currentTrack.parentID, -// IsElectron ? 0 : 1, // Particle type -// elTrack.GetPStepLength(), // Step length -// theTrack->GetEnergyDeposit(), // Total Edep -// ¤tTrack.preStepNavState, // Pre-step point navstate -// ¤tTrack.preStepPos, // Pre-step point position -// ¤tTrack.preStepDir, // Pre-step point momentum direction -// nullptr, // Pre-step point polarization -// currentTrack.preStepEKin, // Pre-step point kinetic energy -// IsElectron ? -1 : 1, // Pre-step point charge -// ¤tTrack.navState, // Post-step point navstate -// ¤tTrack.pos, // Post-step point position -// ¤tTrack.dir, // Post-step point momentum direction -// nullptr, // Post-step point polarization -// currentTrack.eKin, // Post-step point kinetic energy -// IsElectron ? -1 : 1); // Post-step point charge - -// // Save the `number-of-interaction-left` in our track. -// for (int ip = 0; ip < 3; ++ip) { -// double numIALeft = theTrack->GetNumIALeft(ip); -// currentTrack.numIALeft[ip] = numIALeft; -// } - -// if (currentTrack.stopped) { -// if (!IsElectron) { -// // Annihilate the stopped positron into two gammas heading to opposite -// // directions (isotropic). -// Track &gamma1 = secondaries.gammas->NextTrack(); -// Track &gamma2 = secondaries.gammas->NextTrack(); - -// adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 0, /*numPositrons*/ 0, /*numGammas*/ 2); - -// const double cost = 2 * currentTrack.Uniform() - 1; -// const double sint = sqrt(1 - cost * cost); -// const double phi = k2Pi * currentTrack.Uniform(); -// double sinPhi, cosPhi; -// sincos(phi, &sinPhi, &cosPhi); - -// gamma1.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); -// currentTrack.newRNG.Advance(); -// gamma1.parentID = currentTrack.parentID; -// gamma1.rngState = currentTrack.newRNG; -// gamma1.eKin = copcore::units::kElectronMassC2; -// gamma1.dir.Set(sint * cosPhi, sint * sinPhi, cost); - -// gamma2.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); -// // Reuse the RNG state of the dying track. -// gamma2.parentID = currentTrack.parentID; -// gamma2.rngState = currentTrack.rngState; -// gamma2.eKin = copcore::units::kElectronMassC2; -// gamma2.dir = -gamma1.dir; -// } -// // Particles are killed by not enqueuing them into the new activeQueue. -// continue; -// } - -// if (currentTrack.nextState.IsOnBoundary()) { -// // For now, just count that we hit something. - -// // Kill the particle if it left the world. -// if (!currentTrack.nextState.IsOutside()) { - -// // Move to the next boundary. -// currentTrack.navState = currentTrack.nextState; -// // Check if the next volume belongs to the GPU region and push it to the appropriate queue -// #ifndef ADEPT_USE_SURF -// const int nextlvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -// #else -// const int nextlvolID = currentTrack.navState.GetLogicalId(); -// #endif -// VolAuxData const &nextauxData = auxDataArray[nextlvolID]; -// if (nextauxData.fGPUregion > 0) -// survive(); -// else { -// // To be safe, just push a bit the track exiting the GPU region to make sure -// // Geant4 does not relocate it again inside the same region -// currentTrack.pos += kPushOutRegion * currentTrack.dir; -// survive(/*leak*/ true); -// } -// } -// continue; -// } else if (!currentTrack.propagated || currentTrack.restrictedPhysicalStepLength) { -// // Did not yet reach the interaction point due to error in the magnetic -// // field propagation. Try again next time. -// survive(); -// continue; -// } else if (theTrack->GetWinnerProcessIndex() < 0) { -// // No discrete process, move on. -// survive(); -// continue; -// } - -// // Reset number of interaction left for the winner discrete process. -// // (Will be resampled in the next iteration.) -// currentTrack.numIALeft[theTrack->GetWinnerProcessIndex()] = -1.0; - -// // 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. -// survive(); -// continue; -// } - -// // Perform the discrete interaction, make sure the branched RNG state is -// // ready to be used. -// currentTrack.newRNG.Advance(); -// // Also advance the current RNG state to provide a fresh round of random -// // numbers after MSC used up a fair share for sampling the displacement. -// currentTrack.rngState.Advance(); - -// const double theElCut = g4HepEmData.fTheMatCutData->fMatCutData[auxData.fMCIndex].fSecElProdCutE; - -// switch (theTrack->GetWinnerProcessIndex()) { -// case 0: { -// // Invoke ionization (for e-/e+): -// double deltaEkin = -// (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, currentTrack.eKin, &rnge) -// : G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, currentTrack.eKin, &rnge); - -// double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()}; -// double dirSecondary[3]; -// G4HepEmElectronInteractionIoni::SampleDirections(currentTrack.eKin, deltaEkin, dirSecondary, dirPrimary, &rnge); - -// Track &secondary = secondaries.electrons->NextTrack(); - -// adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 1, /*numPositrons*/ 0, /*numGammas*/ 0); - -// secondary.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); -// secondary.parentID = currentTrack.parentID; -// secondary.rngState = currentTrack.newRNG; -// secondary.eKin = deltaEkin; -// secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); - -// currentTrack.eKin -= deltaEkin; -// currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]); -// survive(); -// break; -// } -// case 1: { -// // Invoke model for Bremsstrahlung: either SB- or Rel-Brem. -// double logEnergy = std::log(currentTrack.eKin); -// double deltaEkin = currentTrack.eKin < g4HepEmPars.fElectronBremModelLim -// ? G4HepEmElectronInteractionBrem::SampleETransferSB( -// &g4HepEmData, currentTrack.eKin, logEnergy, auxData.fMCIndex, &rnge, IsElectron) -// : G4HepEmElectronInteractionBrem::SampleETransferRB( -// &g4HepEmData, currentTrack.eKin, logEnergy, auxData.fMCIndex, &rnge, IsElectron); - -// double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()}; -// double dirSecondary[3]; -// G4HepEmElectronInteractionBrem::SampleDirections(currentTrack.eKin, deltaEkin, dirSecondary, dirPrimary, &rnge); - -// Track &gamma = secondaries.gammas->NextTrack(); -// adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 0, /*numPositrons*/ 0, /*numGammas*/ 1); - -// gamma.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); -// gamma.parentID = currentTrack.parentID; -// gamma.rngState = currentTrack.newRNG; -// gamma.eKin = deltaEkin; -// gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); - -// currentTrack.eKin -= deltaEkin; -// currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]); -// survive(); -// break; -// } -// case 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( -// currentTrack.eKin, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge); - -// Track &gamma1 = secondaries.gammas->NextTrack(); -// Track &gamma2 = secondaries.gammas->NextTrack(); -// adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 0, /*numPositrons*/ 0, /*numGammas*/ 2); - -// gamma1.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); -// gamma1.parentID = currentTrack.parentID; -// gamma1.rngState = currentTrack.newRNG; -// gamma1.eKin = theGamma1Ekin; -// gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]); - -// gamma2.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); -// // Reuse the RNG state of the dying track. -// gamma2.parentID = currentTrack.parentID; -// gamma2.rngState = currentTrack.rngState; -// gamma2.eKin = theGamma2Ekin; -// gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]); - -// // The current track is killed by not enqueuing into the next activeQueue. -// break; -// } -// } -// } -// } - - - - - - - - - - - -//////////////////////// - -// // Init G4HepEMTrack again -// G4HepEmElectronTrack elTrack; -// G4HepEmTrack *theTrack = elTrack.GetTrack(); -// theTrack->SetEKin(currentTrack.eKin); -// theTrack->SetMCIndex(auxData.fMCIndex); -// theTrack->SetOnBoundary(currentTrack.navState.IsOnBoundary()); -// theTrack->SetCharge(Charge); -// G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData(); -// mscData->fIsFirstStep = currentTrack.initialRange < 0; -// mscData->fInitialRange = currentTrack.initialRange; -// mscData->fDynamicRangeFactor = currentTrack.dynamicRangeFactor; -// mscData->fTlimitMin = currentTrack.tlimitMin; -// // And add the information we stored -// theTrack->SetGStepLength(currentTrack.geometryStepLength); -// elTrack.SetPStepLength(currentTrack.physicalStepLength); -// for (int i = 0; i < 3; i++) { -// theTrack->SetMFP(currentTrack.preStepMFPs[i], i); -// theTrack->SetNumIALeft(currentTrack.numIALeft[i], i); -// } - -//////////////////////// - -// Init rnge again -// G4HepEmRandomEngine rnge(¤tTrack.rngState); - -//////////////////////// - -// This is only needed for the interactions - -// Prepare a branched RNG state while threads are synchronized. Even if not -// used, this provides a fresh round of random numbers and reduces thread -// divergence because the RNG state doesn't need to be advanced later. -// RanluxppDouble newRNG(currentTrack.rngState.BranchNoAdvance()); - -//////////////////////// \ No newline at end of file diff --git a/include/AdePT/kernels/gammas_experimental.cuh b/include/AdePT/kernels/gammas_experimental.cuh deleted file mode 100644 index 9b85e095..00000000 --- a/include/AdePT/kernels/gammas_experimental.cuh +++ /dev/null @@ -1,411 +0,0 @@ -// SPDX-FileCopyrightText: 2022 CERN -// SPDX-License-Identifier: Apache-2.0 - -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -// Pull in implementation. -#include -#include -#include -#include - -using VolAuxData = adeptint::VolAuxData; - -__global__ void GammaPhysics1(adept::TrackManager *gammas, G4HepEmGammaTrack *hepEMTracks, VolAuxData const *auxDataArray) -{ - int activeSize = gammas->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*gammas->fActiveTracks)[i]; - Track ¤tTrack = (*gammas)[slot]; - // Save current values for scoring - currentTrack.preStepEKin = currentTrack.eKin; - currentTrack.preStepPos = currentTrack.pos; - currentTrack.preStepDir = currentTrack.dir; - currentTrack.preStepNavState = currentTrack.navState; - // the MCC vector is indexed by the logical volume id -#ifndef ADEPT_USE_SURF - int lvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -#else - int lvolID = currentTrack.navState.GetLogicalId(); -#endif - VolAuxData const &auxData = auxDataArray[lvolID]; - - // Init a track with the needed data to call into G4HepEm. - G4HepEmGammaTrack &gammaTrack = hepEMTracks[slot]; - // G4HepEmGammaTrack gammaTrack; - G4HepEmTrack *theTrack = gammaTrack.GetTrack(); - theTrack->SetEKin(currentTrack.eKin); - theTrack->SetMCIndex(auxData.fMCIndex); - - // Sample the `number-of-interaction-left` and put it into the track. - for (int ip = 0; ip < 3; ++ip) { - double numIALeft = currentTrack.numIALeft[ip]; - if (numIALeft <= 0) { - numIALeft = -std::log(currentTrack.Uniform()); - // currentTrack.numIALeft[ip] = numIALeft; - } - theTrack->SetNumIALeft(numIALeft, ip); - } - - // Call G4HepEm to compute the physics step limit. - G4HepEmGammaManager::HowFar(&g4HepEmData, &g4HepEmPars, &gammaTrack); - - // hepEMTracks[slot] = gammaTrack; - - // Save the info we need from the G4HepEM track - // currentTrack.geometricalStepLengthFromPhysics = theTrack->GetGStepLength(); - // currentTrack.winnerProcessIndex = theTrack->GetWinnerProcessIndex(); - // currentTrack.PEmxSec = gammaTrack.GetPEmxSec(); - // currentTrack.preStepMFPs[0] = theTrack->GetMFP(0); - // currentTrack.preStepMFPs[1] = theTrack->GetMFP(1); - // currentTrack.preStepMFPs[2] = theTrack->GetMFP(2); - } -} - -__global__ void GammaTransport1(adept::TrackManager *gammas, G4HepEmGammaTrack *hepEMTracks, VolAuxData const *auxDataArray) -{ -#ifdef VECGEOM_FLOAT_PRECISION - const Precision kPush = 10 * vecgeom::kTolerance; -#else - const Precision kPush = 0.; -#endif - // constexpr Precision kPushOutRegion = 10 * vecgeom::kTolerance; - // constexpr int Pdg = 22; - int activeSize = gammas->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*gammas->fActiveTracks)[i]; - Track ¤tTrack = (*gammas)[slot]; - // Retrieve the HepEM Track - G4HepEmTrack *theTrack = hepEMTracks[slot].GetTrack(); - - // Check if there's a volume boundary in between. - vecgeom::NavigationState nextState; - double geometryStepLength; -#ifdef ADEPT_USE_SURF - long hitsurf_index = -1; - geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume( - currentTrack.pos, currentTrack.dir, theTrack->GetGStepLength(), currentTrack.navState, - nextState, hitsurf_index, kPush); - currentTrack.hitsurfID = hitsurf_index; -#else - geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(currentTrack.pos, currentTrack.dir, - theTrack->GetGStepLength(), - currentTrack.navState, nextState, kPush); -#endif - currentTrack.pos += geometryStepLength * currentTrack.dir; - - // Store the actual geometrical step length traveled - currentTrack.geometryStepLength = geometryStepLength; - - // Set boundary state in navState so the next step and secondaries get the - // correct information (navState = nextState only if relocated - // in case of a boundary; see below) - currentTrack.navState.SetBoundaryState(nextState.IsOnBoundary()); - - // Propagate information from geometrical step to G4HepEm. - theTrack->SetGStepLength(geometryStepLength); - theTrack->SetOnBoundary(nextState.IsOnBoundary()); - - // Update the number-of-interaction-left - G4HepEmGammaManager::UpdateNumIALeft(theTrack); - - // Save the `number-of-interaction-left` in our track. - for (int ip = 0; ip < 3; ++ip) { - double numIALeft = theTrack->GetNumIALeft(ip); - currentTrack.numIALeft[ip] = numIALeft; - } - - // Update the flight times of the particle - double deltaTime = theTrack->GetGStepLength() / copcore::units::kCLight; - currentTrack.globalTime += deltaTime; - currentTrack.localTime += deltaTime; - - // Save the next state in the track - currentTrack.nextState = nextState; - } -} - -__global__ void GammaRelocation(adept::TrackManager *gammas, MParrayTracks *leakedQueue, - VolAuxData const *auxDataArray) -{ - // #ifdef VECGEOM_FLOAT_PRECISION - // const Precision kPush = 10 * vecgeom::kTolerance; - // #else - // const Precision kPush = 0.; - // #endif - constexpr Precision kPushOutRegion = 10 * vecgeom::kTolerance; - constexpr int Pdg = 22; - int activeSize = gammas->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*gammas->fActiveTracks)[i]; - Track ¤tTrack = (*gammas)[slot]; - - auto survive = [&](bool leak = false) { - adeptint::TrackData trackdata; - currentTrack.CopyTo(trackdata, Pdg); - if (leak) - leakedQueue->push_back(trackdata); - else - gammas->fNextTracks->push_back(slot); - }; - - // Last transportation call, in case we are on a boundary we need to relocate to the next volume - // Every track on a boundary will stop being processed here (no interaction) - // For now (experimental kernels) I just mark it in the track, and in the next kernel we just - // return inmediately (but ideally we won't even give those tracks to the kernel) - if (currentTrack.nextState.IsOnBoundary()) { - currentTrack.restrictedPhysicalStepLength = true; - // Kill the particle if it left the world. - if (!currentTrack.nextState.IsOutside()) { -#ifdef ADEPT_USE_SURF - AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.hitsurfID, - currentTrack.nextState); - if (currentTrack.nextState.IsOutside()) continue; -#else - AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.nextState); -#endif - // Move to the next boundary. - currentTrack.navState = currentTrack.nextState; - // Check if the next volume belongs to the GPU region and push it to the appropriate queue -#ifndef ADEPT_USE_SURF - const int nextlvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -#else - const int nextlvolID = currentTrack.navState.GetLogicalId(); -#endif - VolAuxData const &nextauxData = auxDataArray[nextlvolID]; - if (nextauxData.fGPUregion > 0) - survive(); - else { - // To be safe, just push a bit the track exiting the GPU region to make sure - // Geant4 does not relocate it again inside the same region - currentTrack.pos += kPushOutRegion * currentTrack.dir; - survive(/*leak*/ true); - } - } - continue; - } else { - currentTrack.restrictedPhysicalStepLength = false; - } - } -} - -template -__global__ void GammaPhysics2(adept::TrackManager *gammas, G4HepEmGammaTrack *hepEMTracks, Secondaries secondaries, Scoring *userScoring, - VolAuxData const *auxDataArray) -{ - int activeSize = gammas->fActiveTracks->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*gammas->fActiveTracks)[i]; - Track ¤tTrack = (*gammas)[slot]; - // Retrieve the HepEM Track - G4HepEmGammaTrack &gammaTrack = hepEMTracks[slot]; - G4HepEmTrack *theTrack = gammaTrack.GetTrack(); - -#ifndef ADEPT_USE_SURF - int lvolID = currentTrack.navState.Top()->GetLogicalVolume()->id(); -#else - int lvolID = currentTrack.navState.GetLogicalId(); -#endif - VolAuxData const &auxData = auxDataArray[lvolID]; - - auto survive = [&]() { gammas->fNextTracks->push_back(slot); }; - - // Temporary solution - if (currentTrack.restrictedPhysicalStepLength) { - continue; - } - - if(theTrack->GetWinnerProcessIndex() < 0) - { - continue; - } - - // Reset number of interaction left for the winner discrete process. - // (Will be resampled in the next iteration.) - currentTrack.numIALeft[theTrack->GetWinnerProcessIndex()] = -1.0; - - // Perform the discrete interaction. - G4HepEmRandomEngine rnge(¤tTrack.rngState); - // We might need one branched RNG state, prepare while threads are synchronized. - RanluxppDouble newRNG(currentTrack.rngState.Branch()); - - switch (theTrack->GetWinnerProcessIndex()) { - case 0: { - // Invoke gamma conversion to e-/e+ pairs, if the energy is above the threshold. - if (currentTrack.eKin < 2 * copcore::units::kElectronMassC2) { - survive(); - continue; - } - - double logEnergy = std::log(currentTrack.eKin); - double elKinEnergy, posKinEnergy; - G4HepEmGammaInteractionConversion::SampleKinEnergies(&g4HepEmData, currentTrack.eKin, logEnergy, auxData.fMCIndex, - elKinEnergy, posKinEnergy, &rnge); - - double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()}; - double dirSecondaryEl[3], dirSecondaryPos[3]; - G4HepEmGammaInteractionConversion::SampleDirections(dirPrimary, dirSecondaryEl, dirSecondaryPos, elKinEnergy, - posKinEnergy, &rnge); - - Track &electron = secondaries.electrons->NextTrack(); - Track &positron = secondaries.positrons->NextTrack(); - - adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 1, /*numPositrons*/ 1, /*numGammas*/ 0); - - electron.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - electron.parentID = currentTrack.parentID; - electron.rngState = newRNG; - electron.eKin = elKinEnergy; - electron.dir.Set(dirSecondaryEl[0], dirSecondaryEl[1], dirSecondaryEl[2]); - - positron.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - // Reuse the RNG state of the dying track. - positron.parentID = currentTrack.parentID; - positron.rngState = currentTrack.rngState; - positron.eKin = posKinEnergy; - positron.dir.Set(dirSecondaryPos[0], dirSecondaryPos[1], dirSecondaryPos[2]); - - // The current track is killed by not enqueuing into the next activeQueue. - break; - } - case 1: { - // Invoke Compton scattering of gamma. - constexpr double LowEnergyThreshold = 100 * copcore::units::eV; - if (currentTrack.eKin < LowEnergyThreshold) { - survive(); - continue; - } - const double origDirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()}; - double dirPrimary[3]; - const double newEnergyGamma = G4HepEmGammaInteractionCompton::SamplePhotonEnergyAndDirection( - currentTrack.eKin, dirPrimary, origDirPrimary, &rnge); - vecgeom::Vector3D newDirGamma(dirPrimary[0], dirPrimary[1], dirPrimary[2]); - - const double energyEl = currentTrack.eKin - newEnergyGamma; - if (energyEl > LowEnergyThreshold) { - // Create a secondary electron and sample/compute directions. - Track &electron = secondaries.electrons->NextTrack(); - adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 1, /*numPositrons*/ 0, /*numGammas*/ 0); - - electron.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - electron.parentID = currentTrack.parentID; - electron.rngState = newRNG; - electron.eKin = energyEl; - electron.dir = currentTrack.eKin * currentTrack.dir - newEnergyGamma * newDirGamma; - electron.dir.Normalize(); - } else { - if (auxData.fSensIndex >= 0) - adept_scoring::RecordHit(userScoring, - currentTrack.parentID, // Track ID - 2, // Particle type - currentTrack.geometryStepLength, // Step length - 0, // Total Edep - ¤tTrack.preStepNavState, // Pre-step point navstate - ¤tTrack.preStepPos, // Pre-step point position - ¤tTrack.preStepDir, // Pre-step point momentum direction - nullptr, // Pre-step point polarization - currentTrack.preStepEKin, // Pre-step point kinetic energy - 0, // Pre-step point charge - ¤tTrack.navState, // Post-step point navstate - ¤tTrack.pos, // Post-step point position - ¤tTrack.dir, // Post-step point momentum direction - nullptr, // Post-step point polarization - newEnergyGamma, // Post-step point kinetic energy - 0); // Post-step point charge - } - - // Check the new gamma energy and deposit if below threshold. - if (newEnergyGamma > LowEnergyThreshold) { - currentTrack.eKin = newEnergyGamma; - currentTrack.dir = newDirGamma; - survive(); - } else { - if (auxData.fSensIndex >= 0) - adept_scoring::RecordHit(userScoring, - currentTrack.parentID, // Track ID - 2, // Particle type - currentTrack.geometryStepLength, // Step length - 0, // Total Edep - ¤tTrack.preStepNavState, // Pre-step point navstate - ¤tTrack.preStepPos, // Pre-step point position - ¤tTrack.preStepDir, // Pre-step point momentum direction - nullptr, // Pre-step point polarization - currentTrack.preStepEKin, // Pre-step point kinetic energy - 0, // Pre-step point charge - ¤tTrack.navState, // Post-step point navstate - ¤tTrack.pos, // Post-step point position - ¤tTrack.dir, // Post-step point momentum direction - nullptr, // Post-step point polarization - newEnergyGamma, // Post-step point kinetic energy - 0); // Post-step point charge - // The current track is killed by not enqueuing into the next activeQueue. - } - break; - } - case 2: { - // Invoke photoelectric process. - const double theLowEnergyThreshold = 1 * copcore::units::eV; - - ////////////////////////////////// - - // const double bindingEnergy = G4HepEmGammaInteractionPhotoelectric::SelectElementBindingEnergy( - // &g4HepEmData, auxData.fMCIndex, currentTrack.PEmxSec, currentTrack.eKin, &rnge); - const double bindingEnergy = G4HepEmGammaInteractionPhotoelectric::SelectElementBindingEnergy( - &g4HepEmData, auxData.fMCIndex, gammaTrack.GetPEmxSec(), currentTrack.eKin, &rnge); - double edep = bindingEnergy; - - const double photoElecE = currentTrack.eKin - edep; - if (photoElecE > theLowEnergyThreshold) { - // Create a secondary electron and sample directions. - Track &electron = secondaries.electrons->NextTrack(); - adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 1, /*numPositrons*/ 0, /*numGammas*/ 0); - - double dirGamma[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()}; - double dirPhotoElec[3]; - G4HepEmGammaInteractionPhotoelectric::SamplePhotoElectronDirection(photoElecE, dirGamma, dirPhotoElec, &rnge); - - electron.InitAsSecondary(currentTrack.pos, currentTrack.navState, currentTrack.globalTime); - electron.parentID = currentTrack.parentID; - electron.rngState = newRNG; - electron.eKin = photoElecE; - electron.dir.Set(dirPhotoElec[0], dirPhotoElec[1], dirPhotoElec[2]); - } else { - edep = currentTrack.eKin; - } - - if (auxData.fSensIndex >= 0) - adept_scoring::RecordHit(userScoring, - currentTrack.parentID, // Track ID - 2, // Particle type - currentTrack.geometryStepLength, // Step length - edep, // Total Edep - ¤tTrack.preStepNavState, // Pre-step point navstate - ¤tTrack.preStepPos, // Pre-step point position - ¤tTrack.preStepDir, // Pre-step point momentum direction - nullptr, // Pre-step point polarization - currentTrack.preStepEKin, // Pre-step point kinetic energy - 0, // Pre-step point charge - ¤tTrack.navState, // Post-step point navstate - ¤tTrack.pos, // Post-step point position - ¤tTrack.dir, // Post-step point momentum direction - nullptr, // Post-step point polarization - 0, // Post-step point kinetic energy - 0); // Post-step point charge - - ////////////////////////////////// - - // The current track is killed by not enqueuing into the next activeQueue. - break; - } - } - } -}