From c5cce6e9ee72b942e8d797c758aeef045a80425e Mon Sep 17 00:00:00 2001 From: JuanGonzalezCaminero Date: Thu, 14 Mar 2024 16:34:52 +0100 Subject: [PATCH] Remove TestEM3, add test for example1 --- README.md | 13 - examples/Example1/CMakeLists.txt | 5 + examples/Example1/macros/example1.mac.in | 4 +- test/CMakeLists.txt | 1 - test/TestEm3/CMakeLists.txt | 25 -- test/TestEm3/README.md | 53 --- test/TestEm3/TestEm3.cpp | 294 -------------- test/TestEm3/TestEm3.cu | 466 ----------------------- test/TestEm3/TestEm3.cuh | 141 ------- test/TestEm3/TestEm3.h | 38 -- test/TestEm3/electrons.cu | 369 ------------------ test/TestEm3/gammas.cu | 240 ------------ 12 files changed, 7 insertions(+), 1642 deletions(-) delete mode 100644 test/TestEm3/CMakeLists.txt delete mode 100644 test/TestEm3/README.md delete mode 100644 test/TestEm3/TestEm3.cpp delete mode 100644 test/TestEm3/TestEm3.cu delete mode 100644 test/TestEm3/TestEm3.cuh delete mode 100644 test/TestEm3/TestEm3.h delete mode 100644 test/TestEm3/electrons.cu delete mode 100644 test/TestEm3/gammas.cu diff --git a/README.md b/README.md index 998b6554..abec4836 100644 --- a/README.md +++ b/README.md @@ -127,19 +127,6 @@ target_link_libraries(example_target ${AdePT_LIBRARIES}) ``` -### Targets dependent on vecgeomcuda - -When using AdePT_LIBRARIES we are linking against the library AdePT::AdePT_cuda, which in turn links to the static library VecGeom::vecgeomcuda_static. For most purposes this will be fine, however, in case the final executable uses some feature from vecgeomcuda and needs to link against it, this aproach will not work. - -In order to address this we provide AdePT::AdePT_cuda_standalone, which links with the shared VecGeom::vecgeomcuda. This allows the final executable to correctly link with VecGeom::vecgeomcuda_static. This library is provided as part of AdePT_LIBRARIES_EXTRA. - -``` -target_link_libraries(example_target - - VecGeom::vecgeomcuda_static - ${AdePT_LIBRARIES_EXTRA}) -``` - ## Copyright AdePT code is Copyright (C) CERN, 2020, for the benefit of the AdePT project. diff --git a/examples/Example1/CMakeLists.txt b/examples/Example1/CMakeLists.txt index b3e8d282..4351772f 100644 --- a/examples/Example1/CMakeLists.txt +++ b/examples/Example1/CMakeLists.txt @@ -63,7 +63,12 @@ target_link_libraries(example1 ) # Install macros and geometry file +SET(GDML ${PROJECT_BINARY_DIR}/cms2018_sd.gdml) configure_file("macros/example1.mac.in" "${PROJECT_BINARY_DIR}/example1.mac") configure_file("macros/example1_ttbar.mac.in" "${PROJECT_BINARY_DIR}/example1_ttbar.mac") # Tests + +add_test(NAME example1 + COMMAND $ -m ${PROJECT_BINARY_DIR}/example1.mac +) \ No newline at end of file diff --git a/examples/Example1/macros/example1.mac.in b/examples/Example1/macros/example1.mac.in index e3860d35..25074c65 100644 --- a/examples/Example1/macros/example1.mac.in +++ b/examples/Example1/macros/example1.mac.in @@ -16,9 +16,9 @@ ## /adept/setSeed 1 -/detector/filename cms2018_sd.gdml +/detector/filename @GDML@ # Temporary workaround since we don't have a G4 to VecGeom converter -/adept/setVecGeomGDML cms2018_sd.gdml +/adept/setVecGeomGDML @GDML@ /adept/setVerbosity 0 ## Threshold for buffering tracks before sending to GPU /adept/setTransportBufferThreshold 2000 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 53f062a2..0da46e7d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -39,7 +39,6 @@ file(DOWNLOAD https://gitlab.cern.ch/VecGeom/VecGeom/raw/v1.2.0/persistency/gdml # Detailed tests #----------------------------------------------------------------------------# add_subdirectory(testField) -# add_subdirectory(TestEm3) #----------------------------------------------------------------------------# # Link/RDC tests diff --git a/test/TestEm3/CMakeLists.txt b/test/TestEm3/CMakeLists.txt deleted file mode 100644 index b5ebe9c6..00000000 --- a/test/TestEm3/CMakeLists.txt +++ /dev/null @@ -1,25 +0,0 @@ -# SPDX-FileCopyrightText: 2021 CERN -# SPDX-License-Identifier: Apache-2.0 - -if(NOT TARGET G4HepEm::g4HepEm) - message(STATUS "Disabling TestEm3 (needs G4HepEm)") - return() -endif() - -# TestEm3 for validation of particle transportation with GPUs. -add_executable(TestEm3 TestEm3.cu electrons.cu gammas.cu TestEm3.cpp) -target_include_directories(TestEm3 PRIVATE ${PROJECT_SOURCE_DIR}/include) -target_link_libraries(TestEm3 - PRIVATE - CopCore - AdePT_G4_integration - VecGeom::vecgeomcuda_static -) -set_target_properties(TestEm3 - PROPERTIES - CUDA_SEPARABLE_COMPILATION ON - CUDA_RESOLVE_DEVICE_SYMBOLS ON -) -set(CMAKE_EXE_LINKER_FLAGS "-Wl,--disable-new-dtags") - -add_test(NAME TestEm3 COMMAND TestEm3) diff --git a/test/TestEm3/README.md b/test/TestEm3/README.md deleted file mode 100644 index 31fd25cd..00000000 --- a/test/TestEm3/README.md +++ /dev/null @@ -1,53 +0,0 @@ - - -## TestEm3 - -Test application for the validation of particle transportation on GPUs, with the same workflow as Example 11: - - * geometry via VecGeom and optionally a magnetic field with constant Bz, - * physics processes for e-/e+ (including MSC) and gammas using G4HepEm. - -Electrons, positrons, and gammas are stored in separate containers in device memory. -Free positions in the storage are handed out with monotonic slot numbers, slots are not reused. -Active tracks are passed via three queues per particle type (see `struct ParticleQueues`). -Results are reproducible using one RANLUX++ state per track. - -Additionally, the kernels score energy deposit and the charged track length per volume. - -### Kernels - -This example uses one stream per particle type to launch kernels asynchronously. -They are synchronized via a fourth stream using CUDA events. - -#### `TransportElectrons` - -1. Obtain safety unless the track is currently on a boundary. -2. Determine physics step limit, including conversion to geometric step length according to MSC. -3. Query geometry (or optionally magnetic field) to get geometry step length. -4. Convert geometry to true step length according to MSC, apply net direction change and displacement. -5. Apply continuous effects; kill track if stopped. -6. If the particle reaches a boundary, perform relocation. -7. If not, and if there is a discrete process: - 1. Sample the final state. - 2. Update the primary and produce secondaries. - -#### `TransportGammas` - -1. Determine the physics step limit. -2. Query VecGeom to get geometry step length (no magnetic field for neutral particles!). -3. If the particle reaches a boundary, perform relocation. -4. If not, and if there is a discrete process: - 1. Sample the final state. - 2. Update the primary and produce secondaries. - -#### `FinishIteration` - -Clear the queues and return the tracks in flight. -This kernel runs after all secondary particles were produced. - -#### `InitPrimaries` and `InitParticleQueues` - -Used to initialize multiple primary particles with separate seeds. diff --git a/test/TestEm3/TestEm3.cpp b/test/TestEm3/TestEm3.cpp deleted file mode 100644 index d9afd742..00000000 --- a/test/TestEm3/TestEm3.cpp +++ /dev/null @@ -1,294 +0,0 @@ -// SPDX-FileCopyrightText: 2021 CERN -// SPDX-License-Identifier: Apache-2.0 - -#include "TestEm3.h" - -#include -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include - -#include - -#include -#include -#include -#include -#include -#include - -#ifdef ADEPT_USE_SURF -#include -#include -#endif - -#include - -const char *WorldMaterial = "G4_Galactic"; -const char *GapMaterial = "G4_Pb"; -const char *AbsorberMaterial = "G4_lAr"; - -enum MaterialCutCouples { - WorldMC = 0, - GapMC, - AbsorberMC, -}; - -constexpr G4double ProductionCut = 0.7 * mm; - -constexpr double CalorSizeYZ = 40 * copcore::units::cm; -constexpr int NbOfLayers = 50; -constexpr int NbOfAbsorbers = 2; -constexpr double GapThickness = 2.3 * copcore::units::mm; -constexpr double AbsorberThickness = 5.7 * copcore::units::mm; - -constexpr double LayerThickness = GapThickness + AbsorberThickness; -constexpr double CalorThickness = NbOfLayers * LayerThickness; - -constexpr double WorldSizeX = 1.2 * CalorThickness; -constexpr double WorldSizeYZ = 1.2 * CalorSizeYZ; - -static void InitGeant4() -{ - // --- Create materials. - G4Material *world = G4NistManager::Instance()->FindOrBuildMaterial(WorldMaterial); - G4Material *gap = G4NistManager::Instance()->FindOrBuildMaterial(GapMaterial); - G4Material *absorber = G4NistManager::Instance()->FindOrBuildMaterial(AbsorberMaterial); - // - // --- Define a world. - G4double worldDim = 1 * m; - G4Box *worldBox = new G4Box("world", worldDim, worldDim, worldDim); - G4LogicalVolume *worldLog = new G4LogicalVolume(worldBox, world, "world"); - G4PVPlacement *worldPlaced = new G4PVPlacement(nullptr, {}, worldLog, "world", nullptr, false, 0); - // --- Define two boxes to use the materials. - G4double boxDim = 0.25 * m; - G4double gapBoxPos = 0.5 * boxDim; - G4double absorberBoxPos = 1.5 * boxDim; - G4Box *box = new G4Box("", boxDim, boxDim, boxDim); - G4LogicalVolume *gapLog = new G4LogicalVolume(box, gap, "gap"); - G4LogicalVolume *absorberLog = new G4LogicalVolume(box, absorber, "absorber"); - new G4PVPlacement(nullptr, {gapBoxPos, gapBoxPos, gapBoxPos}, gapLog, "gap", worldLog, false, 0); - new G4PVPlacement(nullptr, {absorberBoxPos, absorberBoxPos, absorberBoxPos}, absorberLog, "absorber", worldLog, false, - 0); - // - // --- Create particles that have secondary production threshold. - G4Gamma::Gamma(); - G4Electron::Electron(); - G4Positron::Positron(); - G4Proton::Proton(); - G4ParticleTable *partTable = G4ParticleTable::GetParticleTable(); - partTable->SetReadiness(); - // - // --- Create production - cuts object and set the secondary production threshold. - G4ProductionCuts *productionCuts = new G4ProductionCuts(); - productionCuts->SetProductionCut(ProductionCut); - // - // --- Register a region for the world. - G4Region *reg = new G4Region("default"); - reg->AddRootLogicalVolume(worldLog); - reg->UsedInMassGeometry(true); - reg->SetProductionCuts(productionCuts); - // - // --- Update the couple tables. - G4ProductionCutsTable *theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable(); - theCoupleTable->UpdateCoupleTable(worldPlaced); - // - // --- Set the world volume to fix initialization of G4SafetyHelper (used by G4UrbanMscModel) - G4TransportationManager::GetTransportationManager()->SetWorldForTracking(worldPlaced); - // - // --- Set MSC range factor to match G4HepEm physics lists. - G4EmParameters *param = G4EmParameters::Instance(); - param->SetDefaults(); - param->SetMscRangeFactor(0.04); -} - -const void CreateVecGeomWorld() -{ - auto worldSolid = new vecgeom::UnplacedBox(0.5 * WorldSizeX, 0.5 * WorldSizeYZ, 0.5 * WorldSizeYZ); - auto worldLogic = new vecgeom::LogicalVolume("World", worldSolid); - vecgeom::VPlacedVolume *worldPlaced = worldLogic->Place(); - - // - // Calorimeter - // - auto calorSolid = new vecgeom::UnplacedBox(0.5 * CalorThickness, 0.5 * CalorSizeYZ, 0.5 * CalorSizeYZ); - auto calorLogic = new vecgeom::LogicalVolume("Calorimeter", calorSolid); - vecgeom::Transformation3D origin; - worldLogic->PlaceDaughter(calorLogic, &origin); - - // - // Layers - // - auto layerSolid = new vecgeom::UnplacedBox(0.5 * LayerThickness, 0.5 * CalorSizeYZ, 0.5 * CalorSizeYZ); - - // - // Absorbers - // - auto gapSolid = new vecgeom::UnplacedBox(0.5 * GapThickness, 0.5 * CalorSizeYZ, 0.5 * CalorSizeYZ); - auto gapLogic = new vecgeom::LogicalVolume("Gap", gapSolid); - vecgeom::Transformation3D gapPlacement(-0.5 * LayerThickness + 0.5 * GapThickness, 0, 0); - - auto absorberSolid = new vecgeom::UnplacedBox(0.5 * AbsorberThickness, 0.5 * CalorSizeYZ, 0.5 * CalorSizeYZ); - auto absorberLogic = new vecgeom::LogicalVolume("Absorber", absorberSolid); - vecgeom::Transformation3D absorberPlacement(0.5 * LayerThickness - 0.5 * AbsorberThickness, 0, 0); - - // Create a new LogicalVolume per layer, we need unique IDs for scoring. - double xCenter = -0.5 * CalorThickness + 0.5 * LayerThickness; - for (int i = 0; i < NbOfLayers; i++) { - auto layerLogic = new vecgeom::LogicalVolume("Layer", layerSolid); - vecgeom::Transformation3D placement(xCenter, 0, 0); - calorLogic->PlaceDaughter(layerLogic, &placement); - - layerLogic->PlaceDaughter(gapLogic, &gapPlacement); - layerLogic->PlaceDaughter(absorberLogic, &absorberPlacement); - - xCenter += LayerThickness; - } - - vecgeom::GeoManager::Instance().SetWorldAndClose(worldPlaced); -} - -void PrintDaughters(const vecgeom::VPlacedVolume *placed, int level = 0) -{ - for (auto *daughter : placed->GetDaughters()) { - std::cout << std::setw(level * 2) << ""; - auto *trans = daughter->GetTransformation(); - std::cout << " ID " << daughter->id() << " @ " << trans->Translation() << std::endl; - PrintDaughters(daughter, level + 1); - } -} - -void PrintScoringPerVolume(const vecgeom::VPlacedVolume *placed, const ScoringPerVolume *scoring, int level = 0) -{ - for (auto *daughter : placed->GetDaughters()) { - std::cout << std::setw(level * 2) << ""; - auto id = daughter->id(); - std::cout << " ID " << id << " Charged-TrakL " << scoring->chargedTrackLength[id] / copcore::units::mm - << " mm; Energy-Dep " << scoring->energyDeposit[id] / copcore::units::MeV << " MeV" << std::endl; - PrintScoringPerVolume(daughter, scoring, level + 1); - } -} - -void InitBVH() -{ - vecgeom::cxx::BVHManager::Init(); - vecgeom::cxx::BVHManager::DeviceInit(); -} - -int main(int argc, char *argv[]) -{ -#ifndef VECGEOM_USE_NAVINDEX - std::cout << "### VecGeom must be compiled with USE_NAVINDEX support to run this.\n"; - return 1; -#endif - - OPTION_INT(cache_depth, 0); // 0 = full depth - OPTION_INT(particles, 1000); - OPTION_DOUBLE(energy, 10); // entered in GeV - energy *= copcore::units::GeV; - OPTION_INT(batch, -1); - - InitGeant4(); - CreateVecGeomWorld(); - - const vecgeom::VPlacedVolume *world = vecgeom::GeoManager::Instance().GetWorld(); -#ifdef VERBOSE - std::cout << "World (ID " << world->id() << ")" << std::endl; - PrintDaughters(world); -#endif - -#ifdef ADEPT_USE_SURF - if (!vgbrep::BrepHelper::Instance().Convert()) return 1; - vgbrep::BrepHelper::Instance().PrintSurfData(); -#endif - - // Map VecGeom volume IDs to Geant4 material-cuts couples. - constexpr int NumVolumes = 1 + 1 + NbOfLayers * (1 + NbOfAbsorbers); - int MCIndex[NumVolumes]; - // Fill world and calorimeter. - MCIndex[0] = MCIndex[1] = WorldMC; - for (int i = 2; i < NumVolumes; i += (1 + NbOfAbsorbers)) { - MCIndex[i] = WorldMC; - MCIndex[i + 1] = GapMC; - MCIndex[i + 2] = AbsorberMC; - } - // Place particles between the world boundary and the calorimeter. - double startX = -0.25 * (WorldSizeX + CalorThickness); - double chargedTrackLength[NumVolumes]; - double energyDeposit[NumVolumes]; - ScoringPerVolume scoringPerVolume; - scoringPerVolume.chargedTrackLength = chargedTrackLength; - scoringPerVolume.energyDeposit = energyDeposit; - GlobalScoring globalScoring; - - TestEm3(world, particles, energy, batch, startX, MCIndex, &scoringPerVolume, NumVolumes, &globalScoring); - - std::cout << std::endl; - std::cout << std::endl; - double meanEnergyDeposit = globalScoring.energyDeposit / particles; - std::cout << "Mean energy deposit " << (meanEnergyDeposit / copcore::units::GeV) << " GeV\n" - << "Mean number of gamma " << ((double)globalScoring.numGammas / particles) << "\n" - << "Mean number of e- " << ((double)globalScoring.numElectrons / particles) << "\n" - << "Mean number of e+ " << ((double)globalScoring.numPositrons / particles) << "\n" - << "Mean number of charged steps " << ((double)globalScoring.chargedSteps / particles) << "\n" - << "Mean number of neutral steps " << ((double)globalScoring.neutralSteps / particles) << "\n" - << "Mean number of hits " << ((double)globalScoring.hits / particles) << "\n"; - if (globalScoring.numKilled > 0) { - std::cout << "Total killed particles " << globalScoring.numKilled << "\n"; - } - std::cout << std::endl; - - // Average charged track length and energy deposit per particle. - for (int i = 0; i < NumVolumes; i++) { - chargedTrackLength[i] /= particles; - energyDeposit[i] /= particles; - } - - std::cout << std::scientific; -#ifdef VERBOSE - std::cout << std::endl; - const int id = world->id(); - std::cout << "World (ID " << world->id() << ") Charged-TrakL " << chargedTrackLength[id] / copcore::units::mm - << " Energy-Dep " << energyDeposit[id] / copcore::units::MeV << " MeV" << std::endl; - PrintScoringPerVolume(world, &scoringPerVolume); -#endif - - // Accumulate per material. - double energyDepGap = 0, energyDepAbsorber = 0; - for (int i = 2; i < NumVolumes; i += (1 + NbOfAbsorbers)) { - energyDepGap += energyDeposit[i + 1]; - energyDepAbsorber += energyDeposit[i + 2]; - } - - std::cout << std::endl; - std::cout << " " << GapMaterial << ": " << energyDepGap / copcore::units::GeV << " GeV" << std::endl; - std::cout << " " << AbsorberMaterial << ": " << energyDepAbsorber / copcore::units::GeV << " GeV" << std::endl; - // Accumulate per layer. - std::cout << std::endl; - std::cout << " Layer Charged-TrakL [mm] Energy-Dep [MeV]" << std::endl; - for (int i = 0; i < NbOfLayers; i++) { - int layerVolume = 2 + i * (1 + NbOfAbsorbers); - double chargedTrackLen = chargedTrackLength[layerVolume + 1] + chargedTrackLength[layerVolume + 2]; - double energyDep = energyDeposit[layerVolume + 1] + energyDeposit[layerVolume + 2]; - std::cout << std::setw(5) << i << std::setw(20) << chargedTrackLen / copcore::units::mm << std::setw(20) - << energyDep / copcore::units::MeV << std::endl; - } -} diff --git a/test/TestEm3/TestEm3.cu b/test/TestEm3/TestEm3.cu deleted file mode 100644 index 4ba5d553..00000000 --- a/test/TestEm3/TestEm3.cu +++ /dev/null @@ -1,466 +0,0 @@ -// SPDX-FileCopyrightText: 2021 CERN -// SPDX-License-Identifier: Apache-2.0 - -#include "TestEm3.h" -#include "TestEm3.cuh" - -#include -#include -#include - -#include -#include -#include - -#include -#include -#ifdef VECGEOM_ENABLE_CUDA -#include -#endif -#ifdef ADEPT_USE_SURF -#include -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -__constant__ __device__ struct G4HepEmParameters g4HepEmPars; -__constant__ __device__ struct G4HepEmData g4HepEmData; - -__constant__ __device__ int *MCIndex = nullptr; - -struct G4HepEmState { - G4HepEmData data; - G4HepEmParameters parameters; -}; - -static G4HepEmState *InitG4HepEm() -{ - G4HepEmState *state = new G4HepEmState; - InitG4HepEmData(&state->data); - InitHepEmParameters(&state->parameters); - - InitMaterialAndCoupleData(&state->data, &state->parameters); - - InitElectronData(&state->data, &state->parameters, true); - InitElectronData(&state->data, &state->parameters, false); - InitGammaData(&state->data, &state->parameters); - - G4HepEmMatCutData *cutData = state->data.fTheMatCutData; - std::cout << "fNumG4MatCuts = " << cutData->fNumG4MatCuts << ", fNumMatCutData = " << cutData->fNumMatCutData - << std::endl; - - // Copy to GPU. - CopyG4HepEmDataToGPU(&state->data); - COPCORE_CUDA_CHECK(cudaMemcpyToSymbol(g4HepEmPars, &state->parameters, sizeof(G4HepEmParameters))); - - // Create G4HepEmData with the device pointers. - G4HepEmData dataOnDevice; - dataOnDevice.fTheMatCutData = state->data.fTheMatCutData_gpu; - dataOnDevice.fTheMaterialData = state->data.fTheMaterialData_gpu; - dataOnDevice.fTheElementData = state->data.fTheElementData_gpu; - dataOnDevice.fTheElectronData = state->data.fTheElectronData_gpu; - dataOnDevice.fThePositronData = state->data.fThePositronData_gpu; - dataOnDevice.fTheSBTableData = state->data.fTheSBTableData_gpu; - dataOnDevice.fTheGammaData = state->data.fTheGammaData_gpu; - // The other pointers should never be used. - dataOnDevice.fTheMatCutData_gpu = nullptr; - dataOnDevice.fTheMaterialData_gpu = nullptr; - dataOnDevice.fTheElementData_gpu = nullptr; - dataOnDevice.fTheElectronData_gpu = nullptr; - dataOnDevice.fThePositronData_gpu = nullptr; - dataOnDevice.fTheSBTableData_gpu = nullptr; - dataOnDevice.fTheGammaData_gpu = nullptr; - - COPCORE_CUDA_CHECK(cudaMemcpyToSymbol(g4HepEmData, &dataOnDevice, sizeof(G4HepEmData))); - - return state; -} - -static void FreeG4HepEm(G4HepEmState *state) -{ - FreeG4HepEmData(&state->data); - delete state; -} - -// A bundle of queues per particle type: -// * Two for active particles, one for the current iteration and the second for the next. -struct ParticleQueues { - adept::MParray *currentlyActive; - adept::MParray *nextActive; - - void SwapActive() { std::swap(currentlyActive, nextActive); } -}; - -struct ParticleType { - Track *tracks; - SlotManager *slotManager; - ParticleQueues queues; - cudaStream_t stream; - cudaEvent_t event; - - enum { - Electron = 0, - Positron = 1, - Gamma = 2, - - NumParticleTypes, - }; -}; - -// A bundle of queues for the three particle types. -struct AllParticleQueues { - ParticleQueues queues[ParticleType::NumParticleTypes]; -}; - -// Kernel to initialize the set of queues per particle type. -__global__ void InitParticleQueues(ParticleQueues queues, size_t Capacity) -{ - adept::MParray::MakeInstanceAt(Capacity, queues.currentlyActive); - adept::MParray::MakeInstanceAt(Capacity, queues.nextActive); -} - -// Kernel function to initialize a set of primary particles. -__global__ void InitPrimaries(ParticleGenerator generator, int startEvent, int numEvents, double energy, double startX, - const vecgeom::VPlacedVolume *world, GlobalScoring *globalScoring) -{ - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < numEvents; i += blockDim.x * gridDim.x) { - Track &track = generator.NextTrack(); - - track.rngState.SetSeed(startEvent + i); - track.energy = energy; - track.numIALeft[0] = -1.0; - track.numIALeft[1] = -1.0; - track.numIALeft[2] = -1.0; - - track.initialRange = -1.0; - track.dynamicRangeFactor = -1.0; - track.tlimitMin = -1.0; - - track.pos = {vecgeom::Precision(startX), 0, 0}; - track.dir = {1.0, 0, 0}; - track.navState.Clear(); - BVHNavigator::LocatePointIn(world, track.pos, track.navState, true); - - atomicAdd(&globalScoring->numElectrons, 1); - } -} - -// A data structure to transfer statistics after each iteration. -struct Stats { - int inFlight[ParticleType::NumParticleTypes]; -}; - -// Finish iteration: clear queues and fill statistics. -__global__ void FinishIteration(AllParticleQueues all, Stats *stats) -{ - for (int i = 0; i < ParticleType::NumParticleTypes; i++) { - all.queues[i].currentlyActive->clear(); - stats->inFlight[i] = all.queues[i].nextActive->size(); - } -} - -__global__ void ClearQueue(adept::MParray *queue) -{ - queue->clear(); -} - -void TestEm3(const vecgeom::cxx::VPlacedVolume *world, int numParticles, double energy, int batch, double startX, - const int *MCIndex_host, ScoringPerVolume *scoringPerVolume_host, int numVolumes, - GlobalScoring *globalScoring_host) -{ - auto &cudaManager = vecgeom::cxx::CudaManager::Instance(); - cudaManager.LoadGeometry(world); - cudaManager.Synchronize(); - - const vecgeom::cuda::VPlacedVolume *world_dev = cudaManager.world_gpu(); - - InitBVH(); - -#ifdef ADEPT_USE_SURF - auto &surfdata = vgbrep::SurfData::Instance(); - vgbrep::BrepCudaManager::Instance().TransferSurfData(surfdata); -#endif - - G4HepEmState *state = InitG4HepEm(); - - // Transfer MC indices. - int *MCIndex_dev = nullptr; - COPCORE_CUDA_CHECK(cudaMalloc(&MCIndex_dev, sizeof(int) * numVolumes)); - COPCORE_CUDA_CHECK(cudaMemcpy(MCIndex_dev, MCIndex_host, sizeof(int) * numVolumes, cudaMemcpyHostToDevice)); - COPCORE_CUDA_CHECK(cudaMemcpyToSymbol(MCIndex, &MCIndex_dev, sizeof(int *))); - - // Capacity of the different containers aka the maximum number of particles. - constexpr int Capacity = 1024 * 1024; - - std::cout << "INFO: capacity of containers set to " << Capacity << std::endl; - if (batch == -1) { - // Rule of thumb: at most 1000 particles of one type per GeV primary. - batch = Capacity / ((int)energy / copcore::units::GeV) / 1000; - } else if (batch < 1) { - batch = 1; - } - std::cout << "INFO: batching " << batch << " particles for transport on the GPU" << std::endl; - if (BzFieldValue != 0) { - std::cout << "INFO: running with field Bz = " << BzFieldValue / copcore::units::tesla << " T" << std::endl; - } else { - std::cout << "INFO: running with magnetic field OFF" << std::endl; - } - - // Allocate structures to manage tracks of an implicit type: - // * memory to hold the actual Track elements, - // * objects to manage slots inside the memory, - // * queues of slots to remember active particle and those needing relocation, - // * a stream and an event for synchronization of kernels. - constexpr size_t TracksSize = sizeof(Track) * Capacity; - constexpr size_t ManagerSize = sizeof(SlotManager); - const size_t QueueSize = adept::MParray::SizeOfInstance(Capacity); - - ParticleType particles[ParticleType::NumParticleTypes]; - for (int i = 0; i < ParticleType::NumParticleTypes; i++) { - COPCORE_CUDA_CHECK(cudaMalloc(&particles[i].tracks, TracksSize)); - - COPCORE_CUDA_CHECK(cudaMalloc(&particles[i].slotManager, ManagerSize)); - - COPCORE_CUDA_CHECK(cudaMalloc(&particles[i].queues.currentlyActive, QueueSize)); - COPCORE_CUDA_CHECK(cudaMalloc(&particles[i].queues.nextActive, QueueSize)); - InitParticleQueues<<<1, 1>>>(particles[i].queues, Capacity); - - COPCORE_CUDA_CHECK(cudaStreamCreate(&particles[i].stream)); - COPCORE_CUDA_CHECK(cudaEventCreate(&particles[i].event)); - } - COPCORE_CUDA_CHECK(cudaDeviceSynchronize()); - - ParticleType &electrons = particles[ParticleType::Electron]; - ParticleType &positrons = particles[ParticleType::Positron]; - ParticleType &gammas = particles[ParticleType::Gamma]; - - // Create a stream to synchronize kernels of all particle types. - cudaStream_t stream; - COPCORE_CUDA_CHECK(cudaStreamCreate(&stream)); - - // Allocate memory to score charged track length and energy deposit per volume. - double *chargedTrackLength = nullptr; - COPCORE_CUDA_CHECK(cudaMalloc(&chargedTrackLength, sizeof(double) * numVolumes)); - COPCORE_CUDA_CHECK(cudaMemset(chargedTrackLength, 0, sizeof(double) * numVolumes)); - double *energyDeposit = nullptr; - COPCORE_CUDA_CHECK(cudaMalloc(&energyDeposit, sizeof(double) * numVolumes)); - COPCORE_CUDA_CHECK(cudaMemset(energyDeposit, 0, sizeof(double) * numVolumes)); - - // Allocate and initialize scoring and statistics. - GlobalScoring *globalScoring = nullptr; - COPCORE_CUDA_CHECK(cudaMalloc(&globalScoring, sizeof(GlobalScoring))); - COPCORE_CUDA_CHECK(cudaMemset(globalScoring, 0, sizeof(GlobalScoring))); - - ScoringPerVolume *scoringPerVolume = nullptr; - ScoringPerVolume scoringPerVolume_devPtrs; - scoringPerVolume_devPtrs.chargedTrackLength = chargedTrackLength; - scoringPerVolume_devPtrs.energyDeposit = energyDeposit; - COPCORE_CUDA_CHECK(cudaMalloc(&scoringPerVolume, sizeof(ScoringPerVolume))); - COPCORE_CUDA_CHECK( - cudaMemcpy(scoringPerVolume, &scoringPerVolume_devPtrs, sizeof(ScoringPerVolume), cudaMemcpyHostToDevice)); - - Stats *stats_dev = nullptr; - COPCORE_CUDA_CHECK(cudaMalloc(&stats_dev, sizeof(Stats))); - Stats *stats = nullptr; - COPCORE_CUDA_CHECK(cudaMallocHost(&stats, sizeof(Stats))); - - // Allocate memory to hold a "vanilla" SlotManager to initialize for each batch. - SlotManager slotManagerInit(Capacity); - SlotManager *slotManagerInit_dev = nullptr; - COPCORE_CUDA_CHECK(cudaMalloc(&slotManagerInit_dev, sizeof(SlotManager))); - COPCORE_CUDA_CHECK(cudaMemcpy(slotManagerInit_dev, &slotManagerInit, sizeof(SlotManager), cudaMemcpyHostToDevice)); - - vecgeom::Stopwatch timer; - timer.Start(); - - std::cout << std::endl << "Simulating particles "; - const bool detailed = (numParticles / batch) < 50; - if (!detailed) { - std::cout << "... " << std::flush; - } - - unsigned long long killed = 0; - - for (int startEvent = 1; startEvent <= numParticles; startEvent += batch) { - if (detailed) { - std::cout << startEvent << " ... " << std::flush; - } - int left = numParticles - startEvent + 1; - int chunk = std::min(left, batch); - - for (int i = 0; i < ParticleType::NumParticleTypes; i++) { - COPCORE_CUDA_CHECK(cudaMemcpyAsync(particles[i].slotManager, slotManagerInit_dev, ManagerSize, - cudaMemcpyDeviceToDevice, stream)); - } - - // Initialize primary particles. - constexpr int InitThreads = 32; - int initBlocks = (chunk + InitThreads - 1) / InitThreads; - ParticleGenerator electronGenerator(electrons.tracks, electrons.slotManager, electrons.queues.currentlyActive); - InitPrimaries<<>>(electronGenerator, startEvent, chunk, energy, startX, - world_dev, globalScoring); - COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); - - stats->inFlight[ParticleType::Electron] = chunk; - stats->inFlight[ParticleType::Positron] = 0; - stats->inFlight[ParticleType::Gamma] = 0; - - constexpr int MaxBlocks = 1024; - constexpr int TransportThreads = 32; - int transportBlocks; - - int inFlight; - int loopingNo = 0; - int previousElectrons = -1, previousPositrons = -1; - - do { - Secondaries secondaries = { - .electrons = {electrons.tracks, electrons.slotManager, electrons.queues.nextActive}, - .positrons = {positrons.tracks, positrons.slotManager, positrons.queues.nextActive}, - .gammas = {gammas.tracks, gammas.slotManager, gammas.queues.nextActive}, - }; - - // *** ELECTRONS *** - int numElectrons = stats->inFlight[ParticleType::Electron]; - if (numElectrons > 0) { - transportBlocks = (numElectrons + TransportThreads - 1) / TransportThreads; - transportBlocks = std::min(transportBlocks, MaxBlocks); - - TransportElectrons<<>>( - electrons.tracks, electrons.queues.currentlyActive, secondaries, electrons.queues.nextActive, globalScoring, - scoringPerVolume); - - COPCORE_CUDA_CHECK(cudaEventRecord(electrons.event, electrons.stream)); - COPCORE_CUDA_CHECK(cudaStreamWaitEvent(stream, electrons.event, 0)); - } - - // *** POSITRONS *** - int numPositrons = stats->inFlight[ParticleType::Positron]; - if (numPositrons > 0) { - transportBlocks = (numPositrons + TransportThreads - 1) / TransportThreads; - transportBlocks = std::min(transportBlocks, MaxBlocks); - - TransportPositrons<<>>( - positrons.tracks, positrons.queues.currentlyActive, secondaries, positrons.queues.nextActive, globalScoring, - scoringPerVolume); - - COPCORE_CUDA_CHECK(cudaEventRecord(positrons.event, positrons.stream)); - COPCORE_CUDA_CHECK(cudaStreamWaitEvent(stream, positrons.event, 0)); - } - - // *** GAMMAS *** - int numGammas = stats->inFlight[ParticleType::Gamma]; - if (numGammas > 0) { - transportBlocks = (numGammas + TransportThreads - 1) / TransportThreads; - transportBlocks = std::min(transportBlocks, MaxBlocks); - - TransportGammas<<>>( - gammas.tracks, gammas.queues.currentlyActive, secondaries, gammas.queues.nextActive, globalScoring, - scoringPerVolume); - - COPCORE_CUDA_CHECK(cudaEventRecord(gammas.event, gammas.stream)); - COPCORE_CUDA_CHECK(cudaStreamWaitEvent(stream, gammas.event, 0)); - } - - // *** END OF TRANSPORT *** - - // The events ensure synchronization before finishing this iteration and - // copying the Stats back to the host. - AllParticleQueues queues = {{electrons.queues, positrons.queues, gammas.queues}}; - FinishIteration<<<1, 1, 0, stream>>>(queues, stats_dev); - COPCORE_CUDA_CHECK(cudaMemcpyAsync(stats, stats_dev, sizeof(Stats), cudaMemcpyDeviceToHost, stream)); - - // Finally synchronize all kernels. - COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); - - // Count the number of particles in flight. - inFlight = 0; - for (int i = 0; i < ParticleType::NumParticleTypes; i++) { - inFlight += stats->inFlight[i]; - } - - // Swap the queues for the next iteration. - electrons.queues.SwapActive(); - positrons.queues.SwapActive(); - gammas.queues.SwapActive(); - - // Check if only charged particles are left that are looping. - numElectrons = stats->inFlight[ParticleType::Electron]; - numPositrons = stats->inFlight[ParticleType::Positron]; - numGammas = stats->inFlight[ParticleType::Gamma]; - if (numElectrons == previousElectrons && numPositrons == previousPositrons && numGammas == 0) { - loopingNo++; - } else { - previousElectrons = numElectrons; - previousPositrons = numPositrons; - loopingNo = 0; - } - - } while (inFlight > 0 && loopingNo < 200); - - if (inFlight > 0) { - killed += inFlight; - for (int i = 0; i < ParticleType::NumParticleTypes; i++) { - ParticleType &pType = particles[i]; - int inFlightParticles = stats->inFlight[i]; - if (inFlightParticles == 0) { - continue; - } - - ClearQueue<<<1, 1, 0, stream>>>(pType.queues.currentlyActive); - } - COPCORE_CUDA_CHECK(cudaStreamSynchronize(stream)); - } - } - std::cout << "done!" << std::endl; - - auto time = timer.Stop(); - std::cout << "Run time: " << time << "\n"; - - // Transfer back scoring. - COPCORE_CUDA_CHECK(cudaMemcpy(globalScoring_host, globalScoring, sizeof(GlobalScoring), cudaMemcpyDeviceToHost)); - globalScoring_host->numKilled = killed; - - // Transfer back the scoring per volume (charged track length and energy deposit). - COPCORE_CUDA_CHECK(cudaMemcpy(scoringPerVolume_host->chargedTrackLength, scoringPerVolume_devPtrs.chargedTrackLength, - sizeof(double) * numVolumes, cudaMemcpyDeviceToHost)); - COPCORE_CUDA_CHECK(cudaMemcpy(scoringPerVolume_host->energyDeposit, scoringPerVolume_devPtrs.energyDeposit, - sizeof(double) * numVolumes, cudaMemcpyDeviceToHost)); - - // Free resources. - COPCORE_CUDA_CHECK(cudaFree(MCIndex_dev)); - COPCORE_CUDA_CHECK(cudaFree(chargedTrackLength)); - COPCORE_CUDA_CHECK(cudaFree(energyDeposit)); - - COPCORE_CUDA_CHECK(cudaFree(globalScoring)); - COPCORE_CUDA_CHECK(cudaFree(scoringPerVolume)); - COPCORE_CUDA_CHECK(cudaFree(stats_dev)); - COPCORE_CUDA_CHECK(cudaFreeHost(stats)); - COPCORE_CUDA_CHECK(cudaFree(slotManagerInit_dev)); - - COPCORE_CUDA_CHECK(cudaStreamDestroy(stream)); - - for (int i = 0; i < ParticleType::NumParticleTypes; i++) { - COPCORE_CUDA_CHECK(cudaFree(particles[i].tracks)); - COPCORE_CUDA_CHECK(cudaFree(particles[i].slotManager)); - - COPCORE_CUDA_CHECK(cudaFree(particles[i].queues.currentlyActive)); - COPCORE_CUDA_CHECK(cudaFree(particles[i].queues.nextActive)); - - COPCORE_CUDA_CHECK(cudaStreamDestroy(particles[i].stream)); - COPCORE_CUDA_CHECK(cudaEventDestroy(particles[i].event)); - } - - FreeG4HepEm(state); -} diff --git a/test/TestEm3/TestEm3.cuh b/test/TestEm3/TestEm3.cuh deleted file mode 100644 index a2a33168..00000000 --- a/test/TestEm3/TestEm3.cuh +++ /dev/null @@ -1,141 +0,0 @@ -// SPDX-FileCopyrightText: 2022 CERN -// SPDX-License-Identifier: Apache-2.0 - -#ifndef TESTEM3_CUH -#define TESTEM3_CUH - -#include "TestEm3.h" - -#include -#include -#include - -#include -#include -#include - -#include -#include - -// A data structure to represent a particle track. The particle type is implicit -// by the queue and not stored in memory. -struct Track { - using Precision = vecgeom::Precision; - RanluxppDouble rngState; - double energy; - double numIALeft[3]; - double initialRange; - double dynamicRangeFactor; - double tlimitMin; - - vecgeom::Vector3D pos; - vecgeom::Vector3D dir; - vecgeom::NavStateIndex navState; - - __host__ __device__ double Uniform() { return rngState.Rndm(); } - - __host__ __device__ void InitAsSecondary(const vecgeom::Vector3D &parentPos, - const vecgeom::NavStateIndex &parentNavState) - { - // The caller is responsible to branch a new RNG state and to set the energy. - this->numIALeft[0] = -1.0; - this->numIALeft[1] = -1.0; - this->numIALeft[2] = -1.0; - - this->initialRange = -1.0; - this->dynamicRangeFactor = -1.0; - this->tlimitMin = -1.0; - - // A secondary inherits the position of its parent; the caller is responsible - // to update the directions. - this->pos = parentPos; - this->navState = parentNavState; - } -}; - -#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 -// communicate the intent.) -inline __device__ double G4HepEmRandomEngine::flat() -{ - return ((RanluxppDouble *)fObject)->Rndm(); -} - -inline __device__ void G4HepEmRandomEngine::flatArray(const int size, double *vect) -{ - for (int i = 0; i < size; i++) { - vect[i] = ((RanluxppDouble *)fObject)->Rndm(); - } -} -#endif - -// A data structure to manage slots in the track storage. -class SlotManager { - adept::Atomic_t fNextSlot; - const int fMaxSlot; - -public: - __host__ __device__ SlotManager(int maxSlot) : fMaxSlot(maxSlot) { fNextSlot = 0; } - - __host__ __device__ int NextSlot() - { - int next = fNextSlot.fetch_add(1); - if (next >= fMaxSlot) return -1; - return next; - } -}; - -// A bundle of pointers to generate particles of an implicit type. -class ParticleGenerator { - Track *fTracks; - SlotManager *fSlotManager; - adept::MParray *fActiveQueue; - -public: - __host__ __device__ ParticleGenerator(Track *tracks, SlotManager *slotManager, adept::MParray *activeQueue) - : fTracks(tracks), fSlotManager(slotManager), fActiveQueue(activeQueue) - { - } - - __host__ __device__ Track &NextTrack() - { - int slot = fSlotManager->NextSlot(); - if (slot == -1) { - COPCORE_EXCEPTION("No slot available in ParticleGenerator::NextTrack"); - } - fActiveQueue->push_back(slot); - return fTracks[slot]; - } -}; - -// A bundle of generators for the three particle types. -struct Secondaries { - ParticleGenerator electrons; - ParticleGenerator positrons; - ParticleGenerator gammas; -}; - -// Kernels in different TUs. -__global__ void TransportElectrons(Track *electrons, const adept::MParray *active, Secondaries secondaries, - adept::MParray *activeQueue, GlobalScoring *globalScoring, - ScoringPerVolume *scoringPerVolume); -__global__ void TransportPositrons(Track *positrons, const adept::MParray *active, Secondaries secondaries, - adept::MParray *activeQueue, GlobalScoring *globalScoring, - ScoringPerVolume *scoringPerVolume); - -__global__ void TransportGammas(Track *gammas, const adept::MParray *active, Secondaries secondaries, - adept::MParray *activeQueue, GlobalScoring *globalScoring, - ScoringPerVolume *scoringPerVolume); - -// Constant data structures from G4HepEm accessed by the kernels. -// (defined in TestEm3.cu) -extern __constant__ __device__ struct G4HepEmParameters g4HepEmPars; -extern __constant__ __device__ struct G4HepEmData g4HepEmData; - -extern __constant__ __device__ int *MCIndex; - -// constexpr vecgeom::Precision BzFieldValue = 1 * copcore::units::tesla; -constexpr vecgeom::Precision BzFieldValue = 0; - -#endif diff --git a/test/TestEm3/TestEm3.h b/test/TestEm3/TestEm3.h deleted file mode 100644 index 452542e1..00000000 --- a/test/TestEm3/TestEm3.h +++ /dev/null @@ -1,38 +0,0 @@ -// SPDX-FileCopyrightText: 2021 CERN -// SPDX-License-Identifier: Apache-2.0 - -#ifndef TESTEM3_H -#define TESTEM3_H - -#include -#ifdef VECGEOM_ENABLE_CUDA -#include // forward declares vecgeom::cxx::VPlacedVolume -#endif - -// Data structures for scoring. The accessors must make sure to use atomic operations if needed. -struct GlobalScoring { - double energyDeposit; - // Not int to avoid overflows for more than 100,000 events; unsigned long long - // is the only other data type available for atomicAdd(). - unsigned long long chargedSteps; - unsigned long long neutralSteps; - unsigned long long hits; - unsigned long long numGammas; - unsigned long long numElectrons; - unsigned long long numPositrons; - // Not used on the device, filled in by the host. - unsigned long long numKilled; -}; - -struct ScoringPerVolume { - double *energyDeposit; - double *chargedTrackLength; -}; - -// Interface between C++ and CUDA. -void TestEm3(const vecgeom::cxx::VPlacedVolume *world, int numParticles, double energy, int batch, double startX, - const int *MCIndex, ScoringPerVolume *scoringPerVolume, int numVolumes, GlobalScoring *globalScoring); - -void InitBVH(); - -#endif diff --git a/test/TestEm3/electrons.cu b/test/TestEm3/electrons.cu deleted file mode 100644 index dee8b825..00000000 --- a/test/TestEm3/electrons.cu +++ /dev/null @@ -1,369 +0,0 @@ -// SPDX-FileCopyrightText: 2022 CERN -// SPDX-License-Identifier: Apache-2.0 - -#include "TestEm3.cuh" - -#include -#include -#ifdef ADEPT_USE_SURF -#include -#endif -#include - -#include - -#define NOFLUCTUATION - -#include -#include -#include -#include -#include -#include -// Pull in implementation. -#include -#include -#include -#include -#include -#include -#include - -// Compute the physics and geometry step limit, transport the electrons while -// applying the continuous effects and maybe a discrete process that could -// generate secondaries. -template -static __device__ __forceinline__ void TransportElectrons(Track *electrons, const adept::MParray *active, - Secondaries &secondaries, adept::MParray *activeQueue, - GlobalScoring *globalScoring, - ScoringPerVolume *scoringPerVolume) -{ -#ifdef ADEPT_USE_SURF - using Navigator = SurfNavigator; -#else - using Navigator = BVHNavigator; -#endif - constexpr int Charge = IsElectron ? -1 : 1; - constexpr double Mass = copcore::units::kElectronMassC2; - fieldPropagatorConstBz fieldPropagatorBz(BzFieldValue); - - int activeSize = active->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*active)[i]; - Track ¤tTrack = electrons[slot]; - auto energy = currentTrack.energy; - auto pos = currentTrack.pos; - auto dir = currentTrack.dir; - auto navState = currentTrack.navState; - const int volumeID = navState.Top()->id(); - const int theMCIndex = MCIndex[volumeID]; - - auto survive = [&] { - currentTrack.energy = energy; - currentTrack.pos = pos; - currentTrack.dir = dir; - currentTrack.navState = navState; - activeQueue->push_back(slot); - }; - - // Init a track with the needed data to call into G4HepEm. - G4HepEmElectronTrack elTrack; - G4HepEmTrack *theTrack = elTrack.GetTrack(); - theTrack->SetEKin(energy); - theTrack->SetMCIndex(theMCIndex); - theTrack->SetOnBoundary(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. - RanluxppDouble newRNG(currentTrack.rngState.BranchNoAdvance()); - - // Compute safety, needed for MSC step limit. - double safety = 0; - if (!navState.IsOnBoundary()) { - safety = Navigator::ComputeSafety(pos, navState); - } - theTrack->SetSafety(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); - } - - // Call G4HepEm to compute the physics step limit. - G4HepEmElectronManager::HowFar(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge); - - // Remember MSC values for the next step(s). - currentTrack.initialRange = mscData->fInitialRange; - currentTrack.dynamicRangeFactor = mscData->fDynamicRangeFactor; - currentTrack.tlimitMin = mscData->fTlimitMin; - - // Get result into variables. - double geometricalStepLengthFromPhysics = theTrack->GetGStepLength(); - // The phyiscal step length is the amount that the particle experiences - // which might be longer than the geometrical step length due to MSC. As - // long as we call PerformContinuous in the same kernel we don't need to - // care, but we need to make this available when splitting the operations. - // double physicalStepLength = elTrack.GetPStepLength(); - int winnerProcessIndex = theTrack->GetWinnerProcessIndex(); - // Leave the range and MFP inside the G4HepEmTrack. If we split kernels, we - // also need to carry them over! - - // Check if there's a volume boundary in between. - bool propagated = true; - double geometryStepLength; - vecgeom::NavStateIndex nextState; - if (BzFieldValue != 0) { - geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume( - energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, propagated, safety); - } else { - geometryStepLength = - Navigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, nextState); - pos += geometryStepLength * 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) - navState.SetBoundaryState(nextState.IsOnBoundary()); - - // Propagate information from geometrical step to MSC. - theTrack->SetDirection(dir.x(), dir.y(), dir.z()); - theTrack->SetGStepLength(geometryStepLength); - theTrack->SetOnBoundary(nextState.IsOnBoundary()); - - // Apply continuous effects. - bool stopped = G4HepEmElectronManager::PerformContinuous(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge); - - // Collect the direction change and displacement by MSC. - const double *direction = theTrack->GetDirection(); - dir.Set(direction[0], direction[1], direction[2]); - if (!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. - safety -= geometryStepLength; - constexpr double sFact = 0.99; - double reducedSafety = sFact * safety; - - // Apply displacement, depending on how close we are to a boundary. - // 1a. Far away from geometry boundary: - if (reducedSafety > 0.0 && dispR <= reducedSafety) { - pos += displacement; - } else { - // Recompute safety. - safety = Navigator::ComputeSafety(pos, navState); - reducedSafety = sFact * safety; - - // 1b. Far away from geometry boundary: - if (reducedSafety > 0.0 && dispR <= reducedSafety) { - pos += displacement; - // 2. Push to boundary: - } else if (reducedSafety > kGeomMinLength) { - pos += displacement * (reducedSafety / dispR); - } - // 3. Very small safety: do nothing. - } - } - } - - // Collect the charged step length (might be changed by MSC). - atomicAdd(&globalScoring->chargedSteps, 1); - atomicAdd(&scoringPerVolume->chargedTrackLength[volumeID], elTrack.GetPStepLength()); - - // Collect the changes in energy and deposit. - energy = theTrack->GetEKin(); - double energyDeposit = theTrack->GetEnergyDeposit(); - atomicAdd(&globalScoring->energyDeposit, energyDeposit); - atomicAdd(&scoringPerVolume->energyDeposit[volumeID], energyDeposit); - - // 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 (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(); - atomicAdd(&globalScoring->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(pos, navState); - newRNG.Advance(); - gamma1.rngState = newRNG; - gamma1.energy = copcore::units::kElectronMassC2; - gamma1.dir.Set(sint * cosPhi, sint * sinPhi, cost); - - gamma2.InitAsSecondary(pos, navState); - // Reuse the RNG state of the dying track. - gamma2.rngState = currentTrack.rngState; - gamma2.energy = copcore::units::kElectronMassC2; - gamma2.dir = -gamma1.dir; - } - // Particles are killed by not enqueuing them into the new activeQueue. - continue; - } - - if (nextState.IsOnBoundary()) { - // For now, just count that we hit something. - atomicAdd(&globalScoring->hits, 1); - - // Kill the particle if it left the world. - if (nextState.Top() != nullptr) { - Navigator::RelocateToNextVolume(pos, dir, nextState); - - // Move to the next boundary. - navState = nextState; - survive(); - } - continue; - } else if (!propagated) { - // Did not yet reach the interaction point due to error in the magnetic - // field propagation. Try again next time. - survive(); - continue; - } else if (winnerProcessIndex < 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[winnerProcessIndex] = -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. - 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[theMCIndex].fSecElProdCutE; - - switch (winnerProcessIndex) { - case 0: { - // Invoke ionization (for e-/e+): - double deltaEkin = (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, energy, &rnge) - : G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, energy, &rnge); - - double dirPrimary[] = {dir.x(), dir.y(), dir.z()}; - double dirSecondary[3]; - G4HepEmElectronInteractionIoni::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge); - - Track &secondary = secondaries.electrons.NextTrack(); - atomicAdd(&globalScoring->numElectrons, 1); - - secondary.InitAsSecondary(pos, navState); - secondary.rngState = newRNG; - secondary.energy = deltaEkin; - secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); - - energy -= deltaEkin; - 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(energy); - double deltaEkin = energy < g4HepEmPars.fElectronBremModelLim - ? G4HepEmElectronInteractionBrem::SampleETransferSB(&g4HepEmData, energy, logEnergy, - theMCIndex, &rnge, IsElectron) - : G4HepEmElectronInteractionBrem::SampleETransferRB(&g4HepEmData, energy, logEnergy, - theMCIndex, &rnge, IsElectron); - - double dirPrimary[] = {dir.x(), dir.y(), dir.z()}; - double dirSecondary[3]; - G4HepEmElectronInteractionBrem::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge); - - Track &gamma = secondaries.gammas.NextTrack(); - atomicAdd(&globalScoring->numGammas, 1); - - gamma.InitAsSecondary(pos, navState); - gamma.rngState = newRNG; - gamma.energy = deltaEkin; - gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]); - - energy -= deltaEkin; - dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]); - survive(); - break; - } - case 2: { - // Invoke annihilation (in-flight) for e+ - double dirPrimary[] = {dir.x(), dir.y(), dir.z()}; - double theGamma1Ekin, theGamma2Ekin; - double theGamma1Dir[3], theGamma2Dir[3]; - G4HepEmPositronInteractionAnnihilation::SampleEnergyAndDirectionsInFlight( - energy, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge); - - Track &gamma1 = secondaries.gammas.NextTrack(); - Track &gamma2 = secondaries.gammas.NextTrack(); - atomicAdd(&globalScoring->numGammas, 2); - - gamma1.InitAsSecondary(pos, navState); - gamma1.rngState = newRNG; - gamma1.energy = theGamma1Ekin; - gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]); - - gamma2.InitAsSecondary(pos, navState); - // Reuse the RNG state of the dying track. - gamma2.rngState = currentTrack.rngState; - gamma2.energy = theGamma2Ekin; - gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]); - - // The current track is killed by not enqueuing into the next activeQueue. - break; - } - } - } -} - -// Instantiate kernels for electrons and positrons. -__global__ void TransportElectrons(Track *electrons, const adept::MParray *active, Secondaries secondaries, - adept::MParray *activeQueue, GlobalScoring *globalScoring, - ScoringPerVolume *scoringPerVolume) -{ - TransportElectrons(electrons, active, secondaries, activeQueue, globalScoring, scoringPerVolume); -} -__global__ void TransportPositrons(Track *positrons, const adept::MParray *active, Secondaries secondaries, - adept::MParray *activeQueue, GlobalScoring *globalScoring, - ScoringPerVolume *scoringPerVolume) -{ - TransportElectrons(positrons, active, secondaries, activeQueue, globalScoring, - scoringPerVolume); -} diff --git a/test/TestEm3/gammas.cu b/test/TestEm3/gammas.cu deleted file mode 100644 index dceb902e..00000000 --- a/test/TestEm3/gammas.cu +++ /dev/null @@ -1,240 +0,0 @@ -// SPDX-FileCopyrightText: 2022 CERN -// SPDX-License-Identifier: Apache-2.0 - -#include "TestEm3.cuh" - -#include -#include -#ifdef ADEPT_USE_SURF -#include -#endif - -#include - -#include -#include -#include -#include -#include -#include -// Pull in implementation. -#include -#include -#include -#include - -__global__ void TransportGammas(Track *gammas, const adept::MParray *active, Secondaries secondaries, - adept::MParray *activeQueue, GlobalScoring *globalScoring, - ScoringPerVolume *scoringPerVolume) -{ -#ifdef ADEPT_USE_SURF - using Navigator = SurfNavigator; -#else - using Navigator = BVHNavigator; -#endif - int activeSize = active->size(); - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) { - const int slot = (*active)[i]; - Track ¤tTrack = gammas[slot]; - auto energy = currentTrack.energy; - auto pos = currentTrack.pos; - auto dir = currentTrack.dir; - auto navState = currentTrack.navState; - const int volumeID = navState.Top()->id(); - const int theMCIndex = MCIndex[volumeID]; - - auto survive = [&] { - currentTrack.energy = energy; - currentTrack.pos = pos; - currentTrack.dir = dir; - currentTrack.navState = navState; - activeQueue->push_back(slot); - }; - - // Init a track with the needed data to call into G4HepEm. - G4HepEmGammaTrack gammaTrack; - G4HepEmTrack *theTrack = gammaTrack.GetTrack(); - theTrack->SetEKin(energy); - theTrack->SetMCIndex(theMCIndex); - - // 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); - } - - // Call G4HepEm to compute the physics step limit. - G4HepEmGammaManager::HowFar(&g4HepEmData, &g4HepEmPars, &gammaTrack); - - // Get result into variables. - double geometricalStepLengthFromPhysics = theTrack->GetGStepLength(); - int winnerProcessIndex = theTrack->GetWinnerProcessIndex(); - // Leave the range and MFP inside the G4HepEmTrack. If we split kernels, we - // also need to carry them over! - - // Check if there's a volume boundary in between. - vecgeom::NavStateIndex nextState; - double geometryStepLength = - Navigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, nextState); - pos += geometryStepLength * dir; - atomicAdd(&globalScoring->neutralSteps, 1); - - // 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) - navState.SetBoundaryState(nextState.IsOnBoundary()); - - // Propagate information from geometrical step to G4HepEm. - theTrack->SetGStepLength(geometryStepLength); - theTrack->SetOnBoundary(nextState.IsOnBoundary()); - - 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; - } - - if (nextState.IsOnBoundary()) { - // For now, just count that we hit something. - atomicAdd(&globalScoring->hits, 1); - - // Kill the particle if it left the world. - if (nextState.Top() != nullptr) { - Navigator::RelocateToNextVolume(pos, dir, nextState); - - // Move to the next boundary. - navState = nextState; - survive(); - } - continue; - } else if (winnerProcessIndex < 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[winnerProcessIndex] = -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 (winnerProcessIndex) { - case 0: { - // Invoke gamma conversion to e-/e+ pairs, if the energy is above the threshold. - if (energy < 2 * copcore::units::kElectronMassC2) { - survive(); - continue; - } - - double logEnergy = std::log(energy); - double elKinEnergy, posKinEnergy; - G4HepEmGammaInteractionConversion::SampleKinEnergies(&g4HepEmData, energy, logEnergy, theMCIndex, elKinEnergy, - posKinEnergy, &rnge); - - double dirPrimary[] = {dir.x(), dir.y(), 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(); - atomicAdd(&globalScoring->numElectrons, 1); - atomicAdd(&globalScoring->numPositrons, 1); - - electron.InitAsSecondary(pos, navState); - electron.rngState = newRNG; - electron.energy = elKinEnergy; - electron.dir.Set(dirSecondaryEl[0], dirSecondaryEl[1], dirSecondaryEl[2]); - - positron.InitAsSecondary(pos, navState); - // Reuse the RNG state of the dying track. - positron.rngState = currentTrack.rngState; - positron.energy = 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 (energy < LowEnergyThreshold) { - survive(); - continue; - } - const double origDirPrimary[] = {dir.x(), dir.y(), dir.z()}; - double dirPrimary[3]; - const double newEnergyGamma = - G4HepEmGammaInteractionCompton::SamplePhotonEnergyAndDirection(energy, dirPrimary, origDirPrimary, &rnge); - vecgeom::Vector3D newDirGamma(dirPrimary[0], dirPrimary[1], dirPrimary[2]); - - const double energyEl = energy - newEnergyGamma; - if (energyEl > LowEnergyThreshold) { - // Create a secondary electron and sample/compute directions. - Track &electron = secondaries.electrons.NextTrack(); - atomicAdd(&globalScoring->numElectrons, 1); - - electron.InitAsSecondary(pos, navState); - electron.rngState = newRNG; - electron.energy = energyEl; - electron.dir = energy * dir - newEnergyGamma * newDirGamma; - electron.dir.Normalize(); - } else { - atomicAdd(&globalScoring->energyDeposit, energyEl); - atomicAdd(&scoringPerVolume->energyDeposit[volumeID], energyEl); - } - - // Check the new gamma energy and deposit if below threshold. - if (newEnergyGamma > LowEnergyThreshold) { - energy = newEnergyGamma; - dir = newDirGamma; - survive(); - } else { - atomicAdd(&globalScoring->energyDeposit, newEnergyGamma); - atomicAdd(&scoringPerVolume->energyDeposit[volumeID], newEnergyGamma); - // 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, theMCIndex, gammaTrack.GetPEmxSec(), energy, &rnge); - - double edep = bindingEnergy; - const double photoElecE = energy - edep; - if (photoElecE > theLowEnergyThreshold) { - // Create a secondary electron and sample directions. - Track &electron = secondaries.electrons.NextTrack(); - atomicAdd(&globalScoring->numElectrons, 1); - - double dirGamma[] = {dir.x(), dir.y(), dir.z()}; - double dirPhotoElec[3]; - G4HepEmGammaInteractionPhotoelectric::SamplePhotoElectronDirection(photoElecE, dirGamma, dirPhotoElec, &rnge); - - electron.InitAsSecondary(pos, navState); - electron.rngState = newRNG; - electron.energy = photoElecE; - electron.dir.Set(dirPhotoElec[0], dirPhotoElec[1], dirPhotoElec[2]); - } else { - edep = energy; - } - atomicAdd(&globalScoring->energyDeposit, edep); - atomicAdd(&scoringPerVolume->energyDeposit[volumeID], edep); - // The current track is killed by not enqueuing into the next activeQueue. - break; - } - } - } -}