Field Propagation using Runge-Kutta integration and example15 as demo

examples/CMakeLists.txt

@@ -17,6 +17,7 @@ add_subdirectory(Example11)
 add_subdirectory(Example12)
 add_subdirectory(Example13)
 add_subdirectory(Example14)
+add_subdirectory(Example15)
 add_subdirectory(Example16)
 add_subdirectory(Example17)
 add_subdirectory(Example18)

examples/Example11/README.md

@@ -19,7 +19,7 @@ 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.
+They are synchronized via a fourth stream using CUDA events.
 
 #### `TransportElectrons<bool IsElectron>`
 

examples/Example12/README.md

@@ -22,7 +22,7 @@ Additionally, the kernels score energy deposit and the charged track length per
 ### Kernels
 
 This example uses one stream per particle type to launch kernels asynchronously.
-They are synchronized via a forth stream using CUDA events.
+They are synchronized via a fourth stream using CUDA events.
 
 #### `TransportElectrons<bool IsElectron>`
 

examples/Example13/README.md

@@ -24,7 +24,7 @@ Additionally, the kernels score energy deposit and the charged track length per
 ### Kernels
 
 This example uses one stream per particle type to launch kernels asynchronously.
-They are synchronized via a forth stream using CUDA events.
+They are synchronized via a fourth stream using CUDA events.
 
 #### `TransportElectrons<bool IsElectron>`
 

examples/Example13/electrons.cu

@@ -123,9 +123,12 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
     bool propagated = true;
     double geometryStepLength;
     vecgeom::NavStateIndex nextState;
+    constexpr int max_iters = 10;
+    // printf("max_iters= %3d\n", max_iters);
     if (BzFieldValue != 0) {
       geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume<BVHNavigator>(
-          energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, propagated, safety);
+          energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState,
+          propagated, safety, max_iters );
     } else {
       geometryStepLength = BVHNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState,
                                                                   nextState, kPush);

examples/Example13/example13.cuh

@@ -135,7 +135,7 @@ extern __constant__ __device__ struct G4HepEmData g4HepEmData;
 
 extern __constant__ __device__ int *MCIndex;
 
-// constexpr vecgeom::Precision BzFieldValue = 0.1 * copcore::units::tesla;
-constexpr vecgeom::Precision BzFieldValue = 0;
+constexpr vecgeom::Precision BzFieldValue = 3.8 * copcore::units::tesla;
+// constexpr vecgeom::Precision BzFieldValue = 0;
 
 #endif

examples/Example15/CMakeLists.txt

@@ -0,0 +1,42 @@
+# SPDX-FileCopyrightText: 2021 CERN
+# SPDX-License-Identifier: Apache-2.0
+
+if(NOT TARGET G4HepEm::g4HepEm)
+  message(STATUS "Disabling example15 (needs G4HepEm)")
+  return()
+endif()
+
+if(Geant4_FOUND)
+  if(NOT Geant4_gdml_FOUND)
+    message(STATUS "Disabling example15 (needs Geant4 with GDML support)")
+    return()
+  endif()
+else()
+  message(STATUS "Disabling example15 (needs Geant4)")
+  return()
+endif()
+
+# example15 is the AdePT demo example using material cuts as defined in the input gdml file
+add_executable(example15
+  example15.cpp
+  example15.cu
+  electrons.cu
+  gammas.cu)
+target_link_libraries(example15
+  PRIVATE
+    AdePT
+    CopCore::CopCore
+    VecGeom::vecgeom
+    VecGeom::vecgeomcuda_static
+    VecGeom::vgdml
+    ${Geant4_LIBRARIES}
+    G4HepEm::g4HepEmData
+    G4HepEm::g4HepEmInit
+    G4HepEm::g4HepEmRun
+    CUDA::cudart)
+
+set_target_properties(example15 PROPERTIES CUDA_SEPARABLE_COMPILATION ON CUDA_RESOLVE_DEVICE_SYMBOLS ON)
+
+# Tests
+add_test(NAME example15
+  COMMAND $<TARGET_FILE:example15> -gdml_file ${CMAKE_BINARY_DIR}/cms2018.gdml)

examples/Example15/README.md

@@ -0,0 +1,64 @@
+<!--
+SPDX-FileCopyrightText: 2021 CERN
+SPDX-License-Identifier: CC-BY-4.0
+-->
+
+## Example 14
+
+Example demonstrating particle transportation on GPUs with a non-uniform magnetic field - in arbitrary geometry read from a GDML file.
+
+New feature is
+ - definition of a non-uniform magnetic field
+ - integration of the equation of motion of charged particles using embedded Runge-Kutta method, Dormand Prince 4(5).
+ Caveats (at present) : 
+   - fixed accuracy of integration (constant in source)
+
+The features carried over from Example13 are:
+ * arbitrary geometry via gdml file (tested with cms2018.gdml from VecGeom persistency/gdml/gdmls folder) and optionally a magnetic field with constant Bz,
+ * geometry read as Geant4 geometry, reading in regions and cuts, to initialize G4HepEm data
+ * geometry read then into VecGeom, and synchronized to GPU
+ * G4HepEm material-cuts couple indices mapped to VecGeom logical volume id's
+ * physics processes for e-/e+ (including MSC) and gammas using G4HepEm.
+ * scoring per placed volume, no sensitive detector feature
+
+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<bool IsElectron>`
+
+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 discplacement.
+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.