diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 8243d81b..d900b5fe 100644 --- a/examples/CMakeLists.txt +++ b/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) diff --git a/examples/IntegrationBenchmark/CMakeLists.txt b/examples/IntegrationBenchmark/CMakeLists.txt index 81fbea7d..e00f9117 100644 --- a/examples/IntegrationBenchmark/CMakeLists.txt +++ b/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 $ -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 -# "$" "${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 + "$" "${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 diff --git a/examples/IntegrationBenchmark/ci_tests/python_scripts/check_validation.py b/examples/IntegrationBenchmark/ci_tests/python_scripts/check_validation.py index ddfdacc8..d5f8a7ca 100755 --- a/examples/IntegrationBenchmark/ci_tests/python_scripts/check_validation.py +++ b/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__": diff --git a/examples/IntegrationBenchmark/ci_tests/reproducibility.sh b/examples/IntegrationBenchmark/ci_tests/reproducibility.sh index ede1a05c..a1c729df 100644 --- a/examples/IntegrationBenchmark/ci_tests/reproducibility.sh +++ b/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 diff --git a/include/AdePT/kernels/gammas.cuh b/include/AdePT/kernels/gammas.cuh index 6e1d2c99..787dcabd 100644 --- a/include/AdePT/kernels/gammas.cuh +++ b/include/AdePT/kernels/gammas.cuh @@ -76,24 +76,17 @@ __global__ void TransportGammas(adept::TrackManager *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 *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 *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 *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(); + } } } }