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/examples/Example1/CMakeLists.txt b/examples/Example1/CMakeLists.txt index adc85a54..595d455f 100644 --- a/examples/Example1/CMakeLists.txt +++ b/examples/Example1/CMakeLists.txt @@ -74,8 +74,6 @@ configure_file("macros/example1_ttbar_noadept.mac.in" "${PROJECT_BINARY_DIR}/exa # Tests -string(APPEND CMAKE_CUDA_FLAGS " -Xptxas=-v") - add_test(NAME example1 COMMAND $ -m ${PROJECT_BINARY_DIR}/example1_large_stack.mac ) \ No newline at end of file diff --git a/examples/Example1/macros/example1_ttbar_LHCb.mac.in b/examples/Example1/macros/example1_ttbar_LHCb.mac.in index 3c1497a2..051736cc 100644 --- a/examples/Example1/macros/example1_ttbar_LHCb.mac.in +++ b/examples/Example1/macros/example1_ttbar_LHCb.mac.in @@ -23,10 +23,10 @@ ## Threshold for buffering tracks before sending to GPU /adept/setTransportBufferThreshold 2000 ## Total number of GPU track slots (not per thread) -/adept/setMillionsOfTrackSlots 8 +/adept/setMillionsOfTrackSlots 4 /adept/setMillionsOfHitSlots 1 ## Device stack limit -/adept/setCUDAStackLimit 4096 +# /adept/setCUDAStackLimit 4096 # If true, particles are transported on the GPU across the whole geometry, GPU regions are ignored /adept/setTrackInAllRegions true diff --git a/include/AdePT/core/AdePTTransport.cuh b/include/AdePT/core/AdePTTransport.cuh index 3318fe29..a06326c3 100644 --- a/include/AdePT/core/AdePTTransport.cuh +++ b/include/AdePT/core/AdePTTransport.cuh @@ -9,12 +9,14 @@ #include #include #include -#include -// #include - - -#include +#ifndef USE_SPLIT_KERNELS +#include +#include +#else +#include +#include +#endif #include #ifdef VECGEOM_ENABLE_CUDA @@ -34,6 +36,8 @@ #include #include #include +#include +#include #include #include @@ -263,6 +267,16 @@ GPUstate *InitializeGPU(adeptint::TrackBuffer &buffer, int capacity, int maxbatc COPCORE_CUDA_CHECK(cudaStreamCreate(&gpuState.particles[i].stream)); 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()); @@ -289,6 +303,12 @@ void FreeGPU(GPUstate &gpuState, G4HepEmState *g4hepem_state) COPCORE_CUDA_CHECK(cudaFreeHost(gpuState.stats)); COPCORE_CUDA_CHECK(cudaFree(gpuState.toDevice_dev)); +#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)); for (int i = 0; i < ParticleType::NumParticleTypes; i++) { @@ -381,10 +401,22 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & transportBlocks = (numElectrons + TransportThreads - 1) / TransportThreads; transportBlocks = std::min(transportBlocks, MaxBlocks); #endif +#ifndef USE_SPLIT_KERNELS TransportElectrons<<>>( electrons.trackmgr, secondaries, electrons.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); - +#else + ElectronHowFar<<>>( + electrons.trackmgr, gpuState.hepEMBuffers_d.electronsHepEm, VolAuxArray::GetInstance().fAuxData_dev); + ElectronPropagation<<>>( + electrons.trackmgr, gpuState.hepEMBuffers_d.electronsHepEm); + ElectronMSC<<>>( + electrons.trackmgr, gpuState.hepEMBuffers_d.electronsHepEm); + ElectronRelocation<<>>(electrons.trackmgr); + 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)); } @@ -396,10 +428,22 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & transportBlocks = (numPositrons + TransportThreads - 1) / TransportThreads; transportBlocks = std::min(transportBlocks, MaxBlocks); #endif +#ifndef USE_SPLIT_KERNELS TransportPositrons<<>>( positrons.trackmgr, secondaries, positrons.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); - +#else + ElectronHowFar<<>>( + positrons.trackmgr, gpuState.hepEMBuffers_d.positronsHepEm, VolAuxArray::GetInstance().fAuxData_dev); + ElectronPropagation<<>>( + positrons.trackmgr, gpuState.hepEMBuffers_d.positronsHepEm); + ElectronMSC<<>>( + positrons.trackmgr, gpuState.hepEMBuffers_d.positronsHepEm); + ElectronRelocation<<>>(positrons.trackmgr); + 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)); } @@ -411,22 +455,20 @@ void ShowerGPU(IntegrationLayer &integration, int event, adeptint::TrackBuffer & transportBlocks = (numGammas + TransportThreads - 1) / TransportThreads; transportBlocks = std::min(transportBlocks, MaxBlocks); #endif - Physics1<<>>( - gammas.trackmgr, VolAuxArray::GetInstance().fAuxData_dev - ); - Transport1<<>>( - gammas.trackmgr, VolAuxArray::GetInstance().fAuxData_dev - ); - Relocation<<>>( - gammas.trackmgr, gammas.leakedTracks, VolAuxArray::GetInstance().fAuxData_dev - ); - Physics2<<>>( - gammas.trackmgr, secondaries, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev - ); - - // TransportGammas<<>>( - // gammas.trackmgr, secondaries, gammas.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); - +#ifndef USE_SPLIT_KERNELS + TransportGammas<<>>( + gammas.trackmgr, secondaries, gammas.leakedTracks, scoring_dev, VolAuxArray::GetInstance().fAuxData_dev); +#else + GammaHowFar<<>>( + gammas.trackmgr, gpuState.hepEMBuffers_d.gammasHepEm, VolAuxArray::GetInstance().fAuxData_dev); + GammaPropagation<<>>( + gammas.trackmgr, gpuState.hepEMBuffers_d.gammasHepEm, VolAuxArray::GetInstance().fAuxData_dev); + GammaRelocation<<>>(gammas.trackmgr, gammas.leakedTracks, + VolAuxArray::GetInstance().fAuxData_dev); + GammaInteractions<<>>( + gammas.trackmgr, gpuState.hepEMBuffers_d.gammasHepEm, secondaries, 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 bdd5a5e5..ff80e9f7 100644 --- a/include/AdePT/core/AdePTTransportStruct.cuh +++ b/include/AdePT/core/AdePTTransportStruct.cuh @@ -14,6 +14,11 @@ #include #include +#ifdef USE_SPLIT_KERNELS +#include +#include +#endif + #ifdef __CUDA_ARCH__ // Define inline implementations of the RNG methods for the device. // (nvcc ignores the __device__ attribute in definitions, so this is only to @@ -65,6 +70,14 @@ 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 { adept::TrackManager::Stats mgr_stats[ParticleType::NumParticleTypes]; @@ -78,6 +91,9 @@ struct GPUstate { ParticleType particles[ParticleType::NumParticleTypes]; AllTrackManagers allmgr_h; ///< Host pointers for track managers AllTrackManagers allmgr_d; ///< Device pointers for track managers +#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 49a9251d..f41d6f0b 100644 --- a/include/AdePT/core/Track.cuh +++ b/include/AdePT/core/Track.cuh @@ -8,6 +8,8 @@ #include #include +#include + #include #include @@ -17,11 +19,9 @@ struct Track { using Precision = vecgeom::Precision; int parentID{0}; // Stores the track id of the initial particle given to AdePT - int hitsurfID{0}; RanluxppDouble rngState; double eKin; - double preStepEKin; double numIALeft[3]; double initialRange; double dynamicRangeFactor; @@ -32,22 +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 - bool reachedInteractionPoint{false}; + long hitsurfID{0}; + bool propagated{false}; double geometryStepLength{0}; + double safety{0}; - // Variables used to store physics results from G4HepEM - double geometricalStepLengthFromPhysics{0}; - int winnerProcessIndex{0}; - double preStepMFPs[3]; - double PEmxSec{0}; + // 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_split.cuh b/include/AdePT/kernels/electrons_split.cuh new file mode 100644 index 00000000..c52f4981 --- /dev/null +++ b/include/AdePT/kernels/electrons_split.cuh @@ -0,0 +1,512 @@ +// 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 ElectronHowFar(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()); + } + 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 ElectronPropagation(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 ElectronMSC(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(); + + // 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(); + + 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; + } + } + } +} diff --git a/include/AdePT/kernels/gammas_experimental.cuh b/include/AdePT/kernels/gammas_split.cuh similarity index 66% rename from include/AdePT/kernels/gammas_experimental.cuh rename to include/AdePT/kernels/gammas_split.cuh index e89f6697..1937d56e 100644 --- a/include/AdePT/kernels/gammas_experimental.cuh +++ b/include/AdePT/kernels/gammas_split.cuh @@ -18,17 +18,20 @@ #include #include -__global__ void Physics1(adept::TrackManager *gammas, VolAuxData const *auxDataArray) +using VolAuxData = adeptint::VolAuxData; + +__global__ void GammaHowFar(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; + 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(); @@ -38,8 +41,8 @@ __global__ void Physics1(adept::TrackManager *gammas, VolAuxData const *a VolAuxData const &auxData = auxDataArray[lvolID]; // Init a track with the needed data to call into G4HepEm. - G4HepEmGammaTrack gammaTrack; - G4HepEmTrack *theTrack = gammaTrack.GetTrack(); + G4HepEmGammaTrack &gammaTrack = hepEMTracks[slot]; + G4HepEmTrack *theTrack = gammaTrack.GetTrack(); theTrack->SetEKin(currentTrack.eKin); theTrack->SetMCIndex(auxData.fMCIndex); @@ -54,50 +57,36 @@ __global__ void Physics1(adept::TrackManager *gammas, VolAuxData const *a // Call G4HepEm to compute the physics step limit. G4HepEmGammaManager::HowFar(&g4HepEmData, &g4HepEmPars, &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 Transport1(adept::TrackManager *gammas, VolAuxData const *auxDataArray) +__global__ void GammaPropagation(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(); + 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); - // }; + // 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, currentTrack.geometricalStepLengthFromPhysics, currentTrack.navState, - nextState, hitsurf_index, kPush); + 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, currentTrack.geometricalStepLengthFromPhysics, currentTrack.navState, - nextState, kPush); + geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume( + currentTrack.pos, currentTrack.dir, theTrack->GetGStepLength(), currentTrack.navState, nextState, kPush); #endif currentTrack.pos += geometryStepLength * currentTrack.dir; @@ -109,41 +98,32 @@ __global__ void Transport1(adept::TrackManager *gammas, VolAuxData const // in case of a boundary; see below) currentTrack.navState.SetBoundaryState(nextState.IsOnBoundary()); - // Now init a G4HepEM track with the info it needs to update the the number-of-interaction-left - G4HepEmTrack theTrack; - theTrack.SetGStepLength(geometryStepLength); - theTrack.SetOnBoundary(nextState.IsOnBoundary()); - for (int i = 0; i < 3; i++) { - theTrack.SetMFP(currentTrack.preStepMFPs[i], i); - theTrack.SetNumIALeft(currentTrack.numIALeft[i], i); - } + // Propagate information from geometrical step to G4HepEm. + theTrack->SetGStepLength(geometryStepLength); + theTrack->SetOnBoundary(nextState.IsOnBoundary()); + // Update the number-of-interaction-left - G4HepEmGammaManager::UpdateNumIALeft(&theTrack); + G4HepEmGammaManager::UpdateNumIALeft(theTrack); // Save the `number-of-interaction-left` in our track. for (int ip = 0; ip < 3; ++ip) { - double numIALeft = theTrack.GetNumIALeft(ip); + double numIALeft = theTrack->GetNumIALeft(ip); currentTrack.numIALeft[ip] = numIALeft; } // Update the flight times of the particle - double deltaTime = theTrack.GetGStepLength() / copcore::units::kCLight; + double deltaTime = theTrack->GetGStepLength() / copcore::units::kCLight; currentTrack.globalTime += deltaTime; currentTrack.localTime += deltaTime; - - // Save the hitsurf index and the next state in the track - currentTrack.hitsurfID = hitsurf_index; + + // Save the next state in the track currentTrack.nextState = nextState; } } -__global__ void Relocation(adept::TrackManager *gammas, MParrayTracks *leakedQueue, VolAuxData const *auxDataArray) +__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(); @@ -165,11 +145,12 @@ __global__ void Relocation(adept::TrackManager *gammas, MParrayTracks *le // 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.reachedInteractionPoint = false; + 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); + AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.hitsurfID, + currentTrack.nextState); if (currentTrack.nextState.IsOutside()) continue; #else AdePTNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, currentTrack.nextState); @@ -194,19 +175,22 @@ __global__ void Relocation(adept::TrackManager *gammas, MParrayTracks *le } continue; } else { - currentTrack.reachedInteractionPoint = true; + currentTrack.restrictedPhysicalStepLength = false; } } } template -__global__ void Physics2(adept::TrackManager *gammas, Secondaries secondaries, Scoring *userScoring, - VolAuxData const *auxDataArray) +__global__ void GammaInteractions(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(); @@ -218,20 +202,24 @@ __global__ void Physics2(adept::TrackManager *gammas, Secondaries seconda auto survive = [&]() { gammas->fNextTracks->push_back(slot); }; // Temporary solution - if (!currentTrack.reachedInteractionPoint) { + 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[currentTrack.winnerProcessIndex] = -1.0; + 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 (currentTrack.winnerProcessIndex) { + 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) { @@ -298,22 +286,22 @@ __global__ void Physics2(adept::TrackManager *gammas, Secondaries seconda } 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 + 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. @@ -324,22 +312,22 @@ __global__ void Physics2(adept::TrackManager *gammas, Secondaries seconda } 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 + 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; @@ -349,9 +337,9 @@ __global__ void Physics2(adept::TrackManager *gammas, Secondaries seconda const double theLowEnergyThreshold = 1 * copcore::units::eV; const double bindingEnergy = G4HepEmGammaInteractionPhotoelectric::SelectElementBindingEnergy( - &g4HepEmData, auxData.fMCIndex, currentTrack.PEmxSec, currentTrack.eKin, &rnge); + &g4HepEmData, auxData.fMCIndex, gammaTrack.GetPEmxSec(), currentTrack.eKin, &rnge); + double edep = bindingEnergy; - double edep = bindingEnergy; const double photoElecE = currentTrack.eKin - edep; if (photoElecE > theLowEnergyThreshold) { // Create a secondary electron and sample directions. @@ -370,27 +358,29 @@ __global__ void Physics2(adept::TrackManager *gammas, Secondaries seconda } 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 + 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; } } } -} \ No newline at end of file +}