WIP: Example12 - create a VecGeom geometry from Geant4

examples/Example12/CMakeLists.txt

@@ -0,0 +1,48 @@
+# SPDX-FileCopyrightText: 2021 CERN
+# SPDX-License-Identifier: Apache-2.0
+
+option(BUILTIN_G4VECGEOMNAV "Build G4VecGeomNav from source (requires VecGeom)" ON)
+
+if(NOT TARGET G4HepEm::g4HepEm)
+  message(STATUS "Disabling example12 (needs G4HepEm)")
+  return()
+endif()
+
+# Start  Configuration of G4VecGeomNav - converter from G4 to VecGeom volumes
+# -----  -------------    ------------
+# Minimum version of G4VecGeomNav we need.
+set(G4VecGeomNav_VERSION "") # 2021.06.18") # 0.1")
+
+if(BUILTIN_G4VECGEOMNAV)
+  include(FetchContent)
+
+  FetchContent_Declare(
+    g4vecgeomnav
+    GIT_REPOSITORY https://gitlab.cern.ch/VecGeom/g4vecgeomnav.git
+    GIT_TAG 1a002fbf2bccd81827b6100f86e2ad82433a75e9
+  )
+
+  FetchContent_GetProperties(g4vecgeomnav)
+
+  if(NOT g4vecgeomnav_POPULATED)
+    FetchContent_Populate(g4vecgeomnav)
+    message(STATUS "${g4vecgeomnav_SOURCE_DIR}, ${g4vecgeomnav_BINARY_DIR}")
+    add_subdirectory(${g4vecgeomnav_SOURCE_DIR} ${g4vecgeomnav_BINARY_DIR})
+  endif()
+else()
+  find_package(G4VecGeomNav REQUIRED)
+endif()
+# End   configuration of G4VecGeomNav 
+# ---   -------------    ------------
+
+add_executable(example12 example12.cpp example12.cu electrons.cu gammas.cu)
+target_link_libraries(example12
+   PRIVATE
+     AdePT CopCore::CopCore
+     G4VecGeomNav::g4vecgeomnav
+     VecGeom::vecgeom VecGeom::vecgeomcuda_static VecGeom::vgdml
+     ${Geant4_LIBRARIES}
+     G4HepEm::g4HepEmData G4HepEm::g4HepEmInit G4HepEm::g4HepEmRun)
+set_target_properties(example12 PROPERTIES CUDA_SEPARABLE_COMPILATION ON CUDA_RESOLVE_DEVICE_SYMBOLS ON)
+
+add_test(NAME example12 COMMAND example12 -gdml_name "${TESTING_GDML}")

examples/Example12/README.md

@@ -0,0 +1,50 @@
+<!--
+SPDX-FileCopyrightText: 2021 CERN
+SPDX-License-Identifier: CC-BY-4.0
+-->
+
+## Example 11
+
+Demonstrator of particle transportation on GPUs, with:
+
+ * geometry via VecGeom and a magnetic field with constant Bz,
+ * physics processes for e-/e+ and gammas using G4HepEm.
+ * accelerated navigation using newly added BVH class in VecGeom
+
+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.
+
+### Kernels
+
+This example uses one stream per particle type to launch kernels asynchronously.
+They are synchronized via a forth stream using CUDA events.
+
+#### `TransportElectrons<bool IsElectron>`
+
+1. Determine physics step limit.
+2. Call magnetic field to get geometry step length.
+3. Apply continuous effects; kill track if stopped.
+4. If the particle reaches a boundary, perform relocation.
+5. 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.

examples/Example12/electrons.cu

@@ -0,0 +1,263 @@
+// SPDX-FileCopyrightText: 2021 CERN
+// SPDX-License-Identifier: Apache-2.0
+
+#include "example11.cuh"
+
+#include <fieldPropagatorConstBz.h>
+
+#include <CopCore/PhysicalConstants.h>
+#include <AdePT/BVHNavigator.h>
+
+#include <G4HepEmElectronManager.hh>
+#include <G4HepEmElectronTrack.hh>
+#include <G4HepEmElectronInteractionBrem.hh>
+#include <G4HepEmElectronInteractionIoni.hh>
+#include <G4HepEmPositronInteractionAnnihilation.hh>
+// Pull in implementation.
+#include <G4HepEmRunUtils.icc>
+#include <G4HepEmInteractionUtils.icc>
+#include <G4HepEmElectronManager.icc>
+#include <G4HepEmElectronInteractionBrem.icc>
+#include <G4HepEmElectronInteractionIoni.icc>
+#include <G4HepEmPositronInteractionAnnihilation.icc>
+
+// 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 <bool IsElectron>
+static __device__ __forceinline__ void TransportElectrons(Track *electrons, const adept::MParray *active,
+                                                          Secondaries &secondaries, adept::MParray *activeQueue,
+                                                          GlobalScoring *scoring)
+{
+  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 &currentTrack = electrons[slot];
+
+    // Init a track with the needed data to call into G4HepEm.
+    G4HepEmElectronTrack elTrack;
+    G4HepEmTrack *theTrack = elTrack.GetTrack();
+    theTrack->SetEKin(currentTrack.energy);
+
+    // Fetch the material cuts-couple index from the G4 track
+    auto cur_state= track->currentState();
+    // auto cur_nav_index= cur_state.GetNavIndex(cur_state.GetLevel());
+    VPlacedVolume const * placedVol = cur_state.Top();
+    LogicalVolume const * logicalVol = placedVol->GetLogicalVolume();
+    int theMCIndex = (int) logicalVol->GetMaterialCutsPtr();  // It's not the pointer, it's an index !!
+    std::cout << " Volume (log): " << logicalVol->GetName() << " mat-cuts-cpl = " << theMCindex << "\n";
+
+    theTrack->SetMCIndex(theMCIndex);
+
+    theTrack->SetCharge(Charge);
+
+    // 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.
+    G4HepEmElectronManager::HowFar(&g4HepEmData, &g4HepEmPars, &elTrack);
+
+    // 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.
+    vecgeom::NavStateIndex nextState;
+    double geometryStepLength = fieldPropagatorBz.ComputeStepAndPropagatedState</*Relocate=*/false, BVHNavigator>(
+        currentTrack.energy, Mass, Charge, geometricalStepLengthFromPhysics, currentTrack.pos, currentTrack.dir,
+        currentTrack.navState, nextState);
+
+    if (nextState.IsOnBoundary()) {
+      theTrack->SetGStepLength(geometryStepLength);
+      theTrack->SetOnBoundary(true);
+    }
+
+    // Apply continuous effects.
+    bool stopped = G4HepEmElectronManager::PerformContinuous(&g4HepEmData, &g4HepEmPars, &elTrack);
+    // Collect the changes.
+    currentTrack.energy = theTrack->GetEKin();
+    atomicAdd(&scoring->energyDeposit, theTrack->GetEnergyDeposit());
+
+    // 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(&scoring->secondaries, 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(/*parent=*/currentTrack);
+        gamma1.rngState = currentTrack.rngState.Branch();
+        gamma1.energy = copcore::units::kElectronMassC2;
+        gamma1.dir.Set(sint * cosPhi, sint * sinPhi, cost);
+
+        gamma2.InitAsSecondary(/*parent=*/currentTrack);
+        // 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(&scoring->hits, 1);
+
+      // Kill the particle if it left the world.
+      if (nextState.Top() != nullptr) {
+        activeQueue->push_back(slot);
+
+        // Move to the next boundary.
+        BVHNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, nextState);
+        currentTrack.navState = nextState;
+      }
+      continue;
+    } else if (winnerProcessIndex < 0) {
+      // No discrete process, move on.
+      activeQueue->push_back(slot);
+      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.
+      activeQueue->push_back(slot);
+      continue;
+    }
+
+    // Perform the discrete interaction.
+    RanluxppDoubleEngine rnge(&currentTrack.rngState);
+    // We will need one branched RNG state, prepare while threads are synchronized.
+    RanluxppDouble newRNG(currentTrack.rngState.Branch());
+
+    const double energy   = currentTrack.energy;
+    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[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
+      double dirSecondary[3];
+      G4HepEmElectronInteractionIoni::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
+
+      Track &secondary = secondaries.electrons.NextTrack();
+      atomicAdd(&scoring->secondaries, 1);
+
+      secondary.InitAsSecondary(/*parent=*/currentTrack);
+      secondary.rngState = newRNG;
+      secondary.energy = deltaEkin;
+      secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
+
+      currentTrack.energy = energy - deltaEkin;
+      currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
+      // The current track continues to live.
+      activeQueue->push_back(slot);
+      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[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
+      double dirSecondary[3];
+      G4HepEmElectronInteractionBrem::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
+
+      Track &gamma = secondaries.gammas.NextTrack();
+      atomicAdd(&scoring->secondaries, 1);
+
+      gamma.InitAsSecondary(/*parent=*/currentTrack);
+      gamma.rngState = newRNG;
+      gamma.energy = deltaEkin;
+      gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
+
+      currentTrack.energy = energy - deltaEkin;
+      currentTrack.dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
+      // The current track continues to live.
+      activeQueue->push_back(slot);
+      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(
+          energy, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge);
+
+      Track &gamma1 = secondaries.gammas.NextTrack();
+      Track &gamma2 = secondaries.gammas.NextTrack();
+      atomicAdd(&scoring->secondaries, 2);
+
+      gamma1.InitAsSecondary(/*parent=*/currentTrack);
+      gamma1.rngState = newRNG;
+      gamma1.energy = theGamma1Ekin;
+      gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]);
+
+      gamma2.InitAsSecondary(/*parent=*/currentTrack);
+      // 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 *scoring)
+{
+  TransportElectrons</*IsElectron*/ true>(electrons, active, secondaries, activeQueue, scoring);
+}
+__global__ void TransportPositrons(Track *positrons, const adept::MParray *active, Secondaries secondaries,
+                                   adept::MParray *activeQueue, GlobalScoring *scoring)
+{
+  TransportElectrons</*IsElectron*/ false>(positrons, active, secondaries, activeQueue, scoring);
+}