Update gamma kernel, add reproducibility CI test (#333)

This updates the gamma kernel, which is a leftover from #324.

The `SampleInteraction` only needs to be called if no boundary is hit
(as otherwise no physics process is used). This enables to keep the
`numIALeft` which is stored in the `currentTrack` and needs to be passed
to the G4HepEm theTrack. The previous early exit with setting the
winnerprocess to -1 caused non-reproducibility.

Since the code is reproducible with this fix, the CI test for
reproducibility with 8 ttbar events in CMS (no B field) is activated and
replaced the previous test that just ran one event.
Physics validation results are just as before.

examples/CMakeLists.txt

@@ -9,4 +9,4 @@ configure_file(data/ppttbar.hepmc3 ${PROJECT_BINARY_DIR}/ppttbar.hepmc3)
 # - Subprojects
 add_subdirectory(Example1)
 add_subdirectory(IntegrationBenchmark)
-add_subdirectory(AsyncExample)
+# add_subdirectory(AsyncExample)

examples/IntegrationBenchmark/CMakeLists.txt

@@ -68,18 +68,12 @@ target_link_libraries(integrationBenchmark
 # Scripts
 configure_file("macros/integrationbenchmark.mac.in" "${CMAKE_BINARY_DIR}/integrationbenchmark.mac")
 
-# Tests
-add_test(NAME integrationBenchmark
-  COMMAND $<TARGET_FILE:integrationBenchmark> -m ${PROJECT_BINARY_DIR}/example1_large_stack.mac
-)
-
 # This test checks the reproducibility of AdePT by running 8 ttbar events and checking that the energy deposition is exactly the same.
-# NOTE: Currently not active since the results are currently NOT reproducible!! To replace the integrationBenchmark test above.
-# add_test(NAME reproducibility_cms_ttbar
-#     COMMAND bash ${PROJECT_SOURCE_DIR}/examples/IntegrationBenchmark/ci_tests/reproducibility.sh
-#     "$<TARGET_FILE:integrationBenchmark>" "${PROJECT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}"
-#     WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
-# )
+add_test(NAME reproducibility_cms_ttbar
+    COMMAND bash ${PROJECT_SOURCE_DIR}/examples/IntegrationBenchmark/ci_tests/reproducibility.sh
+    "$<TARGET_FILE:integrationBenchmark>" "${PROJECT_BINARY_DIR}" "${PROJECT_SOURCE_DIR}"
+    WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+)
 
 # test that compares the physics output of a full AdePT run against a high-statistics Geant4 simulation using G4HepEm.
 # The energy deposition per layer error must be < 1% to pass the test

examples/IntegrationBenchmark/ci_tests/python_scripts/check_validation.py

@@ -80,12 +80,12 @@ def compare_csv(file1, file2, n1, n2, tol=0.01, plot_file=None):
 
     # Print results
     if failed_layers:
-        print(f"The results are not reproducible. Relative errors exceed {100*tol}% in the following layers:")
+        print(f"The physics results are not valid. Relative errors exceed {100*tol}% in the following layers:")
         for layer, err, val1, val2 in failed_layers:
             print(f"Layer {layer}: Relative Error = {err:.6f}%, File1 = {val1:.6f}, File2 = {val2:.6f}")
         sys.exit(1)
     else:
-        print(f"The results are reproducible. All layers have relative errors within {100*tol}%.")
+        print(f"The physics results are valid. All layers have relative errors within {100*tol}%.")
 
 
 if __name__ == "__main__":

examples/IntegrationBenchmark/ci_tests/reproducibility.sh

@@ -34,11 +34,11 @@ cleanup
 $CI_TEST_DIR/python_scripts/macro_generator.py \
     --template ${CI_TEST_DIR}/example_template.mac \
     --output ${CI_TEST_DIR}/reproducibility.mac \
-    --gdml_name ${PROJECT_BINARY_DIR}/cms2018_sd.gdml \
-    --num_threads 24 \
-    --num_events 24 \
-    --num_trackslots 16 \
-    --num_hitslots 12 \
+    --gdml_name ${PROJECT_SOURCE_DIR}/examples/data/cms2018_sd.gdml \
+    --num_threads 4 \
+    --num_events 8 \
+    --num_trackslots 10 \
+    --num_hitslots 4 \
     --gun_type hepmc \
     --event_file ${PROJECT_BINARY_DIR}/ppttbar.hepmc3
 

include/AdePT/kernels/gammas.cuh

@@ -76,24 +76,17 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
     theTrack->SetEKin(eKin);
     theTrack->SetMCIndex(auxData.fMCIndex);
 
-    // Sample the `number-of-interaction-left` and put it into the track.
-    // Use index 0 since numIALeft for gammas is based only on the total macroscopic cross section
-    if (theTrack->GetNumIALeft(0) <= 0.0) {
+    // Re-sample the `number-of-interaction-left` (if needed, otherwise use stored numIALeft) and put it into the
+    // G4HepEmTrack. Use index 0 since numIALeft for gammas is based only on the total macroscopic cross section. The
+    // currentTrack.numIALeft[0] are updated later
+    if (currentTrack.numIALeft[0] <= 0.0) {
       theTrack->SetNumIALeft(-std::log(currentTrack.Uniform()), 0);
+    } else {
+      theTrack->SetNumIALeft(currentTrack.numIALeft[0], 0);
     }
 
     // Call G4HepEm to compute the physics step limit.
     G4HepEmGammaManager::HowFar(&g4HepEmData, &g4HepEmPars, &gammaTrack);
-    G4HepEmGammaManager::SampleInteraction(&g4HepEmData, &gammaTrack, currentTrack.Uniform());
-    int winnerProcessIndex = theTrack->GetWinnerProcessIndex();
-
-    // avoid photo-nuclear reaction that would need to be handled by G4 itself
-    if (winnerProcessIndex == 3) {
-      winnerProcessIndex = -1;
-      // NOTE: no simple re-drawing is possible, since HowFar returns now smaller steps due to the gamma-nuclear
-      // reactions in comparison to without gamma-nuclear reactions. Thus, an empty step without a reaction is needed to
-      // compensate for the smaller step size returned by HowFar.
-    }
 
     // Get result into variables.
     double geometricalStepLengthFromPhysics = theTrack->GetGStepLength();
@@ -126,18 +119,20 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
     theTrack->SetGStepLength(geometryStepLength);
     theTrack->SetOnBoundary(nextState.IsOnBoundary());
 
-    G4HepEmGammaManager::UpdateNumIALeft(theTrack);
-
-    // Save the `number-of-interaction-left` in our track.
-    // Use index 0 since numIALeft for gammas is based only on the total macroscopic cross section
-    double numIALeft          = theTrack->GetNumIALeft(0);
-    currentTrack.numIALeft[0] = numIALeft;
-
+    int winnerProcessIndex;
     if (nextState.IsOnBoundary()) {
       // For now, just count that we hit something.
 
       // Kill the particle if it left the world.
       if (!nextState.IsOutside()) {
+
+        G4HepEmGammaManager::UpdateNumIALeft(theTrack);
+
+        // Save the `number-of-interaction-left` in our track.
+        // Use index 0 since numIALeft stores for gammas only the total macroscopic cross section
+        double numIALeft          = theTrack->GetNumIALeft(0);
+        currentTrack.numIALeft[0] = numIALeft;
+
 #ifdef ADEPT_USE_SURF
         AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState);
         if (nextState.IsOutside()) continue;
@@ -164,15 +159,15 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
         }
       }
       continue;
-    } else if (winnerProcessIndex < 0) {
-      // No discrete process, move on.
-      survive();
-      continue;
-    }
+    } else {
+
+      G4HepEmGammaManager::SampleInteraction(&g4HepEmData, &gammaTrack, currentTrack.Uniform());
+      winnerProcessIndex = theTrack->GetWinnerProcessIndex();
 
-    // Reset number of interaction left for the winner discrete process.
-    // (Will be resampled in the next iteration.)
-    currentTrack.numIALeft[winnerProcessIndex] = -1.0;
+      // Reset number of interaction left for the winner discrete process also in the currentTrack (SampleInteraction()
+      // resets it for theTrack), will be resampled in the next iteration.
+      currentTrack.numIALeft[0] = -1.0;
+    }
 
     // Update the flight times of the particle
     double deltaTime = theTrack->GetGStepLength() / copcore::units::kCLight;
@@ -347,6 +342,11 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries
       // The current track is killed by not enqueuing into the next activeQueue.
       break;
     }
+    case 3: {
+      // Invoke gamma nuclear needs to be handled by Geant4 directly, to be implemented
+      // Just keep particle alive
+      survive();
+    }
     }
   }
 }