From c4f2c4f7b5df21c8a85383ba0c78854679fbd47f Mon Sep 17 00:00:00 2001 From: JuanGonzalezCaminero Date: Tue, 24 Sep 2024 16:35:11 +0200 Subject: [PATCH] Add kernels --- include/AdePT/kernels/electrons_split.cuh | 517 ++++++++++++++++++++++ include/AdePT/kernels/gammas_split.cuh | 411 +++++++++++++++++ 2 files changed, 928 insertions(+) create mode 100644 include/AdePT/kernels/electrons_split.cuh create mode 100644 include/AdePT/kernels/gammas_split.cuh diff --git a/include/AdePT/kernels/electrons_split.cuh b/include/AdePT/kernels/electrons_split.cuh new file mode 100644 index 00000000..2a210bac --- /dev/null +++ b/include/AdePT/kernels/electrons_split.cuh @@ -0,0 +1,517 @@ +// 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; + } + } + } +} \ No newline at end of file diff --git a/include/AdePT/kernels/gammas_split.cuh b/include/AdePT/kernels/gammas_split.cuh new file mode 100644 index 00000000..9b85e095 --- /dev/null +++ b/include/AdePT/kernels/gammas_split.cuh @@ -0,0 +1,411 @@ +// 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; + } + } + } +}