diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index aa52dd3f38b47442e0c1d3697735b2306dfc0d22..53f2e7b4ce54b88a9e45df329eb9a7348829330c 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -16,6 +16,7 @@ variables: # _alternatively_ corsika-data can be downloaded as submodule: GIT_SUBMODULE_STRATEGY: normal # none: we get the submodules in before_script, # normal: get submodules automatically + CTEST_OUTPUT_ON_FAILURE: 1 # @@ -32,6 +33,7 @@ stages: - build_test - example - build_test_example + - python - install - optional @@ -251,6 +253,8 @@ test-clang-8: artifacts: when: always expire_in: 3 days + paths: + - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml reports: junit: - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml @@ -311,6 +315,7 @@ build_test-clang-8: - if: $CI_COMMIT_TAG when: manual allow_failure: true + allow_failure: true artifacts: when: always expire_in: 3 days @@ -319,6 +324,7 @@ build_test-clang-8: - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml paths: - ${CI_PROJECT_DIR}/build/build_examples/examples.log.gz + - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml cache: paths: - ${CI_PROJECT_DIR}/build/ @@ -440,6 +446,7 @@ install-clang-8: - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml paths: - ${CI_PROJECT_DIR}/build/build_examples/examples.log.gz + - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml # release for gcc release-full-u-18_04: @@ -536,3 +543,38 @@ sanity: untracked: true policy: pull key: "${CI_COMMIT_REF_SLUG}-gcc" + +########################################################## +# template for all Python jobs +.python: + stage: python + tags: + - corsika + before_script: + - apt-get update && apt-get install -y python3-pip + script: + - cd ${CI_PROJECT_DIR}/python # change into the Python directory + - pip3 install --user -e '.[test]' # install the package + test deps + - python3 --version + - python3 -m mypy corsika + - python3 -m isort --atomic --check-only corsika tests + - python3 -m black -t py37 corsika tests + - python3 -m flake8 corsika tests + - python3 -m pytest --cov=corsika tests + - cd ${CI_PROJECT_DIR} # reset the directory + coverage: '/^TOTAL\s*\d+\s*\d+\s*(.*\%)/' + +# the default Python version Ubuntu 18.04 is Python3.8 +python-3.8: + extends: .python + image: corsika/devel:u-18.04 + dependencies: + - build-u-18_04 + cache: + key: "${CI_COMMIT_REF_SLUG}-gcc" + artifacts: + when: always + expire_in: 1 year + paths: + - ${CI_PROJECT_DIR}/Python/python-test.log + allow_failure: true diff --git a/.gitmodules b/.gitmodules index 1bbe5d4d093457574b332f922f45b43b2d2f8b53..5429f2d75f591d52f001615fddde188e8bd012de 100644 --- a/.gitmodules +++ b/.gitmodules @@ -12,4 +12,4 @@ path = modules/conex url = ../../AirShowerPhysics/cxroot.git branch = master - shallow = true + shallow = true diff --git a/CMakeLists.txt b/CMakeLists.txt index 3cd6f0ba077e877aa7c1802aa807f9481a073c75..d81de38f705fa7f6b5901a5643e74cb2f5f9ab93 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,3 +1,13 @@ +# +# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu +# +# See file AUTHORS for a list of contributors. +# +# This software is distributed under the terms of the GNU General Public +# Licence version 3 (GPL Version 3). See file LICENSE for a full version of +# the license. +# + cmake_minimum_required (VERSION 3.9) #+++++++++++++++++++++++++++++ @@ -32,12 +42,6 @@ if (POLICY CMP0079) cmake_policy (SET CMP0079 NEW) endif () -#+++++++++++++++++++++++++++++ -# as long as there still are modules using it: -# -enable_language (Fortran) -set (CMAKE_Fortran_FLAGS "-std=legacy -Wfunction-elimination") - #+++++++++++++++++++++++++++++ # warn user if system is not UNIX # @@ -50,15 +54,19 @@ endif () # set (CORSIKA8_CMAKE_DIR "${PROJECT_SOURCE_DIR}/cmake") set (CMAKE_MODULE_PATH "${CORSIKA8_CMAKE_DIR}" ${CMAKE_MODULE_PATH}) -include (CorsikaUtilities) # extra cmake function set (CMAKE_VERBOSE_MAKEFILE OFF) # this can be done with `make VERBOSE=1` # ignore many irrelevant Up-to-date messages during install set (CMAKE_INSTALL_MESSAGE LAZY) #+++++++++++++++++++++++++++++ -# Setup hardware and infrastructure dependent defines +# Extra cmake functions for registering and running tests +# +include (corsikaUtilities) # extra cmake functions + +#+++++++++++++++++++++++++++++ +# Setup hardware and infrastructure dependent defines, and general settings # -include (CorsikaDefines) +include (corsikaDefines) #+++++++++++++++++++++++++++++ # check if compiler is C++17 compliant @@ -81,50 +89,14 @@ set (CMAKE_CXX_EXTENSIONS OFF) # debug: O0, relwithdebinfo: 02, release: O3, minsizerel: Os (all defaults) set (CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers") if (CORSIKA_SCL_CXX) - add_definitions(-D_GLIBCXX_USE_CXX11_ABI=0) + add_definitions (-D_GLIBCXX_USE_CXX11_ABI=0) endif() -set (DEFAULT_BUILD_TYPE "Release") # clang produces a lot of unecessary warnings without this: add_compile_options ( $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>>:-Wno-nonportable-include-path> ) -#+++++++++++++++++++++++++++++ -# Build types settings -# -# setup coverage build type -set (CMAKE_CXX_FLAGS_COVERAGE "-g --coverage") -set (CMAKE_EXE_LINKER_FLAGS_COVERAGE "--coverage") -set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage") -# set a flag to inform code that we are in debug mode -set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG") - -#+++++++++++++++++++++++++++++ -# Build type selection -# -# Set the possible values of build type for cmake-gui and command line check -set (ALLOWED_BUILD_TYPES Debug Release MinSizeRel RelWithDebInfo Coverage) -set_property (CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS ${ALLOWED_BUILD_TYPES}) -# Set a default build type if none was specified -# by default: "Debug", if local ".git" directory is found, otherwise "Release" -set (DEFAULT_BUILD_TYPE "Release") -if (EXISTS "${CMAKE_SOURCE_DIR}/.git") - set (DEFAULT_BUILD_TYPE "Debug") -endif () -if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) - message (STATUS "Setting build type to '${DEFAULT_BUILD_TYPE}' as no other was specified.") - set (CMAKE_BUILD_TYPE "${DEFAULT_BUILD_TYPE}" CACHE - STRING "Choose the type of build." FORCE) -else (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) - # Ignore capitalization when build type is selected manually and check for valid setting - string (TOLOWER ${CMAKE_BUILD_TYPE} SELECTED_LOWER) - string (TOLOWER "${ALLOWED_BUILD_TYPES}" BUILD_TYPES_LOWER) - if (NOT SELECTED_LOWER IN_LIST BUILD_TYPES_LOWER) - message (FATAL_ERROR "Unknown build type: ${CMAKE_BUILD_TYPE} [allowed: ${ALLOWED_BUILD_TYPES}]") - endif () - message (STATUS "Build type is: ${CMAKE_BUILD_TYPE}") -endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) #+++++++++++++++++++++++++++++ # Setup external dependencies. @@ -133,6 +105,7 @@ endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) include (conan) # self-provided in 'cmake' directory # # download and build all dependencies +message (STATUS "Installing dependencies from Conan (this may take some time)...") conan_cmake_run (CONANFILE conanfile.txt BASIC_SETUP CMAKE_TARGETS BUILD missing @@ -238,6 +211,8 @@ target_link_libraries ( CONAN_PKG::eigen CONAN_PKG::spdlog CONAN_PKG::boost + CONAN_PKG::yaml-cpp + CONAN_PKG::arrow cnpy # for SaveBoostHistogram ) @@ -314,6 +289,13 @@ configure_package_config_file ( INSTALL_DESTINATION ${CMAKE_INSTALL_PREFIX} ) +# corsikaDefines for build tree +configure_package_config_file ( + cmake/corsikaDefines.cmake + ${PROJECT_BINARY_DIR}/corsikaDefines.cmake + INSTALL_DESTINATION ${CMAKE_INSTALL_PREFIX} + ) + # config for install tree # overwrite with install location (if it was build as part of corsika) set (Pythia8_INCDIR ${Pythia8_INCDIR_INSTALL}) @@ -329,6 +311,7 @@ install ( ${CMAKE_BINARY_DIR}/conanbuildinfo.cmake ${CMAKE_BINARY_DIR}/conaninfo.txt ${CMAKE_BINARY_DIR}/corsikaConfig.cmake + ${CMAKE_BINARY_DIR}/corsikaDefines.cmake ${CMAKE_BINARY_DIR}/corsikaConfigVersion.cmake DESTINATION lib/cmake/corsika ) diff --git a/MAINTAINERS.md b/MAINTAINERS.md index a57e4bfbbbff4526527819a3f2e4840c58a9b799..2da518cf5b6b9a8e0e3a278bf166fd246137554d 100644 --- a/MAINTAINERS.md +++ b/MAINTAINERS.md @@ -23,6 +23,13 @@ Hadron models: - Felix Riehn <friehn@lip.pt>, Santiago/Lisbon - Anatoli Fedynitch <anatoli.fedynitch@icecube.wisc.edu> ICRR Tokyo +Output formats and infrastructure: +- Remy Prechelt <prechelt@hawaii.edu>, UHM +- Ralf Ulrich <ralf.ulrich@kit.edu>, KIT + +Python library: +- Remy Prechelt <prechelt@hawaii.edu>, UHM + Radio: - Remy Prechelt <prechelt@hawaii.edu> - Tim Huege <tim.huege@kit.edu>, KIT diff --git a/cmake/CorsikaDefines.cmake b/cmake/CorsikaDefines.cmake deleted file mode 100644 index 2019ce2f906f9025e2181a0c9eb624e99f9ca44e..0000000000000000000000000000000000000000 --- a/cmake/CorsikaDefines.cmake +++ /dev/null @@ -1,29 +0,0 @@ -# -# Floating point exception support - select implementation to use -# -include (CheckIncludeFileCXX) -CHECK_INCLUDE_FILE_CXX ("fenv.h" HAS_FEENABLEEXCEPT) -if (HAS_FEENABLEEXCEPT) # FLOATING_POINT_ENVIRONMENT - set (CORSIKA_HAS_FEENABLEEXCEPT 1) - set_property (DIRECTORY ${CMAKE_HOME_DIRECTORY} APPEND PROPERTY COMPILE_DEFINITIONS "HAS_FEENABLEEXCEPT") -endif () - - -# -# General OS Detection -# -if (${CMAKE_SYSTEM_NAME} MATCHES "Windows") - set (CORSIKA_OS_WINDOWS TRUE) - set (CORSIKA_OS "Windows") -elseif (${CMAKE_SYSTEM_NAME} MATCHES "Linux") - set (CORSIKA_OS_LINUX TRUE) - set (CORSIKA_OS "Linux") - # check for RedHat/CentOS compiler from the software collections (scl) - string (FIND ${CMAKE_CXX_COMPILER} "/opt/rh/devtoolset-" index) - if (${index} EQUAL 0) - set (CORSIKA_SCL_CXX TRUE) - endif () -elseif (${CMAKE_SYSTEM_NAME} MATCHES "Darwin") - set (CORSIKA_OS_MAC TRUE) - set (CORSIKA_OS "Mac") -endif () diff --git a/CMakeModules/FindCONEX.cmake b/cmake/FindCONEX.cmake similarity index 100% rename from CMakeModules/FindCONEX.cmake rename to cmake/FindCONEX.cmake diff --git a/cmake/FindPhysUnits.cmake b/cmake/FindPhysUnits.cmake index 321ec043e882c5c98d611a6daf5c23438e8b4c81..20fa1aaa4928ec0f37ce8ee573c7a4659756f98b 100644 --- a/cmake/FindPhysUnits.cmake +++ b/cmake/FindPhysUnits.cmake @@ -1,3 +1,13 @@ +# +# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu +# +# See file AUTHORS for a list of contributors. +# +# This software is distributed under the terms of the GNU General Public +# Licence version 3 (GPL Version 3). See file LICENSE for a full version of +# the license. +# + add_library (PhysUnits INTERFACE) target_compile_options (PhysUnits diff --git a/cmake/FindSphinx.cmake b/cmake/FindSphinx.cmake index 9a6705c4840e57d3c7940289be5e498ad0fcce15..62825ee4c3d14139ca371c9dd8d8f028da25181e 100644 --- a/cmake/FindSphinx.cmake +++ b/cmake/FindSphinx.cmake @@ -1,3 +1,13 @@ +# +# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu +# +# See file AUTHORS for a list of contributors. +# +# This software is distributed under the terms of the GNU General Public +# Licence version 3 (GPL Version 3). See file LICENSE for a full version of +# the license. +# + # Look for an executable called sphinx-build find_program (SPHINX_EXECUTABLE NAMES sphinx-build diff --git a/cmake/corsikaConfig.cmake.in b/cmake/corsikaConfig.cmake.in index 546d5141d2c467cf7293bddf4cabef6fa28d481f..3a833b71c2ad78dcd0d8af5b194bd004774407f6 100644 --- a/cmake/corsikaConfig.cmake.in +++ b/cmake/corsikaConfig.cmake.in @@ -1,37 +1,27 @@ +# +# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu +# +# See file AUTHORS for a list of contributors. +# +# This software is distributed under the terms of the GNU General Public +# Licence version 3 (GPL Version 3). See file LICENSE for a full version of +# the license. +# + set (CORSIKA8_VERSION @c8_version@) @PACKAGE_INIT@ +#+++++++++++++++++++++++++++++ +# Setup hardware and infrastructure dependent defines and other setting +# +include (${CMAKE_CURRENT_LIST_DIR}/corsikaDefines.cmake) + #+++++++++++++++++++++++++++ # Options # option (WITH_HISTORY "Flag to switch on/off HISTORY" ON) - - -#+++++++++++++++++++++++++++++ -# Build types settings -# -# setup coverage build type -set (CMAKE_CXX_FLAGS_COVERAGE "-g --coverage") -set (CMAKE_EXE_LINKER_FLAGS_COVERAGE "--coverage") -set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage") -# set a flag to inform code that we are in debug mode -set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG") -set (ALLOWED_BUILD_TYPES Debug Release MinSizeRel RelWithDebInfo Coverage) -set_property (CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS ${ALLOWED_BUILD_TYPES}) -set (DEFAULT_BUILD_TYPE "Release") -if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) - set (CMAKE_BUILD_TYPE "${DEFAULT_BUILD_TYPE}" CACHE - STRING "Choose the type of build." FORCE) -endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) -#+++++++++++++++++++++++++++++ -# as long as there still are modules using it: -# -enable_language (Fortran) -set (CMAKE_Fortran_FLAGS "-std=legacy -Wfunction-elimination") - - #++++++++++++++++++++++++++++ # General config and flags # diff --git a/cmake/corsikaDefines.cmake b/cmake/corsikaDefines.cmake new file mode 100644 index 0000000000000000000000000000000000000000..60cf7da9a930a6bb7b39e35678e8b7064c411f42 --- /dev/null +++ b/cmake/corsikaDefines.cmake @@ -0,0 +1,86 @@ +# +# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu +# +# See file AUTHORS for a list of contributors. +# +# This software is distributed under the terms of the GNU General Public +# Licence version 3 (GPL Version 3). See file LICENSE for a full version of +# the license. +# + + +#+++++++++++++++++++++++++++++ +# as long as there still are modules using it: +# +enable_language (Fortran) +set (CMAKE_Fortran_FLAGS "-std=legacy -Wfunction-elimination") + + +#+++++++++++++++++++++++++++++ +# Build types settings +# +# setup coverage build type +set (CMAKE_CXX_FLAGS_COVERAGE "-g --coverage") +set (CMAKE_EXE_LINKER_FLAGS_COVERAGE "--coverage") +set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage") +# set a flag to inform code that we are in debug mode +set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG") + + +#+++++++++++++++++++++++++++++ +# Build type selection +# +# Set the possible values of build type for cmake-gui and command line check +set (ALLOWED_BUILD_TYPES Debug Release MinSizeRel RelWithDebInfo Coverage) +set_property (CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS ${ALLOWED_BUILD_TYPES}) + +set (DEFAULT_BUILD_TYPE "Release") +if (EXISTS "${CMAKE_SOURCE_DIR}/.git") + set (DEFAULT_BUILD_TYPE "Debug") +endif () + +if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + message (STATUS "Setting build type to '${DEFAULT_BUILD_TYPE}' as no other was specified.") + set (CMAKE_BUILD_TYPE "${DEFAULT_BUILD_TYPE}" CACHE + STRING "Choose the type of build." FORCE) +else (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + # Ignore capitalization when build type is selected manually and check for valid setting + string (TOLOWER ${CMAKE_BUILD_TYPE} SELECTED_LOWER) + string (TOLOWER "${ALLOWED_BUILD_TYPES}" BUILD_TYPES_LOWER) + if (NOT SELECTED_LOWER IN_LIST BUILD_TYPES_LOWER) + message (FATAL_ERROR "Unknown build type: ${CMAKE_BUILD_TYPE} [allowed: ${ALLOWED_BUILD_TYPES}]") + endif () + message (STATUS "Build type is: ${CMAKE_BUILD_TYPE}") +endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + + + +# +# Floating point exception support - select implementation to use +# +include (CheckIncludeFileCXX) +CHECK_INCLUDE_FILE_CXX ("fenv.h" HAS_FEENABLEEXCEPT) +if (HAS_FEENABLEEXCEPT) # FLOATING_POINT_ENVIRONMENT + set (CORSIKA_HAS_FEENABLEEXCEPT 1) + set_property (DIRECTORY ${CMAKE_HOME_DIRECTORY} APPEND PROPERTY COMPILE_DEFINITIONS "CORSIKA_HAS_FEENABLEEXCEPT") +endif () + + +# +# General OS Detection +# +if (${CMAKE_SYSTEM_NAME} MATCHES "Windows") + set (CORSIKA_OS_WINDOWS TRUE) + set (CORSIKA_OS "Windows") +elseif (${CMAKE_SYSTEM_NAME} MATCHES "Linux") + set (CORSIKA_OS_LINUX TRUE) + set (CORSIKA_OS "Linux") + # check for RedHat/CentOS compiler from the software collections (scl) + string (FIND ${CMAKE_CXX_COMPILER} "/opt/rh/devtoolset-" index) + if (${index} EQUAL 0) + set (CORSIKA_SCL_CXX TRUE) + endif () +elseif (${CMAKE_SYSTEM_NAME} MATCHES "Darwin") + set (CORSIKA_OS_MAC TRUE) + set (CORSIKA_OS "Mac") +endif () diff --git a/cmake/CorsikaUtilities.cmake b/cmake/corsikaUtilities.cmake similarity index 100% rename from cmake/CorsikaUtilities.cmake rename to cmake/corsikaUtilities.cmake diff --git a/conanfile.txt b/conanfile.txt index 85c1238b40b290ca181660b0a8bda665b1c88c5b..78cda273ada5bf3c944cb45f6f22df53d0567a8c 100644 --- a/conanfile.txt +++ b/conanfile.txt @@ -1,9 +1,28 @@ [requires] -spdlog/1.8.0 +readline/8.0 # this is only needed to fix a missing dependency in "bison" which is pulled-in by arrow +spdlog/1.8.5 catch2/2.13.3 eigen/3.3.8 -boost/1.74.0 +boost/1.76.0 zlib/1.2.11 +proposal/7.0.2 +yaml-cpp/0.6.3 +arrow/2.0.0 [generators] cmake + +[options] +readline:shared=True +arrow:shared=False +arrow:parquet=True +arrow:fPIC=False +arrow:with_re2=False +arrow:with_protobuf=False +arrow:with_openssl=False +arrow:with_gflags=False +arrow:with_glog=False +arrow:with_grpc=False +arrow:with_utf8proc=False +arrow:with_zstd=False +arrow:with_bz2=False diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl index 28306c590bb7da30d19d70f7ddfa19917ce525a0..654dfc6d1aad0b603baa76e701f5357fa85ffc13 100644 --- a/corsika/detail/framework/core/Cascade.inl +++ b/corsika/detail/framework/core/Cascade.inl @@ -26,11 +26,14 @@ namespace corsika { - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, typename TStackView> - inline void Cascade<TTracking, TProcessList, TStack, TStackView>::run() { + inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::run() { setNodes(); // put each particle on stack in correct environment volume + // start this event (i.e. this shower) + output_.startOfShower(); + while (!stack_.isEmpty()) { while (!stack_.isEmpty()) { CORSIKA_LOG_TRACE("Stack: {}", stack_.asString()); @@ -51,23 +54,27 @@ namespace corsika { // thus, the double loop // doCascadeEquations(); } + + // end this event (i.e. this shower) + output_.endOfShower(); } - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, typename TStackView> - inline void Cascade<TTracking, TProcessList, TStack, TStackView>::forceInteraction() { + inline void + Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::forceInteraction() { CORSIKA_LOG_TRACE("forced interaction!"); setNodes(); auto vParticle = stack_.getNextParticle(); TStackView secondaries(vParticle); - interaction(secondaries); + interaction(secondaries, sequence_.getInverseInteractionLength(vParticle)); sequence_.doSecondaries(secondaries); vParticle.erase(); // primary particle is done } - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, typename TStackView> - inline void Cascade<TTracking, TProcessList, TStack, TStackView>::step( + inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::step( Particle& vParticle) { // determine combined total interaction length (inverse) @@ -243,62 +250,89 @@ namespace corsika { [[maybe_unused]] auto projectile = secondaries.getProjectile(); if (distance_interact < distance_decay) { - interaction(secondaries); + interaction(secondaries, total_inv_lambda); } else { - decay(secondaries); - // make sure particle actually did decay if it should have done so - // \todo this should go to a validation code and not be included here - if (secondaries.getSize() == 1 && - projectile.getPID() == secondaries.getNextParticle().getPID()) - throw std::runtime_error(fmt::format("Particle {} decays into itself!", - get_name(projectile.getPID()))); + if (decay(secondaries, total_inv_lifetime) == ProcessReturn::Decayed) { + if (secondaries.getSize() == 1 && + projectile.getPID() == secondaries.getNextParticle().getPID()) { + throw std::runtime_error(fmt::format("Particle {} decays into itself!", + get_name(projectile.getPID()))); + } + } } sequence_.doSecondaries(secondaries); vParticle.erase(); } - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, typename TStackView> - inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay( - TStackView& view) { + inline ProcessReturn + Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::decay( + TStackView& view, InverseTimeType initial_inv_decay_time) { CORSIKA_LOG_DEBUG("decay"); + +#ifdef DEBUG InverseTimeType const actual_decay_time = sequence_.getInverseLifetime(view.parent()); + if (actual_decay_time * 0.99 > initial_inv_decay_time) { + CORSIKA_LOG_WARN( + "Decay time decreased during step! This leads to un-physical step length. " + "delta_inverse_decay_time={}", + 1 / initial_inv_decay_time - 1 / actual_decay_time); + } +#endif + + // one option is that decay_time is now larger (less + // probability for decay) than it was before the step, thus, + // no decay might actually occur and is allowed - UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time); + UniformRealDistribution<InverseTimeType> uniDist(initial_inv_decay_time); const auto sample_process = uniDist(rng_); auto const returnCode = sequence_.selectDecay(view, sample_process); if (returnCode != ProcessReturn::Decayed) { - CORSIKA_LOG_WARN("Particle did not decay!"); + CORSIKA_LOG_DEBUG("Particle did not decay!"); } setEventType(view, history::EventType::Decay); return returnCode; } - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, typename TStackView> - inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction( - TStackView& view) { + inline ProcessReturn + Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::interaction( + TStackView& view, InverseGrammageType initial_inv_int_length) { CORSIKA_LOG_DEBUG("collide"); - InverseGrammageType const current_inv_length = - sequence_.getInverseInteractionLength(view.parent()); +#ifdef DEBUG + InverseGrammageType const actual_inv_length = sequence_.getInverseInteractionLength( + view.parent()); // 1/lambda_int after step, -dE/dX etc. + + if (actual_inv_length * 0.99 > initial_inv_int_length) { + CORSIKA_LOG_WARN( + "Interaction length decreased during step! This leads to un-physical step " + "length. delta_innverse_interaction_length={}", + 1 / initial_inv_int_length - 1 / actual_inv_length); + } +#endif - UniformRealDistribution<InverseGrammageType> uniDist(current_inv_length); + // one option is that interaction_length is now larger (less + // probability for collision) than it was before the step, thus, + // no interaction might actually occur and is allowed + UniformRealDistribution<InverseGrammageType> uniDist(initial_inv_int_length); const auto sample_process = uniDist(rng_); auto const returnCode = sequence_.selectInteraction(view, sample_process); if (returnCode != ProcessReturn::Interacted) { - CORSIKA_LOG_WARN("Particle did not interace!"); + CORSIKA_LOG_DEBUG("Particle did not interact!"); } setEventType(view, history::EventType::Interaction); return returnCode; } - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, typename TStackView> - inline void Cascade<TTracking, TProcessList, TStack, TStackView>::setNodes() { + inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::setNodes() { std::for_each(stack_.begin(), stack_.end(), [&](auto& p) { auto const* numericalNode = environment_.getUniverse()->getContainingNode(p.getPosition()); @@ -306,9 +340,9 @@ namespace corsika { }); } - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, typename TStackView> - inline void Cascade<TTracking, TProcessList, TStack, TStackView>::setEventType( + inline void Cascade<TTracking, TProcessList, TOutput, TStack, TStackView>::setEventType( TStackView& view, [[maybe_unused]] history::EventType eventType) { if constexpr (TStackView::has_event) { for (auto&& sec : view) { sec.getEvent()->setEventType(eventType); } diff --git a/corsika/detail/framework/geometry/Path.inl b/corsika/detail/framework/geometry/Path.inl new file mode 100644 index 0000000000000000000000000000000000000000..c367531883c0bce37a488fb2901c6cc13c0bad09 --- /dev/null +++ b/corsika/detail/framework/geometry/Path.inl @@ -0,0 +1,66 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <deque> +#include <corsika/framework/geometry/Point.hpp> + +namespace corsika { + + Path::Path(Point const& point) { points_.push_front(point); } + + Path::Path(std::deque<Point> const& points) + : points_(points) { + int dequesize_ = points.size(); + if (dequesize_ == 0 || dequesize_ == 1) { + length_ = LengthType::zero(); + } else if (dequesize_ == 2) { + length_ = (points.back() - points.front()).getNorm(); + } else { + for (auto point = points.begin(); point != points.end() - 1; ++point) { + auto point_next = *(point + 1); + auto point_now = *(point); + length_ += (point_next - point_now).getNorm(); + } + } + } + + inline void Path::addToEnd(Point const& point) { + length_ += (point - points_.back()).getNorm(); + points_.push_back(point); + } + + inline void Path::removeFromEnd() { + auto lastpoint_ = points_.back(); + points_.pop_back(); + int dequesize_ = points_.size(); + if (dequesize_ == 0 || dequesize_ == 1) { + length_ = LengthType::zero(); + } else if (dequesize_ == 2) { + length_ = (points_.back() - points_.front()).getNorm(); + } else { + length_ -= (lastpoint_ - points_.back()).getNorm(); + } + } + + inline LengthType Path::getLength() const { return length_; } + + inline Point Path::getStart() const { return points_.front(); } + + inline Point Path::getEnd() const { return points_.back(); } + + inline Point Path::getPoint(std::size_t const index) const { return points_.at(index); } + + inline auto Path::begin() { return points_.begin(); } + + inline auto Path::end() { return points_.end(); } + + inline int Path::getNSegments() const { return points_.size() - 1; } + +} // namespace corsika \ No newline at end of file diff --git a/corsika/detail/framework/geometry/Point.inl b/corsika/detail/framework/geometry/Point.inl index c5c0a98a64491d3c6c752b647d60cdb1c2090e73..fe6e8b162383352523f7ba2399ef8685b0b5a2cb 100644 --- a/corsika/detail/framework/geometry/Point.inl +++ b/corsika/detail/framework/geometry/Point.inl @@ -102,4 +102,8 @@ namespace corsika { return os; } -} // namespace corsika + inline LengthType distance(Point const& p1, Point const& p2) { + return (p1 - p2).getNorm(); + } + +} // namespace corsika \ No newline at end of file diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl index 6c9bd789289a6fa24204bfe10df34e5413b5b6fa..5eb3d20ae785b7645101b362cb4697282a29ecdd 100644 --- a/corsika/detail/framework/process/ProcessSequence.inl +++ b/corsika/detail/framework/process/ProcessSequence.inl @@ -221,7 +221,7 @@ namespace corsika { IndexProcess2>::doStack(TStack& stack) { if constexpr (is_process_v<process1_type>) { // to protect from further compiler // errors if process1_type is invalid - if constexpr (std::is_base_of_v<StackProcess<process1_type>, process1_type> || + if constexpr (is_stack_process_v<process1_type> || (process1_type::is_process_sequence && !process1_type::is_switch_process_sequence)) { @@ -237,7 +237,7 @@ namespace corsika { } if constexpr (is_process_v<process2_type>) { // to protect from further compiler // errors if process2_type is invalid - if constexpr (std::is_base_of_v<StackProcess<process2_type>, process2_type> || + if constexpr (is_stack_process_v<process2_type> || (process2_type::is_process_sequence && !process2_type::is_switch_process_sequence)) { @@ -316,14 +316,14 @@ namespace corsika { if constexpr (is_process_v<process1_type>) { // to protect from further compiler // errors if process1_type is invalid - if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> || + if constexpr (is_interaction_process_v<process1_type> || process1_type::is_process_sequence) { tot += A_.getInverseInteractionLength(particle); } } if constexpr (is_process_v<process2_type>) { // to protect from further compiler // errors if process2_type is invalid - if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> || + if constexpr (is_interaction_process_v<process2_type> || process2_type::is_process_sequence) { tot += B_.getInverseInteractionLength(particle); } @@ -350,9 +350,7 @@ namespace corsika { A_.selectInteraction(view, lambda_inv_select, lambda_inv_sum); // if A_ did succeed, stop routine. Not checking other static branch B_. if (ret != ProcessReturn::Ok) { return ret; } - } else if constexpr (is_process_v<process1_type> && - std::is_base_of_v<InteractionProcess<process1_type>, - process1_type>) { + } else if constexpr (is_interaction_process_v<process1_type>) { // if this is not a ContinuousProcess --> evaluate probability lambda_inv_sum += A_.getInverseInteractionLength(view.parent()); // check if we should execute THIS process and then EXIT @@ -376,9 +374,7 @@ namespace corsika { if constexpr (process2_type::is_process_sequence) { // if B_ is a process sequence --> check inside return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum); - } else if constexpr (is_process_v<process2_type> && - std::is_base_of_v<InteractionProcess<process2_type>, - process2_type>) { + } else if constexpr (is_interaction_process_v<process2_type>) { // if this is not a ContinuousProcess --> evaluate probability lambda_inv_sum += B_.getInverseInteractionLength(view.parent()); // soon as SecondaryView::parent() is migrated! diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl index 2f87ade6d329ab249f5a61f21b79c72d07d315fd..634395f55e0a59524263df2f9d72ebe80cac8be6 100644 --- a/corsika/detail/framework/process/SwitchProcessSequence.inl +++ b/corsika/detail/framework/process/SwitchProcessSequence.inl @@ -36,13 +36,11 @@ namespace corsika { typename TParticle::node_type const& from, typename TParticle::node_type const& to) { if (select_(particle)) { - if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>, - process1_type> || + if constexpr (is_boundary_process_v<process1_type> || process1_type::is_process_sequence) { // interface checking on TSequence - if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>, - process1_type>) { + if constexpr (is_boundary_process_v<process1_type>) { static_assert( has_method_doBoundaryCrossing_v<TSequence, ProcessReturn, TParticle&>, @@ -56,13 +54,11 @@ namespace corsika { } } else { - if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>, - process2_type> || + if constexpr (is_boundary_process_v<process2_type> || process2_type::is_process_sequence) { // interface checking on USequence - if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>, - process2_type>) { + if constexpr (is_boundary_process_v<process2_type>) { static_assert( has_method_doBoundaryCrossing_v<USequence, ProcessReturn, TParticle>, @@ -136,7 +132,7 @@ namespace corsika { IndexProcess2>::doSecondaries(TSecondaries& vS) { const auto& particle = vS.parent(); if (select_(particle)) { - if constexpr (std::is_base_of_v<SecondariesProcess<process1_type>, process1_type> || + if constexpr (is_secondaries_process_v<process1_type> || process1_type::is_process_sequence) { // interface checking on TSequence @@ -150,7 +146,7 @@ namespace corsika { A_.doSecondaries(vS); } } else { - if constexpr (std::is_base_of_v<SecondariesProcess<process2_type>, process2_type> || + if constexpr (is_secondaries_process_v<process2_type> || process2_type::is_process_sequence) { // interface checking on USequence @@ -219,14 +215,14 @@ namespace corsika { IndexProcess2>::getInverseInteractionLength(TParticle&& particle) { if (select_(particle)) { - if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> || + if constexpr (is_interaction_process_v<process1_type> || process1_type::is_process_sequence) { return A_.getInverseInteractionLength(particle); } } else { - if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> || + if constexpr (is_interaction_process_v<process2_type> || process2_type::is_process_sequence) { return B_.getInverseInteractionLength(particle); } @@ -249,8 +245,7 @@ namespace corsika { A_.selectInteraction(view, lambda_inv_select, lambda_inv_sum); // if A_ did succeed, stop routine. Not checking other static branch B_. if (ret != ProcessReturn::Ok) { return ret; } - } else if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, - process1_type>) { + } else if constexpr (is_interaction_process_v<process1_type>) { // if this is not a ContinuousProcess --> evaluate probability lambda_inv_sum += A_.getInverseInteractionLength(view.parent()); // check if we should execute THIS process and then EXIT @@ -272,8 +267,7 @@ namespace corsika { if constexpr (process2_type::is_process_sequence) { // if B_ is a process sequence --> check inside return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum); - } else if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, - process2_type>) { + } else if constexpr (is_interaction_process_v<process2_type>) { // if this is not a ContinuousProcess --> evaluate probability lambda_inv_sum += B_.getInverseInteractionLength(view.parent()); // check if we should execute THIS process and then EXIT @@ -301,14 +295,14 @@ namespace corsika { IndexProcess2>::getInverseLifetime(TParticle&& particle) { if (select_(particle)) { - if constexpr (std::is_base_of_v<DecayProcess<process1_type>, process1_type> || + if constexpr (is_decay_process_v<process1_type> || process1_type::is_process_sequence) { return A_.getInverseLifetime(particle); } } else { - if constexpr (std::is_base_of_v<DecayProcess<process2_type>, process2_type> || + if constexpr (is_decay_process_v<process2_type> || process2_type::is_process_sequence) { return B_.getInverseLifetime(particle); } @@ -331,8 +325,7 @@ namespace corsika { ProcessReturn const ret = A_.selectDecay(view, decay_inv_select, decay_inv_sum); // if A_ did succeed, stop routine here (not checking other static branch B_) if (ret != ProcessReturn::Ok) { return ret; } - } else if constexpr (std::is_base_of_v<DecayProcess<process1_type>, - process1_type>) { + } else if constexpr (is_decay_process_v<process1_type>) { // if this is not a ContinuousProcess --> evaluate probability decay_inv_sum += A_.getInverseLifetime(view.parent()); // check if we should execute THIS process and then EXIT @@ -355,8 +348,7 @@ namespace corsika { if constexpr (process2_type::is_process_sequence) { // if B_ is a process sequence --> check inside return B_.selectDecay(view, decay_inv_select, decay_inv_sum); - } else if constexpr (std::is_base_of_v<DecayProcess<process2_type>, - process2_type>) { + } else if constexpr (is_decay_process_v<process2_type>) { // if this is not a ContinuousProcess --> evaluate probability decay_inv_sum += B_.getInverseLifetime(view.parent()); // check if we should execute THIS process and then EXIT diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl new file mode 100644 index 0000000000000000000000000000000000000000..c15594e02a02bef2683e4d3ae4d4ad02895ef747 --- /dev/null +++ b/corsika/detail/framework/utility/CubicSolver.inl @@ -0,0 +1,294 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/utility/CubicSolver.hpp> +#include <cmath> + +namespace corsika { + + namespace andre { + + //--------------------------------------------------------------------------- + // solve cubic equation A x^3 + B*x^2 + C*x + D = 0 + inline std::vector<double> solve_cubic_real_analytic(long double A, long double B, + long double C, long double D, + double const epsilon) { + + if (std::abs(A) < epsilon) { return solve_quadratic_real(B, C, epsilon); } + + long double a = B / A; + long double b = C / A; + long double c = D / A; + + long double a2 = a * a; + long double q = (a2 - 3 * b) / 9; + long double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54; + long double q3 = q * q * q; + + // disc = q**3 + r**2 + // w:e = r*r exactly + long double w = r * r; + long double e = std::fma(r, r, -w); + // s:t = q*q exactly + long double s = -q * q; + long double t = std::fma(-q, q, -s); + // s:t * q + w:e = s*q + w + t*q +e = s*q+w + u:v + e = f + u:v + e + long double f = std::fma(s, q, w); + long double u = t * q; + long double v = std::fma(t, q, -u); + // error-free sum f+u. See Knuth, TAOCP vol. 2 + long double uf = u + f; + long double au = uf - u; + long double ab = uf - au; + au = f - au; + ab = u - ab; + // sum all terms into final result + long double const disc = -(((e + uf) + au) + ab) + v; + + if (disc >= 0) { + long double t = r / std::sqrt(q3); + if (t < -1) t = -1; + if (t > 1) t = 1; + t = std::acos(t); + a /= 3; + q = -2 * std::sqrt(q); + return {double(q * std::cos(t / 3) - a), + double(q * std::cos((t + 2 * M_PI) / 3) - a), + double(q * std::cos((t - 2 * M_PI) / 3) - a)}; + } else { + long double term1 = -cbrt(std::fabs(r) + std::sqrt(-disc)); + if (r < 0) term1 = -term1; + long double term2 = (0 == term1 ? 0 : q / term1); + + a /= 3; + long double test = 0.5 * std::sqrt(3.) * (term1 - term2); + if (std::fabs(test) < epsilon) { + return {double((term1 + term2) - 1), double(-0.5 * (term1 + term2) - a)}; + } + + return {double((term1 + term2) - a)}; + } + } + } // namespace andre + + inline std::vector<double> solve_cubic_depressed_disciminant_real( + long double p, long double q, long double const disc, double const epsilon) { + + CORSIKA_LOG_TRACE("p={}, q={}, disc={}", p, q, disc); + + if (std::abs(disc) < epsilon) { // disc==0 multiple roots ! + if (std::abs(p) < epsilon) { // tripple root + return {0}; + } + // double root, single root + CORSIKA_LOG_TRACE("cubic double root"); + return {double(-3 * q / (2 * p)), double(3 * q / p)}; + } + + if (std::abs(p) < epsilon) { // p==0 --> x^3 + q = 0 + return {double(std::cbrt(-q))}; + } + + if (disc > 0) { // three real roots + CORSIKA_LOG_TRACE("cubic three roots"); + long double const p_third = p / 3; + long double const sqrt_minus_p_third = std::sqrt(-p_third); + + long double const arg = std::acos(q / (2 * p_third) / sqrt_minus_p_third) / 3; + + long double constexpr pi = M_PI; + return {double(2 * sqrt_minus_p_third * std::cos(arg)), + double(2 * sqrt_minus_p_third * std::cos(arg - 2 * pi / 3)), + double(2 * sqrt_minus_p_third * std::cos(arg - 4 * pi / 3))}; + } + + // thus disc < 0 -> one real root + if (p < 0) { + CORSIKA_LOG_TRACE("cubic cosh"); + long double const abs_q = std::abs(q); + long double const p_third = p / 3; + long double const sqrt_minus_p_third = std::sqrt(-p_third); + CORSIKA_LOG_TRACE("sqrt_minus_p_third={}, arcosh({})={}", sqrt_minus_p_third, + -abs_q / (2 * p_third) / sqrt_minus_p_third, + std::acosh(-abs_q / (2 * p_third) / sqrt_minus_p_third)); + CORSIKA_LOG_TRACE( + "complex: {}", + -2 * abs_q / q * sqrt_minus_p_third * + std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3)); + double const z = + -2 * abs_q / q * sqrt_minus_p_third * + std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3); + return {z}; + } else { // p > 0 + CORSIKA_LOG_TRACE("cubic sinh"); + long double const p_third = p / 3; + long double const sqrt_p_third = std::sqrt(p_third); + return {double(-2 * sqrt_p_third * + std::sinh(std::asinh(q / (2 * p_third * sqrt_p_third)) / 3))}; + } + } + + inline std::vector<double> solve_cubic_depressed_real(long double p, long double q, + double const epsilon) { + + // thanks! + // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic + // long double const disc = -(4 * p * p * p + 27 * q * q); + // disc = (p/3)**3 + (q/2)**2 + long double p_third = p / 3; + long double q_half = q / 2; + // w:e = (q/2)*(q/2) exactly + long double w = q_half * q_half; + long double e = std::fma(q_half, q_half, -w); + // s:t = (p/3)*(p/3) exactly + long double s = p_third * p_third; + long double t = std::fma(p_third, p_third, -s); + // s:t * p + w:e = s*p + w + t*p +e = s*p+w + u:v + e = f + u:v + e + long double f = std::fma(s, p_third, w); + long double u = t * p_third; + long double v = std::fma(t, p_third, -u); + // error-free sum f+u. See Knuth, TAOCP vol. 2 + long double a = u + f; + long double b = a - u; + long double c = a - b; + b = f - b; + c = u - c; + // sum all terms into final result + long double const disc = -(((e + a) + b) + c) + v; + return solve_cubic_depressed_disciminant_real(p, q, disc, epsilon); + } + + /** + * Analytical approach. Not very stable in all conditions. + */ + inline std::vector<double> solve_cubic_real_analytic(long double a, long double b, + long double c, long double d, + double const epsilon) { + + CORSIKA_LOG_TRACE("cubic: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b, c, + d, epsilon, (std::abs(a - 1) < epsilon), (std::abs(b) < epsilon)); + + if (std::abs(a) < epsilon) { // this is just a quadratic + return solve_quadratic_real(b, c, d, epsilon); + } + + if ((std::abs(a - 1) < epsilon) && + (std::abs(b) < epsilon)) { // this is a depressed cubic + return solve_cubic_depressed_real(c, d, epsilon); + } + + // p = (3ac - b^2) / (3a^2) = 3 * ( c/(3*a) - b^2/(9*a^2) ) + long double b_over_a = b / a; + long double const p_third = std::fma(-b_over_a, b_over_a / 9, c / (a * 3)); + + // p = std::fma(a * 3, c, -b2) / (3 * a2); + // q = (2*b^3 - 9*abc + 27*a^2*d) / (27a^3) = 2 * ( b^3/(27*a^3) - bc/(6*a^2) + + // d/(2*a) ) + long double const q_half_term1 = std::fma(b_over_a, b_over_a / 27, -c / (a * 6)); + long double const q_half = std::fma(b_over_a, q_half_term1, d / (a * 2)); + + std::vector<double> zs = solve_cubic_depressed_real(3 * p_third, 2 * q_half, epsilon); + CORSIKA_LOG_TRACE("cubic: solve_depressed={}, b/3a={}", fmt::join(zs, ", "), + b / (a * 3)); + for (auto& z : zs) { + z -= b / (a * 3); + double const b1 = z + b / a; + double const b0 = b1 * z + c / a; + std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon); + CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "), + static_pow<3>(z) * a + static_pow<2>(z) * b + z * c + d); + } + CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(zs, ", ")); + return zs; + } + + template <typename T> // T must be floating point type + inline T cubic_function(T x, T a, T b, T c, T d) { + T x2 = x * x; + return x2 * x * a + x2 * b + x * c + d; + } + + template <typename T> // T must be floating point type + inline T cubic_function_dfdx(T x, T a, T b, T c) { + T x2 = x * x; + return x2 * a * 3 + x * b * 2 + c; + } + + template <typename T> // T must be floating point type + inline T cubic_function_d2fd2x(T x, T a, T b) { + return x * a * 6 + b * 2; + } + + /** + * Iterative approach. + */ + + inline std::vector<double> solve_cubic_real(long double a, long double b, long double c, + long double d, double const epsilon) { + + CORSIKA_LOG_TRACE("cubic_2: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b, + c, d, epsilon, (std::abs(a - 1) < epsilon), + (std::abs(b) < epsilon)); + + if (std::abs(a) < epsilon) { // this is just a quadratic + return solve_quadratic_real(b, c, d, epsilon); + } + + long double const dist = std::fma(b / a, b / a, -3 * c / a); + long double const xinfl = -b / (a * 3); + + long double x1 = xinfl; + long double f_x1 = cubic_function(xinfl, a, b, c, d); + + if (std::abs(f_x1) > epsilon) { + if (std::abs(dist) < epsilon) { + x1 = xinfl - std::cbrt(f_x1); + } else if (dist > 0) { + if (f_x1 > 0) + x1 = xinfl - 2 / 3 * std::sqrt(dist); + else + x1 = xinfl + 2 / 3 * std::sqrt(dist); + } + + int niter = 0; + const int maxiter = 100; + do { + long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c); + long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b); + // if (potential) saddle point... avoid + if (std::abs(f_prime_x1) < epsilon) { + x1 -= std::cbrt(f_x1); + } else { + x1 -= + f_x1 * f_prime_x1 / (static_pow<2>(f_prime_x1) - 0.5 * f_x1 * f_prime2_x1); + } + f_x1 = cubic_function(x1, a, b, c, d); + CORSIKA_LOG_TRACE("niter={} x1={} f_x1={} f_prime={} f_prime2={} eps={}", niter, + x1, f_x1, f_prime_x1, f_prime2_x1, epsilon); + } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon)); + + CORSIKA_LOG_TRACE("niter={}", niter); + } + + CORSIKA_LOG_TRACE("x1={} f_x1={}", x1, f_x1); + + double const b1 = x1 + b / a; + double const b0 = b1 * x1 + c / a; + std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon); + CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "), + cubic_function(x1, a, b, c, d)); + + quad_check.push_back(x1); + CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(quad_check, ", ")); + return quad_check; + } + +} // namespace corsika diff --git a/corsika/detail/framework/utility/LinearSolver.inl b/corsika/detail/framework/utility/LinearSolver.inl new file mode 100644 index 0000000000000000000000000000000000000000..6eccee91ce4614d5cf95cd0948c6497824741213 --- /dev/null +++ b/corsika/detail/framework/utility/LinearSolver.inl @@ -0,0 +1,35 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <vector> +#include <cmath> +#include <complex> + +namespace corsika { + + inline std::vector<double> solve_linear_real(double a, double b) { + + if (a == 0) { + return {}; // no (b!=0), or infinite number (b==0) of solutions.... + } + + return {-b / a}; + } + + inline std::vector<std::complex<double>> solve_linear(double a, double b) { + + if (std::abs(a) == 0) { + return {}; // no (b!=0), or infinite number (b==0) of solutions.... + } + + return {std::complex<double>(-b / a, 0)}; + } + +} // namespace corsika diff --git a/corsika/detail/framework/utility/QuadraticSolver.inl b/corsika/detail/framework/utility/QuadraticSolver.inl new file mode 100644 index 0000000000000000000000000000000000000000..9c9759e975ae6d322dee7a478f935c095bdf0ff9 --- /dev/null +++ b/corsika/detail/framework/utility/QuadraticSolver.inl @@ -0,0 +1,83 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <corsika/framework/core/PhysicalUnits.hpp> + +namespace corsika { + + inline std::vector<std::complex<double>> solve_quadratic(long double a, long double b, + long double c, + double const epsilon) { + if (std::abs(a) < epsilon) { return solve_linear(b, c); } + + if (std::abs(c) < epsilon) { + std::vector<std::complex<double>> lin_result = solve_linear(a, b); + lin_result.push_back({0.}); + return lin_result; + } + + long double const radicant = static_pow<2>(b) - a * c * 4; + + if (radicant < -epsilon) { // two complex solutions + double const rpart = -b / 2 * a; + double const ipart = std::sqrt(-radicant); + return {{rpart, ipart}, {rpart, -ipart}}; + } + + if (radicant < epsilon) { // just one real solution + return {{double(-b / 2 * a), 0}}; + } + + // two real solutions, use Citardauq formula and actively avoid cancellation in the + // denominator + + const long double x1 = + 2 * c / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant)); + + return {{double(x1), 0}, {double(c / (a * x1)), 0}}; + } + + inline std::vector<double> solve_quadratic_real(long double a, long double b, + long double c, double const epsilon) { + + CORSIKA_LOG_TRACE("quadratic: a={} b={} c={}", a, b, c); + + if (std::abs(a) < epsilon) { return solve_linear_real(b, c); } + if (std::abs(c) < epsilon) { + std::vector<double> lin_result = solve_linear_real(a, b); + lin_result.push_back(0.); + return lin_result; + } + + // long double const radicant = std::pow(b, 2) - a * c * 4; + // Thanks! + // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic + long double w = a * 4 * c; + long double e = std::fma(a * 4, c, -w); + long double f = std::fma(b, b, -w); + long double radicant = f + e; + + CORSIKA_LOG_TRACE(" radicant={} {} ", radicant, b * b - a * c * 4); + + if (std::abs(radicant) < epsilon) { // just one real solution + return {double(-b / (2 * a))}; + } + + if (radicant < 0) { return {}; } + + // two real solutions, use Citardauq formula and actively avoid cancellation in the + // denominator + + long double const x1 = + c * 2 / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant)); + + CORSIKA_LOG_TRACE("quadratic x1={} x2={}", double(x1), double(c / (a * x1))); + + return {double(x1), double(c / (a * x1))}; + } +} // namespace corsika diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl index a82abf7ac022c3c99f6e342f88eb35f2ddaa6425..988bbe80b315d59344c8c71f8960ac5833e270aa 100644 --- a/corsika/detail/framework/utility/QuarticSolver.inl +++ b/corsika/detail/framework/utility/QuarticSolver.inl @@ -8,128 +8,154 @@ #pragma once +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/utility/CubicSolver.hpp> #include <cmath> -namespace corsika::quartic_solver { - - //--------------------------------------------------------------------------- - // solve cubic equation x^3 + a*x^2 + b*x + c = 0 - // x - array of size 3 - // In case 3 real roots: => x[0], x[1], x[2], return 3 - // 2 real roots: x[0], x[1], return 2 - // 1 real root : x[0], x[1] ± i*x[2], return 1 - inline unsigned int solveP3(double* x, double a, double b, double c) { - double a2 = a * a; - double q = (a2 - 3 * b) / 9; - double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54; - double r2 = r * r; - double q3 = q * q * q; - double A, B; - if (r2 < q3) { - double t = r / sqrt(q3); - if (t < -1) t = -1; - if (t > 1) t = 1; - t = acos(t); - a /= 3; - q = -2 * sqrt(q); - x[0] = q * cos(t / 3) - a; - x[1] = q * cos((t + M_2PI) / 3) - a; - x[2] = q * cos((t - M_2PI) / 3) - a; - return 3; - } else { - A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3); - if (r < 0) A = -A; - B = (0 == A ? 0 : q / A); - - a /= 3; - x[0] = (A + B) - a; - x[1] = -0.5 * (A + B) - a; - x[2] = 0.5 * sqrt(3.) * (A - B); - if (fabs(x[2]) < eps) { - x[2] = x[1]; - return 2; +namespace corsika { + + namespace andre { + + inline std::vector<double> solve_quartic_real(long double a, long double b, + long double c, long double d, + long double e, double const epsilon) { + + if (std::abs(a) < epsilon) { return solve_cubic_real(b, c, d, e, epsilon); } + + b /= a; + c /= a; + d /= a; + e /= a; + + long double a3 = -c; + long double b3 = b * d - 4. * e; + long double c3 = -b * b * e - d * d + 4. * c * e; + + // cubic resolvent + // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0 + + std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon); + long double y = x3[0]; // there is always at least one solution + // The essence - choosing Y with maximal absolute value. + if (x3.size() == 3) { + if (fabs(x3[1]) > fabs(y)) y = x3[1]; + if (fabs(x3[2]) > fabs(y)) y = x3[2]; + } + + long double q1, q2, p1, p2; + // h1+h2 = y && h1*h2 = e <=> h^2 -y*h + e = 0 (h === q) + + long double Det = y * y - 4 * e; + CORSIKA_LOG_TRACE("Det={}", Det); + if (fabs(Det) < epsilon) // in other words - D==0 + { + q1 = q2 = y * 0.5; + // g1+g2 = b && g1+g2 = c-y <=> g^2 - b*g + c-y = 0 (p === g) + Det = b * b - 4 * (c - y); + if (fabs(Det) < epsilon) { // in other words - D==0 + p1 = p2 = b * 0.5; + } else { + if (Det < 0) return {}; + long double sqDet = sqrt(Det); + p1 = (b + sqDet) * 0.5; + p2 = (b - sqDet) * 0.5; + } + } else { + if (Det < 0) return {}; + long double sqDet1 = sqrt(Det); + q1 = (y + sqDet1) * 0.5; + q2 = (y - sqDet1) * 0.5; + // g1+g2 = b && g1*h2 + g2*h1 = c ( && g === p ) Krammer + p1 = (b * q1 - d) / (q1 - q2); + p2 = (d - b * q2) / (q1 - q2); } - return 1; + // solving quadratic eqs. x^2 + p1*x + q1 = 0 + // x^2 + p2*x + q2 = 0 + + std::vector<double> quad1 = solve_quadratic_real(1, p1, q1); + std::vector<double> quad2 = solve_quadratic_real(1, p2, q2); + if (quad2.size() > 0) { + for (auto val : quad2) quad1.push_back(val); + } + return quad1; } - } + } // namespace andre - //--------------------------------------------------------------------------- - // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d - // Attention - this function returns dynamically allocated array. It has to be released - // afterwards. - inline DComplex* solve_quartic(double a, double b, double c, double d) { - double a3 = -b; - double b3 = a * c - 4. * d; - double c3 = -a * a * d - c * c + 4. * b * d; + inline std::vector<double> solve_quartic_depressed_real(long double p, long double q, + long double r, + double const epsilon) { - // cubic resolvent - // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0 + CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f}, epsilon={}", p, q, r, + epsilon); - double x3[3]; - unsigned int iZeroes = solveP3(x3, a3, b3, c3); + long double const p2 = static_pow<2>(p); + long double const q2 = static_pow<2>(q); - double q1, q2, p1, p2, D, sqD, y; + std::vector<double> const resolve_cubic = + solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon); - y = x3[0]; - // The essence - choosing Y with maximal absolute value. - if (iZeroes != 1) { - if (fabs(x3[1]) > fabs(y)) y = x3[1]; - if (fabs(x3[2]) > fabs(y)) y = x3[2]; - } + CORSIKA_LOG_TRACE("resolve_cubic: N={}, m=[{}]", resolve_cubic.size(), + fmt::join(resolve_cubic, ", ")); - // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q) - - D = y * y - 4 * d; - if (fabs(D) < eps) // in other words - D==0 - { - q1 = q2 = y * 0.5; - // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g) - D = a * a - 4 * (b - y); - if (fabs(D) < eps) // in other words - D==0 - p1 = p2 = a * 0.5; - - else { - sqD = sqrt(D); - p1 = (a + sqD) * 0.5; - p2 = (a - sqD) * 0.5; - } - } else { - sqD = sqrt(D); - q1 = (y + sqD) * 0.5; - q2 = (y - sqD) * 0.5; - // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer - p1 = (a * q1 - c) / (q1 - q2); - p2 = (c - a * q2) / (q1 - q2); + if (!resolve_cubic.size()) return {}; + + long double m = 0; + for (auto const& v : resolve_cubic) { + CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p + + v * (p2 / 4 - r) - q2 / 8)); + if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; } } + CORSIKA_LOG_TRACE("check m={}", m); + if (m == 0) { return {0}; } + + CORSIKA_LOG_TRACE("check m={}", m); + + long double const quad_term1 = p / 2 + m; + long double const quad_term2 = std::sqrt(2 * m); + long double const quad_term3 = q / (2 * quad_term2); - DComplex* retval = new DComplex[4]; - - // solving quadratic eq. - x^2 + p1*x + q1 = 0 - D = p1 * p1 - 4 * q1; - if (D < 0.0) { - retval[0].real(-p1 * 0.5); - retval[0].imag(sqrt(-D) * 0.5); - retval[1] = std::conj(retval[0]); - } else { - sqD = sqrt(D); - retval[0].real((-p1 + sqD) * 0.5); - retval[1].real((-p1 - sqD) * 0.5); + std::vector<double> z_quad1 = + solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, epsilon); + std::vector<double> z_quad2 = + solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, epsilon); + for (auto const& z : z_quad2) z_quad1.push_back(z); + return z_quad1; + } + + inline std::vector<double> solve_quartic_real(long double a, long double b, + long double c, long double d, + long double e, double const epsilon) { + + CORSIKA_LOG_TRACE("quartic: a={:f}, b={:f}, c={:f}, d={:f}, e={:f}, epsilon={}", a, b, + c, d, e, epsilon); + + if (std::abs(a) < epsilon) { // this is just a quadratic + return solve_cubic_real(b, c, d, e, epsilon); } - // solving quadratic eq. - x^2 + p2*x + q2 = 0 - D = p2 * p2 - 4 * q2; - if (D < 0.0) { - retval[2].real(-p2 * 0.5); - retval[2].imag(sqrt(-D) * 0.5); - retval[3] = std::conj(retval[2]); - } else { - sqD = sqrt(D); - retval[2].real((-p2 + sqD) * 0.5); - retval[3].real((-p2 - sqD) * 0.5); + if ((std::abs(a - 1) < epsilon) && + (std::abs(b) < epsilon)) { // this is a depressed quartic + return solve_quartic_depressed_real(c, d, e, epsilon); } - return retval; + long double const b2 = static_pow<2>(b); + long double const b3 = static_pow<3>(b); + long double const b4 = static_pow<4>(b); + long double const a2 = static_pow<2>(a); + long double const a3 = static_pow<3>(a); + long double const a4 = static_pow<4>(a); + + long double const p = (c * a * 8 - b2 * 3) / (a4 * 8); + long double const q = (b3 - b * c * a * 4 + d * a2 * 8) / (a4 * 8); + long double const r = + (-b4 * 3 + e * a3 * 256 - b * d * a2 * 64 + b2 * c * a * 16) / (a4 * 256); + + std::vector<double> zs = solve_quartic_depressed_real(p, q, r, epsilon); + CORSIKA_LOG_TRACE("quartic: solve_depressed={}, b/4a={}", fmt::join(zs, ", "), + b / (4 * a)); + for (auto& z : zs) { z -= b / (4 * a); } + CORSIKA_LOG_TRACE("quartic: solve_quartic_real returns={}", fmt::join(zs, ", ")); + return zs; } - -} // namespace corsika::quartic_solver \ No newline at end of file +} // namespace corsika diff --git a/corsika/detail/media/ExponentialRefractiveIndex.inl b/corsika/detail/media/ExponentialRefractiveIndex.inl new file mode 100644 index 0000000000000000000000000000000000000000..ac60050c50147ad8a58672ae3459bb1d1192eeb6 --- /dev/null +++ b/corsika/detail/media/ExponentialRefractiveIndex.inl @@ -0,0 +1,28 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/media/IRefractiveIndexModel.hpp> + +namespace corsika { + + template <typename T> + template <typename... Args> + ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex( + double const n0, InverseLengthType const lambda, Args&&... args) + : T(std::forward<Args>(args)...) + , n_0(n0) + , lambda_(lambda) {} + + template <typename T> + double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const { + return n_0 * exp((-lambda_) * point.getCoordinates().getZ()); + } + +} // namespace corsika diff --git a/corsika/detail/media/ShowerAxis.inl b/corsika/detail/media/ShowerAxis.inl index a42989e9222f0c4ec6131d42fe37110ff62faf5a..4cabe21b5a3e4cee53cfedcaf6b211f9b93a47fa 100644 --- a/corsika/detail/media/ShowerAxis.inl +++ b/corsika/detail/media/ShowerAxis.inl @@ -84,6 +84,7 @@ namespace corsika { l / max_length_); throw std::runtime_error(err.c_str()); } + return getMaximumX(); } CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}", fraction, fractionalBin, lower, upper); diff --git a/corsika/detail/modules/LongitudinalProfile.inl b/corsika/detail/modules/LongitudinalProfile.inl index 800d88240f06f5a02b683f934f731352eeef14fd..643ea017313a7358fdfd22a7e02643975bcebfec 100644 --- a/corsika/detail/modules/LongitudinalProfile.inl +++ b/corsika/detail/modules/LongitudinalProfile.inl @@ -35,15 +35,15 @@ namespace corsika { GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0)); GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1)); - CORSIKA_LOG_TRACE("pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2", + CORSIKA_LOG_DEBUG("longprof: pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2", vTrack.getPosition(0).getCoordinates() / 1_m, vTrack.getPosition(1).getCoordinates() / 1_m, grammageStart / 1_g * square(1_cm), grammageEnd / 1_g * square(1_cm)); // Note: particle may go also "upward", thus, grammageEnd<grammageStart - const int binStart = std::ceil(grammageStart / dX_); - const int binEnd = std::floor(grammageEnd / dX_); + int const binStart = std::ceil(grammageStart / dX_); + int const binEnd = std::floor(grammageEnd / dX_); for (int b = binStart; b <= binEnd; ++b) { if (pid == Code::Photon) { @@ -66,6 +66,7 @@ namespace corsika { inline void LongitudinalProfile::save(std::string const& filename, const int width, const int precision) { + CORSIKA_LOG_DEBUG("Write longprof to {}", filename); std::ofstream f{filename}; f << "# X / g·cm¯², photon, e+, e-, mu+, mu-, all hadrons" << std::endl; for (size_t b = 0; b < profiles_.size(); ++b) { diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl index c4da50e940f88e3b8a3d3e9f36ac0018eeb09466..12614573f13888396cc28e3655b9119fa7af2bec 100644 --- a/corsika/detail/modules/ObservationPlane.inl +++ b/corsika/detail/modules/ObservationPlane.inl @@ -6,44 +6,51 @@ * the license. */ -#include <corsika/modules/ObservationPlane.hpp> -#include <corsika/framework/core/Logging.hpp> +#include <corsika/setup/SetupStack.hpp> #include <corsika/setup/SetupTrajectory.hpp> -#include <fstream> - namespace corsika { - inline ObservationPlane::ObservationPlane(Plane const& obsPlane, - DirectionVector const& x_axis, - std::string const& filename, bool deleteOnHit) + template <typename TOutput> + ObservationPlane<TOutput>::ObservationPlane(Plane const& obsPlane, + DirectionVector const& x_axis, + bool deleteOnHit) : plane_(obsPlane) - , outputStream_(filename) , deleteOnHit_(deleteOnHit) , energy_ground_(0_GeV) , count_ground_(0) , xAxis_(x_axis.normalized()) - , yAxis_(obsPlane.getNormal().cross(xAxis_)) { - outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" - << std::endl; - } + , yAxis_(obsPlane.getNormal().cross(xAxis_)) {} - inline ProcessReturn ObservationPlane::doContinuous( + template <typename TOutput> + inline ProcessReturn ObservationPlane<TOutput>::doContinuous( corsika::setup::Stack::particle_type& particle, corsika::setup::Trajectory&, bool const stepLimit) { - /* The current step did not yet reach the ObservationPlane, do nothing now and wait: */ - if (!stepLimit) { return ProcessReturn::Ok; } + if (!stepLimit) { +#ifdef DEBUG + if (deleteOnHit_) { + LengthType const check = + (particle.getPosition() - plane_.getCenter()).dot(plane_.getNormal()); + if (check < 0_m) { + CORSIKA_LOG_DEBUG("PARTICLE AVOIDED OBSERVATIONPLANE {}", check); + } + } +#endif + return ProcessReturn::Ok; + } HEPEnergyType const energy = particle.getEnergy(); Point const pointOfIntersection = particle.getPosition(); Vector const displacement = pointOfIntersection - plane_.getCenter(); - outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV - << ' ' << displacement.dot(xAxis_) / 1_m << ' ' - << displacement.dot(yAxis_) / 1_m << '\n'; + // add our particles to the output file stream + this->write(particle.getPID(), energy, displacement.dot(xAxis_), + displacement.dot(yAxis_)); + + CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_); if (deleteOnHit_) { count_ground_++; @@ -54,17 +61,19 @@ namespace corsika { } } - inline LengthType ObservationPlane::getMaxStepLength( + template <typename TOutput> + inline LengthType ObservationPlane<TOutput>::getMaxStepLength( corsika::setup::Stack::particle_type const& particle, corsika::setup::Trajectory const& trajectory) { + CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}", particle.asString(), + particle.getPosition(), particle.getDirection(), plane_.asString()); + Intersections const intersection = setup::Tracking::intersect<corsika::setup::Stack::particle_type>(particle, plane_); TimeType const timeOfIntersection = intersection.getEntry(); - CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}, timeOfIntersection={}", - particle.asString(), particle.getPosition(), - particle.getDirection(), plane_.asString(), timeOfIntersection); + CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection); if (timeOfIntersection < TimeType::zero()) { return std::numeric_limits<double>::infinity() * 1_m; } @@ -78,17 +87,62 @@ namespace corsika { return dist; } - inline void ObservationPlane::showResults() const { + template <typename TOutput> + inline void ObservationPlane<TOutput>::showResults() const { CORSIKA_LOG_INFO( " ******************************\n" " ObservationPlane: \n" - " energy in ground (GeV) : {}\n" - " no. of particles in ground : {}\n" + " energy an ground (GeV) : {}\n" + " no. of particles at ground : {}\n" " ******************************", energy_ground_ / 1_GeV, count_ground_); } - inline void ObservationPlane::reset() { + template <typename TOutput> + inline YAML::Node ObservationPlane<TOutput>::getConfig() const { + using namespace units::si; + + // construct the top-level node + YAML::Node node; + + // basic info + node["type"] = "ObservationPlane"; + node["units"] = "m"; // add default units for values + + // the center of the plane + auto const center{plane_.getCenter()}; + + // save each component in its native coordinate system + auto const center_coords{center.getCoordinates(center.getCoordinateSystem())}; + node["plane"]["center"].push_back(center_coords.getX() / 1_m); + node["plane"]["center"].push_back(center_coords.getY() / 1_m); + node["plane"]["center"].push_back(center_coords.getZ() / 1_m); + + // the normal vector of the plane + auto const normal{plane_.getNormal().getComponents()}; + node["plane"]["normal"].push_back(normal.getX().magnitude()); + node["plane"]["normal"].push_back(normal.getY().magnitude()); + node["plane"]["normal"].push_back(normal.getZ().magnitude()); + + // the x-axis vector + auto const xAxis_coords{xAxis_.getComponents(xAxis_.getCoordinateSystem())}; + node["x-axis"].push_back(xAxis_coords.getX().magnitude()); + node["x-axis"].push_back(xAxis_coords.getY().magnitude()); + node["x-axis"].push_back(xAxis_coords.getZ().magnitude()); + + // the y-axis vector + auto const yAxis_coords{yAxis_.getComponents(yAxis_.getCoordinateSystem())}; + node["y-axis"].push_back(yAxis_coords.getX().magnitude()); + node["y-axis"].push_back(yAxis_coords.getY().magnitude()); + node["y-axis"].push_back(yAxis_coords.getZ().magnitude()); + + node["delete_on_hit"] = deleteOnHit_; + + return node; + } + + template <typename TOutput> + inline void ObservationPlane<TOutput>::reset() { energy_ground_ = 0_GeV; count_ground_ = 0; } diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl index caae6c59c8cb43b01cd1e9789366b51db2f9ec14..2d9c8238a4a5bb74b6ae80d6b218df0f1ef31221 100644 --- a/corsika/detail/modules/StackInspector.inl +++ b/corsika/detail/modules/StackInspector.inl @@ -72,13 +72,13 @@ namespace corsika { CORSIKA_LOG_INFO( "StackInspector: " - " time= {}" - ", running= {} seconds" + " time={}" + ", running={} seconds" " ( {}%)" - ", nStep= {}" - ", stackSize= {}" - ", Estack= {} GeV" - ", ETA=", + ", nStep={}" + ", stackSize={}" + ", Estack={} GeV" + ", ETA={}", std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(), int(progress * 100), getStep(), vS.getSize(), Etot / 1_GeV, std::put_time(std::localtime(&eta_time), "%T")); diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl index ebd52339e72eddde23f1f50dc66e709712165963..8af3c31dc45b718d8213346e3896e2dffb02812d 100644 --- a/corsika/detail/modules/TrackWriter.inl +++ b/corsika/detail/modules/TrackWriter.inl @@ -8,54 +8,48 @@ #pragma once -#include <corsika/modules/TrackWriter.hpp> - #include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> -#include <corsika/setup/SetupStack.hpp> -#include <corsika/setup/SetupTrajectory.hpp> - -#include <iomanip> #include <limits> namespace corsika { - inline TrackWriter::TrackWriter(std::string const& filename) - : filename_(filename) { - using namespace std::string_literals; - - file_.open(filename_); - file_ - << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s - << '\n'; - } + template <typename TOutputWriter> + inline TrackWriter<TOutputWriter>::TrackWriter() {} + template <typename TOutputWriter> template <typename TParticle, typename TTrack> - inline ProcessReturn TrackWriter::doContinuous(TParticle const& vP, TTrack const& vT, - bool const) { + inline ProcessReturn TrackWriter<TOutputWriter>::doContinuous(TParticle const& vP, + TTrack const& vT, + bool const) { + auto const start = vT.getPosition(0).getCoordinates(); - auto const delta = vT.getPosition(1).getCoordinates() - start; - auto const pdg = static_cast<int>(get_PDG(vP.getPID())); - - // clang-format off - file_ << std::setw(7) << pdg - << std::setw(width_) << std::scientific << std::setprecision(precision_) << vP.getEnergy() / 1_eV - << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[0] / 1_m - << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[1] / 1_m - << std::setw(width_) << std::scientific << std::setprecision(precision_) << start[2] / 1_m - << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[0] / 1_m - << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[1] / 1_m - << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[2] / 1_m - << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta.getNorm() / 1_m - << '\n'; - // clang-format on + auto const end = vT.getPosition(1).getCoordinates(); + + // write the track to the file + this->write(vP.getPID(), vP.getEnergy(), start, end); return ProcessReturn::Ok; } + template <typename TOutputWriter> template <typename TParticle, typename TTrack> - inline LengthType TrackWriter::getMaxStepLength(TParticle const&, TTrack const&) { + inline LengthType TrackWriter<TOutputWriter>::getMaxStepLength(TParticle const&, + TTrack const&) { return meter * std::numeric_limits<double>::infinity(); } + template <typename TOutputWriter> + YAML::Node TrackWriter<TOutputWriter>::getConfig() const { + using namespace units::si; + + YAML::Node node; + + // add default units for values + node["type"] = "TrackWriter"; + node["units"] = "GeV | m"; + + return node; + } } // namespace corsika diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl index 73fbf351707c35c894dd8cf48e0d250589909527..310613d24d2d8fcbd75402822773891b8374dfca 100644 --- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl +++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl @@ -184,9 +184,9 @@ namespace corsika { energy * 0.9, // either 10% relative loss max., or get_energy_threshold(vParticle.getPID()) // energy thresholds globally defined for // individual particles - * - 0.99 // need to go 1% below global e-cut to assure removal in ParticleCut. The - // 1% does not matter since at cut-time the entire energy is removed. + * 0.99999 // need to go slightly below global e-cut to assure removal in + // ParticleCut. This does not matter since at cut-time the entire + // energy is removed. ); auto const maxGrammage = (vParticle.getEnergy() - energy_lim) / dEdX; diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl index 217fe30b72c3fa27d1ad3794a4a6b1ec8852ffd5..9a86aa880e2333d9247bd2e18d347a83d9184541 100644 --- a/corsika/detail/modules/qgsjetII/Interaction.inl +++ b/corsika/detail/modules/qgsjetII/Interaction.inl @@ -27,20 +27,14 @@ namespace corsika::qgsjetII { - inline Interaction::Interaction(const std::string& dataPath) - : data_path_(dataPath) { - if (dataPath == "") { - if (std::getenv("CORSIKA_DATA")) { - data_path_ = std::string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/"; - CORSIKA_LOG_DEBUG("Searching for QGSJetII data tables in {}", data_path_); - } - } + inline Interaction::Interaction(boost::filesystem::path dataPath) { + CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath); // initialize QgsjetII static bool initialized = false; if (!initialized) { qgset_(); - datadir DIR(data_path_); + datadir DIR(dataPath.string() + "/"); qgaini_(DIR.data); initialized = true; } diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl index 5770fe439dbc309588e990da8669c9f85ceb619c..a6553a71dd3646534645bd8a1ee90c9fa29b08ca 100644 --- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl +++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl @@ -34,15 +34,14 @@ namespace corsika { return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV, double(0)); } // charge of the particle - int const chargeNumber = particle.getChargeNumber(); + ElectricChargeType const charge = particle.getCharge(); auto const* currentLogicalVolumeNode = particle.getNode(); MagneticFieldVector const& magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField( particle.getPosition()); - VelocityVector velocity = - particle.getMomentum() / particle.getEnergy() * constants::c; + VelocityVector velocity = particle.getVelocity(); decltype(meter / (second * volt)) k = - chargeNumber * constants::cSquared * 1_eV / + charge * constants::cSquared * 1_eV / (velocity.getNorm() * particle.getEnergy() * 1_V); DirectionVector direction = velocity.normalized(); auto position = particle.getPosition(); // First Movement @@ -58,23 +57,22 @@ namespace corsika { template <typename TParticle> inline auto Tracking::getTrack(TParticle const& particle) { - VelocityVector const initialVelocity = - particle.getMomentum() / particle.getEnergy() * constants::c; + VelocityVector const initialVelocity = particle.getVelocity(); - auto const position = particle.getPosition(); + auto const& position = particle.getPosition(); CORSIKA_LOG_DEBUG( - "Tracking pid: {}" - " , E = {} GeV", - particle.getPID(), particle.getEnergy() / 1_GeV); - CORSIKA_LOG_DEBUG("Tracking pos: {}", position.getCoordinates()); - CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV); - CORSIKA_LOG_DEBUG("Tracking p: {} GeV", - particle.getMomentum().getComponents() / 1_GeV); - CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents()); + "TrackingLeapfrog_Curved pid: {}" + " , E = {} GeV \n" + "\tTracking pos: {} \n" + "\tTracking p: {} GeV \n" + "\tTracking v: {}", + particle.getPID(), particle.getEnergy() / 1_GeV, position.getCoordinates(), + particle.getMomentum().getComponents() / 1_GeV, + initialVelocity.getComponents()); typedef typename std::remove_reference<decltype(*particle.getNode())>::type node_type; - node_type& volumeNode = *particle.getNode(); + node_type const& volumeNode = *particle.getNode(); // for the event of magnetic fields and curved trajectories, we need to limit // maximum step-length since we need to follow curved @@ -83,25 +81,28 @@ namespace corsika { MagneticFieldVector const& magneticfield = volumeNode.getModelProperties().getMagneticField(position); MagneticFluxType const magnitudeB = magneticfield.getNorm(); - int const chargeNumber = particle.getChargeNumber(); - bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T; + ElectricChargeType const charge = particle.getCharge(); + bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T; - if (no_deflection) { return getLinearTrajectory(particle); } + if (no_deflection) { + CORSIKA_LOG_TRACE("no_deflection"); + return getLinearTrajectory(particle); + } - HEPMomentumType const pAlongB_delta = + HEPMomentumType const p_perp = (particle.getMomentum() - particle.getMomentum().getParallelProjectionOnto(magneticfield)) .getNorm(); - if (pAlongB_delta == 0_GeV) { + if (p_perp < 1_eV) { // particle travel along, parallel to magnetic field. Rg is // "0", but for purpose of step limit we return infinity here. - CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel"); + CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel"); return getLinearTrajectory(particle); } - LengthType const gyroradius = - (pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV)); + LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) * + constants::c / (abs(charge) * magnitudeB)); double const maxRadians = 0.01; LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius; @@ -113,8 +114,15 @@ namespace corsika { // intersection auto [minTime, minNode] = nextIntersect(particle, steplimit_time); - auto const k = - chargeNumber * constants::cSquared * 1_eV / (particle.getEnergy() * 1_V); + auto const p_norm = + constants::c * convert_HEP_to_SI<MassType::dimension_type>( + particle.getMomentum().getNorm()); // kg *m /s + // k = q/|p| + decltype(1 / (tesla * second)) const k = + charge / p_norm * + initialVelocity.getNorm(); // since we use steps in time and not length + // units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1/(T s) + return std::make_tuple( LeapFrogTrajectory(position, initialVelocity, magneticfield, k, minTime), // trajectory @@ -129,99 +137,104 @@ namespace corsika { return Intersections(); } - int const chargeNumber = particle.getChargeNumber(); + ElectricChargeType const charge = particle.getCharge(); auto const& position = particle.getPosition(); auto const* currentLogicalVolumeNode = particle.getNode(); MagneticFieldVector const& magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(position); - VelocityVector const velocity = - particle.getMomentum() / particle.getEnergy() * constants::c; - DirectionVector const directionBefore = - velocity.normalized(); // determine steplength to next volume + VelocityVector const velocity = particle.getVelocity(); + DirectionVector const directionBefore = velocity.normalized(); auto const projectedDirection = directionBefore.cross(magneticfield); auto const projectedDirectionSqrNorm = projectedDirection.getSquaredNorm(); bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T)); - if (chargeNumber == 0 || magneticfield.getNorm() == 0_T || isParallel) { + if ((charge == 0 * constants::e) || magneticfield.getNorm() == 0_T || isParallel) { return tracking_line::Tracking::intersect<TParticle>(particle, sphere); } bool const numericallyInside = sphere.contains(particle.getPosition()); CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside); - auto const absVelocity = velocity.getNorm(); - auto const energy = particle.getEnergy(); + SpeedType const absVelocity = velocity.getNorm(); + auto const p_norm = + constants::c * convert_HEP_to_SI<MassType::dimension_type>( + particle.getMomentum().getNorm()); // km * m /s // this is: k = q/|p| - auto const k = - chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V); + auto const k = charge / p_norm; - auto const direction_B_perp = directionBefore.cross(magneticfield); - auto const denom = 4. / (direction_B_perp.getSquaredNorm() * k * k); + MagneticFieldVector const direction_x_B = directionBefore.cross(magneticfield); + auto const denom = 4. / (direction_x_B.getSquaredNorm() * k * k); Vector<length_d> const deltaPos = position - sphere.getCenter(); - double const b = (direction_B_perp.dot(deltaPos) * k + 1) * denom / (1_m * 1_m); + double const b = (direction_x_B.dot(deltaPos) * k + 1) * denom / (1_m * 1_m); double const c = directionBefore.dot(deltaPos) * 2 * denom / (1_m * 1_m * 1_m); LengthType const deltaPosLength = deltaPos.getNorm(); double const d = (deltaPosLength + sphere.getRadius()) * (deltaPosLength - sphere.getRadius()) * denom / (1_m * 1_m * 1_m * 1_m); CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d); - std::complex<double> const* solutions = quartic_solver::solve_quartic(0, b, c, d); + std::vector<double> solutions = andre::solve_quartic_real(1, 0, b, c, d); LengthType d_enter, d_exit; int first = 0, first_entry = 0, first_exit = 0; - for (int i = 0; i < 4; i++) { - if (solutions[i].imag() == 0) { - LengthType const dist = solutions[i].real() * 1_m; - CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist); - if (numericallyInside) { - // there must be an entry (negative) and exit (positive) solution - if (dist < -0.0001_m) { // security margin to assure transfer to next - // logical volume - if (first_entry == 0) { - d_enter = dist; - } else { - d_enter = std::max(d_enter, dist); // closest negative to zero (-1e-4) m - } - first_entry++; - - } else { // thus, dist > -0.0001_m - - if (first_exit == 0) { - d_exit = dist; - } else { - d_exit = std::min(d_exit, dist); // closest positive to zero (-1e-4) m - } - first_exit++; + for (auto solution : solutions) { + LengthType const dist = solution * 1_m; + CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist); + if (numericallyInside) { + // there must be an entry (negative) and exit (positive) solution + if (dist < 0.0001_m) { // security margin to assure + // transfer to next logical volume + // (even if dist suggest marginal + // entry already, which we + // classify as numerical artifact) + if (first_entry == 0) { + d_enter = dist; + } else { + d_enter = std::max(d_enter, dist); // closest negative to zero >1e-4 m } - first = int(first_exit > 0) + int(first_entry > 0); + first_entry++; - } else { // thus, numericallyInside == false + } else { // thus, dist > +0.0001_m - // both physical solutions (entry, exit) must be positive, and as small as - // possible - if (dist < -0.0001_m) { // need small numerical margin, to assure transport - // into next logical volume - continue; + if (first_exit == 0) { + d_exit = dist; + } else { + d_exit = std::min(d_exit, dist); // closest positive to zero >1e-4 m } - if (first == 0) { + first_exit++; + } + first = int(first_exit > 0) + int(first_entry > 0); + + } else { // thus, numericallyInside == false + + // both physical solutions (entry, exit) must be positive, and as small as + // possible + if (dist < -0.0001_m) { // need small numerical margin, to + // assure transport. We consider + // begin marginally already inside + // next volume (besides + // numericallyInside=false) as numerical glitch. + // into next logical volume + continue; + } + if (first == 0) { + d_enter = dist; + } else { + if (dist < d_enter) { + d_exit = d_enter; d_enter = dist; } else { - if (dist < d_enter) { - d_exit = d_enter; - d_enter = dist; - } else { - d_exit = dist; - } + d_exit = dist; } - first++; } - } // loop over solutions - } - delete[] solutions; + first++; + } + } // loop over solutions - if (first != 2) { // entry and exit points found - CORSIKA_LOG_DEBUG("no intersection! count={}", first); + if (first == 0) { // entry and exit points found + CORSIKA_LOG_DEBUG( + "no intersections found: count={}, first_entry={}, first_exit={}", first, + first_entry, first_exit); return Intersections(); } return Intersections(d_enter / absVelocity, d_exit / absVelocity); @@ -231,79 +244,82 @@ namespace corsika { inline Intersections Tracking::intersect(TParticle const& particle, Plane const& plane) { - int chargeNumber; - if (is_nucleus(particle.getPID())) { - chargeNumber = particle.getNuclearZ(); - } else { - chargeNumber = get_charge_number(particle.getPID()); - } - auto const* currentLogicalVolumeNode = particle.getNode(); - VelocityVector const velocity = - particle.getMomentum() / particle.getEnergy() * constants::c; - auto const absVelocity = velocity.getNorm(); - DirectionVector const direction = - velocity.normalized(); // determine steplength to next volume - Point const position = particle.getPosition(); - - auto const magneticfield = - currentLogicalVolumeNode->getModelProperties().getMagneticField(position); + CORSIKA_LOG_TRACE("intersection particle with plane"); - if (chargeNumber != 0 && abs(plane.getNormal().dot(velocity.cross(magneticfield))) > - 1e-6_T * 1_m / 1_s) { + ElectricChargeType const charge = particle.getCharge(); + + if (charge != 0 * constants::e) { auto const* currentLogicalVolumeNode = particle.getNode(); + VelocityVector const velocity = particle.getVelocity(); + auto const absVelocity = velocity.getNorm(); + DirectionVector const direction = velocity.normalized(); + Point const& position = particle.getPosition(); + auto const magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(position); - auto const k = - chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm(); - auto const direction_B_perp = direction.cross(magneticfield); - auto const denom = plane.getNormal().dot(direction_B_perp) * k; - auto const sqrtArg = - direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) - - (plane.getNormal().dot(position - plane.getCenter()) * denom * 2); + // solve: denom x^2 + p x + q =0 for x = delta-l + + auto const direction_x_B = direction.cross(magneticfield); + double const denom = charge * + plane.getNormal().dot(direction_x_B) // unit: C*T = kg/s + / 1_kg * 1_s; + + CORSIKA_LOG_TRACE("denom={}", denom); + + auto const p_norm = + constants::c * convert_HEP_to_SI<MassType::dimension_type>( + particle.getMomentum().getNorm()); // unit: kg * m/s + + double const p = (2 * p_norm * direction.dot(plane.getNormal())) // unit: kg*m/s + / (1_m * 1_kg) * 1_s; + double const q = + (2 * p_norm * + plane.getNormal().dot(position - plane.getCenter())) // unit: kg*m/s *m + / (1_m * 1_m * 1_kg) * 1_s; + + std::vector<double> const deltaLs = solve_quadratic_real(denom, p, q); - if (sqrtArg < 0) { + CORSIKA_LOG_TRACE("deltaLs=[{}]", fmt::join(deltaLs, ", ")); + + if (deltaLs.size() == 0) { return Intersections(std::numeric_limits<double>::infinity() * 1_s); } - double const sqrSqrtArg = sqrt(sqrtArg); - auto const norm_projected = - direction.dot(plane.getNormal()) / direction.getNorm(); - LengthType const MaxStepLength1 = (sqrSqrtArg - norm_projected) / denom; - LengthType const MaxStepLength2 = (-sqrSqrtArg - norm_projected) / denom; - CORSIKA_LOG_TRACE("MaxStepLength1={}, MaxStepLength2={}", MaxStepLength1, - MaxStepLength2); + // select smallest but positive solution + bool first = true; + LengthType maxStepLength = 0_m; + for (auto const& deltaL : deltaLs) { + if (deltaL < 0) continue; + if (first) { + first = false; + maxStepLength = deltaL * meter; + } else if (maxStepLength > deltaL * meter) { + maxStepLength = deltaL * meter; + } + } - // check: both intersections in past - if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) { + // check: both intersections in past, or no valid intersection + if (first) { return Intersections(std::numeric_limits<double>::infinity() * 1_s); - - // check: next intersection is MaxStepLength2 - } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) { - CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2); - return Intersections( - MaxStepLength2 * - (direction + direction_B_perp * MaxStepLength2 * k / 2).getNorm() / - absVelocity); - - // check: next intersections is MaxStepLength1 - } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) { - CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1); - return Intersections( - MaxStepLength1 * - (direction + direction_B_perp * MaxStepLength1 * k / 2).getNorm() / - absVelocity); } - CORSIKA_LOG_WARN( - "Particle wasn't tracked with curved trajectory -> straight (is this an " - "ERROR?)"); + CORSIKA_LOG_TRACE("maxStepLength={}", maxStepLength); + + // with final length correction + auto const corr = + (direction + direction_x_B * maxStepLength * charge / (p_norm * 2)) + .getNorm() / + absVelocity; // unit: s/m + + return Intersections(maxStepLength * corr); // unit: s - } // end if curved-tracking + } // no charge - CORSIKA_LOG_TRACE("straight tracking with chargeNumber={}, B={}", chargeNumber, - magneticfield); + CORSIKA_LOG_TRACE("(plane) straight tracking with charge={}, B={}", charge, + particle.getNode()->getModelProperties().getMagneticField( + particle.getPosition())); return tracking_line::Tracking::intersect(particle, plane); } @@ -330,7 +346,7 @@ namespace corsika { straightTrajectory.getLine().getVelocity(), MagneticFieldVector(particle.getPosition().getCoordinateSystem(), 0_T, 0_T, 0_T), - square(0_m) / (square(1_s) * 1_V), + 0 * square(meter) / (square(second) * volt), straightTrajectory.getDuration()), // trajectory minNode); // next volume node } diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl index 61ad3793d0537c47c0efb1b870d11281c663072c..4e7efed775468a6525a236fe7b5a72b45896e518 100644 --- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl +++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl @@ -29,15 +29,17 @@ namespace corsika { particle.getMomentum() / particle.getEnergy() * constants::c; Point const initialPosition = particle.getPosition(); + CORSIKA_LOG_DEBUG( - "TrackingB pid: {}" - " , E = {} GeV", - particle.getPID(), particle.getEnergy() / 1_GeV); - CORSIKA_LOG_DEBUG("TrackingB pos: {}", initialPosition.getCoordinates()); - CORSIKA_LOG_DEBUG("TrackingB E: {} GeV", particle.getEnergy() / 1_GeV); - CORSIKA_LOG_DEBUG("TrackingB p: {} GeV", - particle.getMomentum().getComponents() / 1_GeV); - CORSIKA_LOG_DEBUG("TrackingB v: {} ", initialVelocity.getComponents()); + "TrackingLeapFrogStraight pid: {}" + " , E = {} GeV \n" + "\tTracking pos: {} \n" + "\tTracking p: {} GeV \n" + "\tTracking v: {} ", + particle.getPID(), particle.getEnergy() / 1_GeV, + initialPosition.getCoordinates(), + particle.getMomentum().getComponents() / 1_GeV, + initialVelocity.getComponents()); typedef decltype(particle.getNode()) node_type; node_type const volumeNode = particle.getNode(); @@ -50,14 +52,13 @@ namespace corsika { } // charge of the particle, and magnetic field - const int chargeNumber = particle.getChargeNumber(); + auto const charge = particle.getCharge(); auto magneticfield = volumeNode->getModelProperties().getMagneticField(initialPosition); auto const magnitudeB = magneticfield.getNorm(); - CORSIKA_LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT", - magneticfield.getComponents() / 1_uT, chargeNumber, - magnitudeB / 1_T); - bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T; + CORSIKA_LOG_DEBUG("field={} uT, charge={}, magnitudeB={} uT", + magneticfield.getComponents() / 1_uT, charge, magnitudeB / 1_T); + bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T; // check, where the first halve-step direction has geometric intersections auto const [initialTrack, initialTrackNextVolume] = @@ -75,26 +76,26 @@ namespace corsika { return std::make_tuple(initialTrack, initialTrackNextVolume); } - HEPMomentumType const pAlongB_delta = + HEPMomentumType const p_perp = (particle.getMomentum() - particle.getMomentum().getParallelProjectionOnto(magneticfield)) .getNorm(); - if (pAlongB_delta == 0_GeV) { + if (p_perp == 0_GeV) { // particle travel along, parallel to magnetic field. Rg is // "0", but for purpose of step limit we return infinity here. - CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel"); + CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel"); return std::make_tuple(initialTrack, initialTrackNextVolume); } - LengthType const gyroradius = - (pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV)); + LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) * + constants::c / (abs(charge) * magnitudeB)); // we need to limit maximum step-length since we absolutely // need to follow strongly curved trajectories segment-wise, // at least if we don't employ concepts as "Helix // Trajectories" or similar - double const maxRadians = 0.01; + double const maxRadians = 0.001; LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius; CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit); @@ -111,9 +112,9 @@ namespace corsika { firstHalveSteplength, steplimit, initialTrackLength); // perform the first halve-step Point const position_mid = initialPosition + direction * firstHalveSteplength; - auto const k = - chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm(); - auto const new_direction = + auto const k = charge / (constants::c * convert_HEP_to_SI<MassType::dimension_type>( + particle.getMomentum().getNorm())); + DirectionVector const new_direction = direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k; auto const new_direction_norm = new_direction.getNorm(); // by design this is >1 CORSIKA_LOG_DEBUG( diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl index 9a4659250286414cc88f837e96c3ac5898952006..cad8ac39d99ebf67af35c19eea7974897806b3f9 100644 --- a/corsika/detail/modules/tracking/TrackingStraight.inl +++ b/corsika/detail/modules/tracking/TrackingStraight.inl @@ -20,6 +20,7 @@ #include <type_traits> #include <utility> +#include <cmath> namespace corsika::tracking_line { @@ -28,16 +29,15 @@ namespace corsika::tracking_line { VelocityVector const initialVelocity = particle.getMomentum() / particle.getEnergy() * constants::c; - auto const initialPosition = particle.getPosition(); + auto const& initialPosition = particle.getPosition(); CORSIKA_LOG_DEBUG( - "Tracking pid: {}" - " , E = {} GeV", - particle.getPID(), particle.getEnergy() / 1_GeV); - CORSIKA_LOG_DEBUG("Tracking pos: {}", initialPosition.getCoordinates()); - CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV); - CORSIKA_LOG_DEBUG("Tracking p: {} GeV", - particle.getMomentum().getComponents() / 1_GeV); - CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents()); + "TrackingStraight pid: {}" + " , E = {} GeV \n" + "\tTracking pos: {} \n" + "\tTracking p: {} GeV \n" + "\tTracking v: {}", + particle.getPID(), particle.getEnergy() / 1_GeV, initialPosition.getCoordinates(), + particle.getMomentum().getComponents() / 1_GeV, initialVelocity.getComponents()); // traverse the environment volume tree and find next // intersection @@ -51,7 +51,8 @@ namespace corsika::tracking_line { template <typename TParticle> inline Intersections Tracking::intersect(TParticle const& particle, Sphere const& sphere) { - auto const delta = particle.getPosition() - sphere.getCenter(); + auto const& position = particle.getPosition(); + auto const delta = position - sphere.getCenter(); auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c; auto const vSqNorm = velocity.getSquaredNorm(); auto const R = sphere.getRadius(); @@ -63,6 +64,8 @@ namespace corsika::tracking_line { if (discriminant.magnitude() > 0) { auto const sqDisc = sqrt(discriminant); auto const invDenom = 1 / vSqNorm; + + CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position)); return Intersections((-vDotDelta - sqDisc) * invDenom, (-vDotDelta + sqDisc) * invDenom); } diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl index cd5a479bab9a2eda886519e0c75ffe4fce8b7dcc..e6572b37bff6bfcd7e7115a4dbb9ce215f244e35 100644 --- a/corsika/detail/modules/urqmd/UrQMD.inl +++ b/corsika/detail/modules/urqmd/UrQMD.inl @@ -29,7 +29,7 @@ namespace corsika::urqmd { - inline UrQMD::UrQMD(boost::filesystem::path const& xs_file) { + inline UrQMD::UrQMD(boost::filesystem::path xs_file) { readXSFile(xs_file); ::urqmd::iniurqmdc8_(); } @@ -268,10 +268,8 @@ namespace corsika::urqmd { 2 * ::urqmd::options_.CTParam[30 - 1]; ::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV); - if (projectileCode == Code::K0Long) { + if (projectileCode == Code::K0Long || projectileCode == Code::K0Short) { projectileCode = booleanDist_(RNG_) ? Code::K0 : Code::K0Bar; - } else if (projectileCode == Code::K0Short) { - throw std::runtime_error("K0Short should not interact"); } auto const [ityp, iso3] = convertToUrQMD(projectileCode); @@ -396,7 +394,7 @@ namespace corsika::urqmd { return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code))); } - inline void UrQMD::readXSFile(boost::filesystem::path const& filename) { + inline void UrQMD::readXSFile(boost::filesystem::path const filename) { boost::filesystem::ifstream file(filename, std::ios::in); if (!file.is_open()) { @@ -431,6 +429,7 @@ namespace corsika::urqmd { std::getline(file, line); } } + file.close(); } } // namespace corsika::urqmd diff --git a/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl b/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl new file mode 100644 index 0000000000000000000000000000000000000000..60752c56a9fed06b326449530c2d24d5837d233e --- /dev/null +++ b/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl @@ -0,0 +1,54 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +namespace corsika { + + inline ObservationPlaneWriterParquet::ObservationPlaneWriterParquet() + : output_() {} + + inline void ObservationPlaneWriterParquet::startOfLibrary( + boost::filesystem::path const& directory) { + + // setup the streamer + output_.initStreamer((directory / "particles.parquet").string()); + + // enable compression with the default level + output_.enableCompression(); + + // build the schema + output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32, + parquet::ConvertedType::INT_32); + output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + + // and build the streamer + output_.buildStreamer(); + } + + inline void ObservationPlaneWriterParquet::endOfShower() { ++shower_; } + + inline void ObservationPlaneWriterParquet::endOfLibrary() { output_.closeStreamer(); } + + inline void ObservationPlaneWriterParquet::write(Code const& pid, + HEPEnergyType const& energy, + LengthType const& x, + LengthType const& y) { + // write the next row - we must write `shower_` first. + *(output_.getWriter()) << shower_ << static_cast<int>(get_PDG(pid)) + << static_cast<float>(energy / 1_GeV) + << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m) + << parquet::EndRow; + } + +} // namespace corsika diff --git a/corsika/detail/modules/writers/TrackWriterParquet.inl b/corsika/detail/modules/writers/TrackWriterParquet.inl new file mode 100644 index 0000000000000000000000000000000000000000..26c62bdf0e5e8260c431dc114aaed25e4aab0a82 --- /dev/null +++ b/corsika/detail/modules/writers/TrackWriterParquet.inl @@ -0,0 +1,68 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +namespace corsika { + + inline TrackWriterParquet::TrackWriterParquet() + : output_() {} + + inline void TrackWriterParquet::startOfLibrary( + boost::filesystem::path const& directory) { + + // setup the streamer + output_.initStreamer((directory / "tracks.parquet").string()); + + // build the schema + output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32, + parquet::ConvertedType::INT_32); + output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("start_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("start_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("start_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("end_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("end_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("end_z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + + // and build the streamer + output_.buildStreamer(); + } + + inline void TrackWriterParquet::endOfShower() { ++shower_; } + + inline void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); } + + inline void TrackWriterParquet::write(Code const& pid, HEPEnergyType const& energy, + QuantityVector<length_d> const& start, + QuantityVector<length_d> const& end) { + + // write the next row - we must write `shower_` first. + // clang-format off + *(output_.getWriter()) + << shower_ + << static_cast<int>(get_PDG(pid)) + << static_cast<float>(energy / 1_GeV) + << static_cast<float>(start[0] / 1_m) + << static_cast<float>(start[1] / 1_m) + << static_cast<float>(start[2] / 1_m) + << static_cast<float>(end[0] / 1_m) + << static_cast<float>(end[1] / 1_m) + << static_cast<float>(end[2] / 1_m) + << parquet::EndRow; + // clang-format on + } + +} // namespace corsika diff --git a/corsika/detail/output/BaseOutput.inl b/corsika/detail/output/BaseOutput.inl new file mode 100644 index 0000000000000000000000000000000000000000..24ec861ac9c2b3083ea834261f7fa0a57fd8559e --- /dev/null +++ b/corsika/detail/output/BaseOutput.inl @@ -0,0 +1,15 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +namespace corsika { + + inline BaseOutput::BaseOutput() + : shower_(0) {} + +} // namespace corsika diff --git a/corsika/detail/output/DummyOutputManager.inl b/corsika/detail/output/DummyOutputManager.inl new file mode 100644 index 0000000000000000000000000000000000000000..acf2b69f987dbebeb69c9a77438141afc36a036d --- /dev/null +++ b/corsika/detail/output/DummyOutputManager.inl @@ -0,0 +1,27 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +namespace corsika { + + inline DummyOutputManager::DummyOutputManager() {} + + inline DummyOutputManager::~DummyOutputManager() {} + + template <typename TOutput> + inline void DummyOutputManager::add(std::string const&, TOutput&) {} + + inline void DummyOutputManager::startOfLibrary() {} + + inline void DummyOutputManager::startOfShower() {} + + inline void DummyOutputManager::endOfShower() {} + + inline void DummyOutputManager::endOfLibrary() {} + +} // namespace corsika diff --git a/corsika/detail/output/OutputManager.inl b/corsika/detail/output/OutputManager.inl new file mode 100644 index 0000000000000000000000000000000000000000..6a4ec62045acddbf9028f27d9991239c63012553 --- /dev/null +++ b/corsika/detail/output/OutputManager.inl @@ -0,0 +1,238 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +#include <algorithm> +#include <fstream> +#include <functional> + +#include <iomanip> +#include <ctime> +#include <sstream> + +#include <boost/filesystem.hpp> + +#include <fmt/core.h> +#include <fmt/chrono.h> + +namespace corsika { + + inline void OutputManager::writeNode(YAML::Node const& node, + boost::filesystem::path const& path) const { + + // construct a YAML emitter for this config file + YAML::Emitter out; + + // and write the node to the output + out << node; + + // open the output file - this is <output name>.yaml + boost::filesystem::ofstream file(path); + + // dump the YAML to the file + file << out.c_str() << std::endl; + } + + inline void OutputManager::writeTopLevelConfig() const { + + YAML::Node config; + + // some basic info + config["name"] = name_; // the simulation name + config["creator"] = "CORSIKA8"; // a tag to identify C8 libraries + config["version"] = "8.0.0-prealpha"; // the current version + + // write the node to a file + writeNode(config, root_ / ("config.yaml")); + } + + inline void OutputManager::writeTopLevelSummary() const { + + YAML::Node config; + + // the total number of showers contained in the library + config["showers"] = count_; + + // this next section handles writing some time and duration information + + // create a quick lambda function to convert a time-instance to a string + auto timeToString = [&](auto const time) -> std::string { + // ISO 8601 time format + auto format{"%FT%T%z"}; + + // convert the clock to a time_t + auto time_tc{std::chrono::system_clock::to_time_t(time)}; + + // create the string and push the time onto it + std::ostringstream oss; + oss << std::put_time(std::localtime(&time_tc), format); + + return oss.str(); + }; + + auto end_time{std::chrono::system_clock::now()}; + + // now let's construct an estimate of the runtime + auto runtime{end_time - start_time}; + + // add the time and duration info + config["start time"] = timeToString(start_time); + config["end time"] = timeToString(end_time); + config["runtime"] = fmt::format("{:%H:%M:%S}", runtime); + + // write the node to a file + writeNode(config, root_ / ("summary.yaml")); + } + + inline void OutputManager::initOutput(std::string const& name) const { + // construct the path to this directory + auto const path{root_ / name}; + + // create the directory for this process. + boost::filesystem::create_directory(path); + + // get the config for this output + auto config = outputs_.at(name).get().getConfig(); + + // and assign the name for this output + config["name"] = name; + + // write the config for this output to the file + writeNode(config, path / "config.yaml"); + } + + inline OutputManager::OutputManager( + std::string const& name, + boost::filesystem::path const& dir = boost::filesystem::current_path()) + : name_(name) + , root_(dir / name) { + + // check if this directory already exists + if (boost::filesystem::exists(root_)) { + logger->error("Output directory '{}' already exists! Do not overwrite!.", + root_.string()); + throw std::runtime_error("Output directory already exists."); + } + + // construct the directory for this library + boost::filesystem::create_directories(root_); + + // write the top level config file + writeTopLevelConfig(); + } + + inline OutputManager::~OutputManager() { + + if (state_ == OutputState::ShowerInProgress) { + // if this the destructor is called before the shower has been explicitly + // ended, print a warning and end the shower before continuing. + logger->warn( + "OutputManager was destroyed before endOfShower() called." + " The last shower in this libray may be incomplete."); + endOfShower(); + } + + // write the top level summary file (summary.yaml) + writeTopLevelSummary(); + + // if we are being destructed but EndOfLibrary() has not been called, + // make sure that we gracefully close all the outputs. This is a supported + // method of operation so we don't issue a warning here + if (state_ == OutputState::LibraryReady) { endOfLibrary(); } + } + + template <typename TOutput> + inline void OutputManager::add(std::string const& name, TOutput& output) { + + // check if that name is already in the map + if (outputs_.count(name) > 0) { + logger->error("'{}' is already registered. All outputs must have unique names.", + name); + throw std::runtime_error("Output already exists. Do not overwrite!"); + } + + // if we get here, the name is not already in the map + // so we create the output and register it into the map + outputs_.insert(std::make_pair(name, std::ref(output))); + + // and initialize this output + initOutput(name); + } + + inline void OutputManager::startOfLibrary() { + + // this is only valid when we haven't started a library + // or have already finished a library + if (!(state_ == OutputState::NoInit || state_ == OutputState::LibraryFinished)) { + + throw std::runtime_error("startOfLibrary() called in invalid state."); + } + + // we now forward this signal to all of our outputs + for (auto& [name, output] : outputs_) { + + // construct the path to this output subdirectory + auto const path{root_ / name}; + + // and start the library + output.get().startOfLibrary(path); + } + + // we have now started running + state_ = OutputState::LibraryReady; + } + + inline void OutputManager::startOfShower() { + + // if this is called and we still in the "no init" state, then + // this is the first shower in the library so make sure we start it + if (state_ == OutputState::NoInit) { startOfLibrary(); } + + // now start the event for all the outputs + for (auto& [name, output] : outputs_) { output.get().startOfShower(); } + + // increment our shower count + ++count_; + + // and transition to the in progress state + state_ = OutputState::ShowerInProgress; + } + + inline void OutputManager::endOfShower() { + + for (auto& [name, output] : outputs_) { output.get().endOfShower(); } + + // switch back to the initialized state + state_ = OutputState::LibraryReady; + } + + inline void OutputManager::endOfLibrary() { + + // we can only call endOfLibrary when we have already started + if (state_ == OutputState::NoInit) { + throw std::runtime_error("endOfLibrary() called in invalid state."); + } + + // write the summary for each output and forward the endOfLibrary call() + for (auto& [name, output] : outputs_) { + + // we get the summary for each output as a YAML node + auto summary{outputs_.at(name).get().getSummary()}; + + // write the summary for this output to the file + writeNode(summary, root_ / name / "summary.yaml"); + + // and forward the end of library call + output.get().endOfLibrary(); + } + + // and the library has finished + state_ = OutputState::LibraryFinished; + } + +} // namespace corsika diff --git a/corsika/detail/output/ParquetStreamer.inl b/corsika/detail/output/ParquetStreamer.inl new file mode 100644 index 0000000000000000000000000000000000000000..d0bb2f9e01a42fef34c00fd2d41c189e60f252e2 --- /dev/null +++ b/corsika/detail/output/ParquetStreamer.inl @@ -0,0 +1,71 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +namespace corsika { + + inline ParquetStreamer::ParquetStreamer() + : isInit_(false) {} + + inline void ParquetStreamer::initStreamer(std::string const& filepath) { + + // open the file and connect it to our pointer + PARQUET_ASSIGN_OR_THROW(outfile_, arrow::io::FileOutputStream::Open(filepath)); + + // the default builder settings + builder_.created_by("CORSIKA8"); + + // add run and event tags to the file + addField("shower", parquet::Repetition::REQUIRED, parquet::Type::INT32, + parquet::ConvertedType::INT_32); + } + + template <typename... TArgs> + inline void ParquetStreamer::addField(TArgs&&... args) { + fields_.push_back(parquet::schema::PrimitiveNode::Make(args...)); + } + + inline void ParquetStreamer::enableCompression(int const /*level*/) { + // builder_.compression(parquet::Compression::ZSTD); + // builder_.compression_level(level); + } + + inline void ParquetStreamer::buildStreamer() { + + // build the top level schema + schema_ = std::static_pointer_cast<parquet::schema::GroupNode>( + parquet::schema::GroupNode::Make("schema", parquet::Repetition::REQUIRED, + fields_)); + + // and build the writer + writer_ = std::make_shared<parquet::StreamWriter>( + parquet::ParquetFileWriter::Open(outfile_, schema_, builder_.build())); + + // only now this object is ready to stream + isInit_ = true; + } + + inline void ParquetStreamer::closeStreamer() { + writer_.reset(); + [[maybe_unused]] auto status = outfile_->Close(); + isInit_ = false; + } + + inline std::shared_ptr<parquet::StreamWriter> ParquetStreamer::getWriter() { + if (!isInit()) { + throw std::runtime_error( + "ParquetStreamer not initialized. Either 1) add the " + "corresponding module to " + "the OutputManager, or 2) declare the module to write no output using " + "NoOutput."); + } + return writer_; + } + +} // namespace corsika diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl index 728647ca0bfc5ea89738946ccdef9c36cda31ce2..5dab28ab2506181580fd510e28164c507b70b878 100644 --- a/corsika/detail/stack/NuclearStackExtension.inl +++ b/corsika/detail/stack/NuclearStackExtension.inl @@ -113,6 +113,14 @@ namespace corsika::nuclear_stack { return super_type::getMass(); } + template <template <typename> class InnerParticleInterface, + typename StackIteratorInterface> + inline ElectricChargeType NuclearParticleInterface< + InnerParticleInterface, StackIteratorInterface>::getCharge() const { + if (super_type::getPID() == Code::Nucleus) return getNuclearZ() * constants::e; + return super_type::getCharge(); + } + template <template <typename> class InnerParticleInterface, typename StackIteratorInterface> inline int16_t NuclearParticleInterface< diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp index 7ee12186ace6f189b2c375b8e1a04e502a3cf06f..1eadbce1160a405f6cc1738d62997a20923aab3a 100644 --- a/corsika/framework/core/Cascade.hpp +++ b/corsika/framework/core/Cascade.hpp @@ -57,7 +57,7 @@ namespace corsika { * * */ - template <typename TTracking, typename TProcessList, typename TStack, + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack, /* TStackView is needed as explicit template parameter because of issue 161 and the @@ -84,10 +84,11 @@ namespace corsika { ~Cascade() = default; Cascade& operator=(Cascade const&) = default; Cascade(Environment<MediumInterface> const& env, TTracking& tr, TProcessList& pl, - TStack& stack) + TOutput& out, TStack& stack) : environment_(env) , tracking_(tr) , sequence_(pl) + , output_(out) , stack_(stack) { CORSIKA_LOG_INFO(c8_ascii_); CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(), @@ -130,14 +131,16 @@ namespace corsika { */ void step(Particle& vParticle); - ProcessReturn decay(TStackView& view); - ProcessReturn interaction(TStackView& view); + ProcessReturn decay(TStackView& view, InverseTimeType initial_inv_decay_time); + ProcessReturn interaction(TStackView& view, + InverseGrammageType initial_inv_int_length); void setEventType(TStackView& view, history::EventType); // data members Environment<MediumInterface> const& environment_; TTracking& tracking_; TProcessList& sequence_; + TOutput& output_; TStack& stack_; default_prng_type& rng_ = RNGManager::getInstance().getRandomStream("cascade"); unsigned int count_ = 0; diff --git a/corsika/framework/core/Logging.hpp b/corsika/framework/core/Logging.hpp index 67f059c2fe6d1dc5d4b228e551f4ebbbcd1acf9e..d127040b036dab150ce06a40644050f20a3fab2a 100644 --- a/corsika/framework/core/Logging.hpp +++ b/corsika/framework/core/Logging.hpp @@ -49,7 +49,8 @@ namespace corsika { /* * The default pattern for CORSIKA8 loggers. */ - const std::string default_pattern{"[%n:%^%-8l%$] %v"}; + const std::string minimal_pattern{"[%n:%^%-8l%$] %v"}; + const std::string default_pattern{"[%n:%^%-8l%$(%s:%#)] %v"}; const std::string source_pattern{"[%n:%^%-8l%$(%s:%!:%#)] %v"}; /** @@ -91,7 +92,6 @@ namespace corsika { */ static inline std::shared_ptr<spdlog::logger> corsika_logger = get_logger("corsika", true); - // corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); // many of these free functions are special to the logging // infrastructure so we hide them in the corsika::logging namespace. diff --git a/corsika/framework/geometry/LeapFrogTrajectory.hpp b/corsika/framework/geometry/LeapFrogTrajectory.hpp index 6b67bf6345a5c79400d72d9ec16e13afa380e13f..11e3534a19ae6c8ead194077ae515afb194bbd4c 100644 --- a/corsika/framework/geometry/LeapFrogTrajectory.hpp +++ b/corsika/framework/geometry/LeapFrogTrajectory.hpp @@ -39,7 +39,7 @@ namespace corsika { LeapFrogTrajectory(Point const& pos, VelocityVector const& initialVelocity, MagneticFieldVector const& Bfield, - decltype(square(meter) / (square(second) * volt)) const k, + decltype(1 / (tesla * second)) const k, TimeType const timeStep) // leap-from total length : initialPosition_(pos) , initialVelocity_(initialVelocity) @@ -74,7 +74,7 @@ namespace corsika { VelocityVector initialVelocity_; DirectionVector initialDirection_; MagneticFieldVector magneticfield_; - decltype(square(meter) / (square(second) * volt)) k_; + decltype(1 / (tesla * second)) k_; TimeType timeStep_; }; diff --git a/corsika/framework/geometry/Path.hpp b/corsika/framework/geometry/Path.hpp new file mode 100644 index 0000000000000000000000000000000000000000..620cbeeab2210141e8a51f423e9422688313e4d2 --- /dev/null +++ b/corsika/framework/geometry/Path.hpp @@ -0,0 +1,85 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <deque> +#include <corsika/framework/geometry/Point.hpp> + +namespace corsika { + + /** + * This class represents a (potentially) curved path between two + * points using N >= 1 straight-line segments. + */ + class Path { + std::deque<Point> points_; ///< The points that make up this path. + LengthType length_ = LengthType::zero(); ///< The length of the path. + public: + /** + * Create a Path with a given starting Point. + */ + Path(Point const& point); + + /** + * Initialize a Path from an existing collection of Points. + */ + Path(std::deque<Point> const& points); + + /** + * Add a new Point to the end of the path. + */ + inline void addToEnd(Point const& point); + + /** + * Remove a point from the end of the path. + */ + inline void removeFromEnd(); + + /** + * Get the total length of the path. + */ + inline LengthType getLength() const; + + /** + * Get the starting point of the path. + */ + inline Point getStart() const; + + /** + * Get the end point of the path. + */ + inline Point getEnd() const; + + /** + * Get a specific point of the path. + */ + inline Point getPoint(std::size_t const index) const; + + /** + * Return an iterator to the start of the Path. + */ + inline auto begin(); + + /** + * Return an iterator to the end of the Path. + */ + inline auto end(); + + /** + * Get the number of steps in the path. + * This is one less than the number of points that + * defines the path. + */ + inline int getNSegments() const; + + }; // class Path + +} // namespace corsika + +#include <corsika/detail/framework/geometry/Path.inl> \ No newline at end of file diff --git a/corsika/framework/geometry/Point.hpp b/corsika/framework/geometry/Point.hpp index 2c684438438455985d0bebfbb95ed6bdafa4c86b..6eb8786c9df6accf7e654585a5bd4a8d8a957313 100644 --- a/corsika/framework/geometry/Point.hpp +++ b/corsika/framework/geometry/Point.hpp @@ -31,8 +31,8 @@ namespace corsika { /** \todo TODO: this should be private or protected, we don NOT want to expose numbers * without reference to outside: */ - QuantityVector<length_d> const& getCoordinates() const; - QuantityVector<length_d>& getCoordinates(); + inline QuantityVector<length_d> const& getCoordinates() const; + inline QuantityVector<length_d>& getCoordinates(); /** this always returns a QuantityVector as triple @@ -40,7 +40,7 @@ namespace corsika { \returns A value type QuantityVector, since it may have to create a temporary object to transform to pCS. **/ - QuantityVector<length_d> getCoordinates(CoordinateSystemPtr const& pCS) const; + inline QuantityVector<length_d> getCoordinates(CoordinateSystemPtr const& pCS) const; /** * this always returns a QuantityVector as triple @@ -49,7 +49,7 @@ namespace corsika { * is actually transformed to pCS, if needed. Thus, there may be an implicit call to * \ref rebase. **/ - QuantityVector<length_d>& getCoordinates(CoordinateSystemPtr const& pCS); + inline QuantityVector<length_d>& getCoordinates(CoordinateSystemPtr const& pCS); /** * \name access coordinate components @@ -60,25 +60,30 @@ namespace corsika { * created and destroyed each call. This can be avoided by using * \ref rebase first. **/ - LengthType getX(CoordinateSystemPtr const& pCS) const; - LengthType getY(CoordinateSystemPtr const& pCS) const; - LengthType getZ(CoordinateSystemPtr const& pCS) const; + inline LengthType getX(CoordinateSystemPtr const& pCS) const; + inline LengthType getY(CoordinateSystemPtr const& pCS) const; + inline LengthType getZ(CoordinateSystemPtr const& pCS) const; /** \} **/ /*! * transforms the Point into another CoordinateSystem by changing its * coordinates interally */ - void rebase(CoordinateSystemPtr const& pCS); + inline void rebase(CoordinateSystemPtr const& pCS); - Point operator+(Vector<length_d> const& pVec) const; + inline Point operator+(Vector<length_d> const& pVec) const; /*! * returns the distance Vector between two points */ - Vector<length_d> operator-(Point const& pB) const; + inline Vector<length_d> operator-(Point const& pB) const; }; + /* + * calculates the distance between two points + */ + inline LengthType distance(Point const& p1, Point const& p2); + } // namespace corsika -#include <corsika/detail/framework/geometry/Point.inl> +#include <corsika/detail/framework/geometry/Point.inl> \ No newline at end of file diff --git a/corsika/framework/process/ContinuousProcess.hpp b/corsika/framework/process/ContinuousProcess.hpp index c549bc714da8e4c6b655afbb8690568a2f4c4e24..51a7ea7109439401c23e39b722fdd731ff8b1ddb 100644 --- a/corsika/framework/process/ContinuousProcess.hpp +++ b/corsika/framework/process/ContinuousProcess.hpp @@ -18,47 +18,46 @@ namespace corsika { /** - @ingroup Processes - @{ - - Processes with continuous effects along a particle Trajectory - - Create a new ContinuousProcess, e.g. for XYModel, via - @code{.cpp} - class XYModel : public ContinuousProcess<XYModel> {}; - @endcode - - and provide two necessary interface methods: - @code{.cpp} - template <typename TParticle, typename TTrack> - LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const; - @endcode - - which allows any ContinuousProcess to tell to CORSIKA a maximum - allowed step length. Such step-length limitation, if it turns out - to be smaller/sooner than any other limit (decay length, - interaction length, other continuous processes, geometry, etc.) - will lead to a limited step length. - - @code{.cpp} - template <typename TParticle, typename TTrack> - ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit) - const; - @endcode - - which applied any continuous effects on Particle p along - Trajectory t. The particle in all typical scenarios will be - altered by a doContinuous. The flag stepLimit will be true if the - preious evaluation of getMaxStepLength resulted in this - particular ContinuousProcess to be responsible for the step - length limit on the current track t. This information can be - expoited and avoid e.g. any uncessary calculations. - - Particle and Track are the valid classes to - access particles and track (Trajectory) data on the Stack. Those two methods - do not need to be templated, they could use the types - e.g. corsika::setup::Stack::particle_type -- but by the cost of - loosing all flexibility otherwise provided. + * @ingroup Processes + * @{ + * Processes with continuous effects along a particle Trajectory. + * + * Create a new ContinuousProcess, e.g. for XYModel, via: + * @code{.cpp} + * class XYModel : public ContinuousProcess<XYModel> {}; + * @endcode + * + * and provide two necessary interface methods: + * @code{.cpp} + * template <typename TParticle, typename TTrack> + * LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const; + * @endcode + + * which allows any ContinuousProcess to tell to CORSIKA a maximum + * allowed step length. Such step-length limitation, if it turns out + * to be smaller/sooner than any other limit (decay length, + * interaction length, other continuous processes, geometry, etc.) + * will lead to a limited step length. + + * @code{.cpp} + * template <typename TParticle, typename TTrack> + * ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit) + * const; + * @endcode + + * which applied any continuous effects on Particle p along + * Trajectory t. The particle in all typical scenarios will be + * altered by a doContinuous. The flag stepLimit will be true if the + * preious evaluation of getMaxStepLength resulted in this + * particular ContinuousProcess to be responsible for the step + * length limit on the current track t. This information can be + * expoited and avoid e.g. any uncessary calculations. + + * Particle and Track are the valid classes to + * access particles and track (Trajectory) data on the Stack. Those two methods + * do not need to be templated, they could use the types + * e.g. corsika::setup::Stack::particle_type -- but by the cost of + * loosing all flexibility otherwise provided. */ @@ -68,8 +67,8 @@ namespace corsika { }; /** - * ProcessTraits specialization to flag ContinuousProcess objects - **/ + * ProcessTraits specialization to flag ContinuousProcess objects. + */ template <typename TProcess> struct is_continuous_process< TProcess, std::enable_if_t< diff --git a/corsika/framework/process/ContinuousProcessIndex.hpp b/corsika/framework/process/ContinuousProcessIndex.hpp index e36bf2812a33c6f03e4058033585f846a871b897..7192b86bf968ffdc74b5420940ea549befa12c73 100644 --- a/corsika/framework/process/ContinuousProcessIndex.hpp +++ b/corsika/framework/process/ContinuousProcessIndex.hpp @@ -11,12 +11,10 @@ namespace corsika { /** - @ingroup Processes - - To index individual processes (continuous processes) inside a - ProcessSequence. - - **/ + * @ingroup Processes + * To index individual processes (continuous processes) inside a + * ProcessSequence. + */ class ContinuousProcessIndex { public: diff --git a/corsika/framework/process/ContinuousProcessStepLength.hpp b/corsika/framework/process/ContinuousProcessStepLength.hpp index 575be9063d49a32b5853320dc9d25bac40a7b9d8..55c2f417a470b8e4a475de8269c5e19fe8247ff5 100644 --- a/corsika/framework/process/ContinuousProcessStepLength.hpp +++ b/corsika/framework/process/ContinuousProcessStepLength.hpp @@ -14,12 +14,10 @@ namespace corsika { /** - @ingroup Processes - - To store step length in LengthType and unique index in ProcessSequence of shortest - step ContinuousProcess. - - **/ + * @ingroup Processes + * To store step length in LengthType and unique index in ProcessSequence of shortest + * step ContinuousProcess. + */ class ContinuousProcessStepLength { public: diff --git a/corsika/framework/utility/CubicSolver.hpp b/corsika/framework/utility/CubicSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..0959a4fa84c3bee6aad856ccb57c887d21e30c9c --- /dev/null +++ b/corsika/framework/utility/CubicSolver.hpp @@ -0,0 +1,93 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <vector> +#include <cmath> +#include <algorithm> +#include <numeric> +#include <complex> + +#include <corsika/framework/utility/QuadraticSolver.hpp> +#include <corsika/framework/core/Logging.hpp> + +#include <boost/multiprecision/cpp_bin_float.hpp> + +/** + * @file CubicSolver.hpp + */ + +namespace corsika { + + /** + * @ingroup Utility + * @{ + */ + + namespace andre { + /** + Solve a x^3 + b x^2 + c x + d = 0 + */ + std::vector<double> solveP3(long double a, long double b, long double c, + long double d, double const epsilon = 1e-12); + } // namespace andre + + /** + Solve depressed cubic: x^3 + p x + q = 0 + */ + inline std::vector<double> solve_cubic_depressed_disciminant_real( + long double p, long double q, long double const disc, double const epsilon = 1e-12); + + std::vector<double> solve_cubic_depressed_real(long double p, long double q, + double const epsilon = 1e-12); + + /** + Solve a x^3 + b x^2 + c x + d = 0 + + Analytical approach. Not very stable in all conditions. + */ + std::vector<double> solve_cubic_real_analytic(long double a, long double b, + long double c, long double d, + double const epsilon = 1e-12); + + /** + * Cubic function. + * + * T must be a floating point type. + * + * @return a x^3 + b x^2 + c x + d + */ + template <typename T> + T cubic_function(T x, T a, T b, T c, T d); + + /** + @return 3 a x^2 + 2 b x + c + */ + template <typename T> + T cubic_function_dfdx(T x, T a, T b, T c); + + /** + @return 6 a x + 2 b + */ + template <typename T> + T cubic_function_d2fd2x(T x, T a, T b); + + /** + Iterative approach to solve: a x^3 + b x^2 + c x + d = 0 + + This is often fastest and most precise. + */ + std::vector<double> solve_cubic_real(long double a, long double b, long double c, + long double d, double const epsilon = 1e-12); + + //! @} + +} // namespace corsika + +#include <corsika/detail/framework/utility/CubicSolver.inl> diff --git a/corsika/framework/utility/LinearSolver.hpp b/corsika/framework/utility/LinearSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..754039571287a1a816c77ec741b7d850b87c46da --- /dev/null +++ b/corsika/framework/utility/LinearSolver.hpp @@ -0,0 +1,23 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <vector> +#include <cmath> +#include <complex> + +namespace corsika { + + std::vector<double> solve_linear_real(double a, double b); + + std::vector<std::complex<double>> solve_linear(double a, double b); + +} // namespace corsika + +#include <corsika/detail/framework/utility/LinearSolver.inl> diff --git a/corsika/framework/utility/QuadraticSolver.hpp b/corsika/framework/utility/QuadraticSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..968c01bf13e3680b91b27005084c909b4e473ea7 --- /dev/null +++ b/corsika/framework/utility/QuadraticSolver.hpp @@ -0,0 +1,27 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <vector> +#include <cmath> +#include <complex> + +#include <corsika/framework/utility/LinearSolver.hpp> + +namespace corsika { + + std::vector<std::complex<double>> solve_quadratic(long double a, long double b, + long double c, + double const epsilon = 1e-12); + + std::vector<double> solve_quadratic_real(long double a, long double b, long double c, + double const epsilon = 1e-12); +} // namespace corsika + +#include <corsika/detail/framework/utility/QuadraticSolver.inl> diff --git a/corsika/framework/utility/QuarticSolver.hpp b/corsika/framework/utility/QuarticSolver.hpp index ad9484939d10691e38c7e74d9e2632b3d2e60463..d6e20f02fc3d6d594a64475d9ce1117b4dec159c 100644 --- a/corsika/framework/utility/QuarticSolver.hpp +++ b/corsika/framework/utility/QuarticSolver.hpp @@ -8,59 +8,46 @@ #pragma once +#include <corsika/framework/utility/CubicSolver.hpp> +#include <corsika/framework/utility/QuadraticSolver.hpp> + #include <complex> +#include <algorithm> +#include <cmath> +#include <numeric> /** * @file QuarticSolver.hpp - * @ingroup Utilities - * @{ - * \todo convert to class */ -namespace corsika::quartic_solver { - - const double PI = 3.141592653589793238463L; - const double M_2PI = 2 * PI; - const double eps = 1e-12; - - typedef std::complex<double> DComplex; - - //--------------------------------------------------------------------------- - // useful for testing - inline DComplex polinom_2(DComplex x, double a, double b) { - // Horner's scheme for x*x + a*x + b - return x * (x + a) + b; - } +namespace corsika { - //--------------------------------------------------------------------------- - // useful for testing - inline DComplex polinom_3(DComplex x, double a, double b, double c) { - // Horner's scheme for x*x*x + a*x*x + b*x + c; - return x * (x * (x + a) + b) + c; - } + /** + * @ingroup Utilities + * @{ + */ - //--------------------------------------------------------------------------- - // useful for testing - inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) { - // Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d; - return x * (x * (x * (x + a) + b) + c) + d; - } + namespace andre { - //--------------------------------------------------------------------------- - // x - array of size 3 - // In case 3 real roots: => x[0], x[1], x[2], return 3 - // 2 real roots: x[0], x[1], return 2 - // 1 real root : x[0], x[1] ± i*x[2], return 1 - unsigned int solveP3(double* x, double a, double b, double c); + /** + solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e + Attention - this function returns dynamically allocated array. It has to be + released afterwards. + */ + std::vector<double> solve_quartic_real(long double a, long double b, long double c, + long double d, long double e, + double const epsilon = 1e-12); + } // namespace andre - //--------------------------------------------------------------------------- - // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d - // Attention - this function returns dynamically allocated array. It has to be released - // afterwards. - DComplex* solve_quartic(double a, double b, double c, double d); + /** + solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e + */ + std::vector<double> solve_quartic_real(long double a, long double b, long double c, + long double d, long double e, + double const epsilon = 1e-12); //! @} -} // namespace corsika::quartic_solver +} // namespace corsika #include <corsika/detail/framework/utility/QuarticSolver.inl> diff --git a/corsika/media/ExponentialRefractiveIndex.hpp b/corsika/media/ExponentialRefractiveIndex.hpp new file mode 100644 index 0000000000000000000000000000000000000000..441f7b33fbb14c7ef9568aa5d942845c7591f0cb --- /dev/null +++ b/corsika/media/ExponentialRefractiveIndex.hpp @@ -0,0 +1,54 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/media/IRefractiveIndexModel.hpp> + +namespace corsika { + + /** + * An exponential refractive index. + * + * This class returns the value of an exponential refractive index + * for all evaluated locations. + * + */ + template <typename T> + class ExponentialRefractiveIndex : public T { + + double n_0; ///< n0 constant. + InverseLengthType lambda_; ///< lambda parameter. + + public: + /** + * Construct an ExponentialRefractiveIndex. + * + * This is initialized with two parameters n_0 and lambda + * and returns the value of the exponential refractive index + * at a given point location. + * + * @param field The refractive index to return to a given point. + */ + template <typename... Args> + ExponentialRefractiveIndex(double const n0, InverseLengthType const lambda, + Args&&... args); + + /** + * Evaluate the refractive index at a given location using its z-coordinate. + * + * @param point The location to evaluate at. + * @returns The refractive index at this point. + */ + double getRefractiveIndex(Point const& point) const final override; + + }; // END: class ExponentialRefractiveIndex + +} // namespace corsika + +#include <corsika/detail/media/ExponentialRefractiveIndex.inl> diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp index 8276c0d1235f44539bfc17468f20aff5eda902b5..613708d1090b19f6cb314f09f62e837ee2e55e20 100644 --- a/corsika/modules/ObservationPlane.hpp +++ b/corsika/modules/ObservationPlane.hpp @@ -13,22 +13,31 @@ #include <corsika/framework/process/ContinuousProcess.hpp> #include <corsika/setup/SetupStack.hpp> #include <corsika/setup/SetupTrajectory.hpp> - -#include <fstream> -#include <string> +#include <corsika/modules/writers/ObservationPlaneWriterParquet.hpp> namespace corsika { /** - * The ObservationPlane writes PDG codes, energies, and distances of particles to the - * central point of the plane into its output file. The particles are considered - * "absorbed" afterwards. + @ingroup Modules + @{ + + The ObservationPlane writes PDG codes, energies, and distances of particles to the + central point of the plane into its output file. The particles are considered + "absorbed" afterwards. + + **Note/Limitation:** as discussed in + https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397 + you cannot put two ObservationPlanes exactly on top of each + other. Even if one of them is "permeable". You have to put a + small gap in between the two plane in such a scenario, or develop + another more specialized output class. */ - class ObservationPlane : public ContinuousProcess<ObservationPlane> { + template <typename TOutputWriter = ObservationPlaneWriterParquet> + class ObservationPlane : public ContinuousProcess<ObservationPlane<TOutputWriter>>, + public TOutputWriter { public: - ObservationPlane(Plane const&, DirectionVector const&, std::string const&, - bool = true); + ObservationPlane(Plane const&, DirectionVector const&, bool = true); ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle, corsika::setup::Trajectory& vTrajectory, @@ -40,16 +49,17 @@ namespace corsika { void showResults() const; void reset(); HEPEnergyType getEnergyGround() const { return energy_ground_; } + YAML::Node getConfig() const; private: Plane const plane_; - std::ofstream outputStream_; bool const deleteOnHit_; HEPEnergyType energy_ground_; unsigned int count_ground_; DirectionVector const xAxis_; DirectionVector const yAxis_; }; + //! @} } // namespace corsika #include <corsika/detail/modules/ObservationPlane.inl> diff --git a/corsika/modules/TrackWriter.hpp b/corsika/modules/TrackWriter.hpp index 5bf9856b389fe80774ea629b27a4b4dc10f33c85..79e1dfc91bbb52d08300699fc3e6ee1e098feb8a 100644 --- a/corsika/modules/TrackWriter.hpp +++ b/corsika/modules/TrackWriter.hpp @@ -8,18 +8,17 @@ #pragma once -#include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/process/ContinuousProcess.hpp> - -#include <fstream> -#include <string> +#include <corsika/modules/writers/TrackWriterParquet.hpp> namespace corsika { - class TrackWriter : public ContinuousProcess<TrackWriter> { + template <typename TOutputWriter = TrackWriterParquet> + class TrackWriter : public ContinuousProcess<TrackWriter<TOutputWriter>>, + public TOutputWriter { public: - TrackWriter(std::string const& filename); + TrackWriter(); template <typename TParticle, typename TTrack> ProcessReturn doContinuous(TParticle const&, TTrack const&, bool const limitFlag); @@ -27,12 +26,7 @@ namespace corsika { template <typename TParticle, typename TTrack> LengthType getMaxStepLength(TParticle const&, TTrack const&); - private: - std::string const filename_; - std::ofstream file_; - - int width_ = 14; - int precision_ = 6; + YAML::Node getConfig() const; }; } // namespace corsika diff --git a/corsika/modules/qgsjetII/Interaction.hpp b/corsika/modules/qgsjetII/Interaction.hpp index a45282f58ae3073237adf19a3c0c739d6153047c..63c878df5d4cea4feadbe2ff028760a79643e741 100644 --- a/corsika/modules/qgsjetII/Interaction.hpp +++ b/corsika/modules/qgsjetII/Interaction.hpp @@ -13,6 +13,10 @@ #include <corsika/framework/random/RNGManager.hpp> #include <corsika/framework/process/InteractionProcess.hpp> #include <corsika/modules/qgsjetII/ParticleConversion.hpp> +#include <corsika/framework/utility/CorsikaData.hpp> + +#include <boost/filesystem/path.hpp> + #include <qgsjet-II-04.hpp> #include <string> @@ -21,14 +25,8 @@ namespace corsika::qgsjetII { class Interaction : public corsika::InteractionProcess<Interaction> { - std::string data_path_; - int count_ = 0; - bool initialized_ = false; - QgsjetIIHadronType alternate_ = - QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles - public: - Interaction(const std::string& dataPath = ""); + Interaction(boost::filesystem::path dataPath = corsika_data("QGSJetII")); ~Interaction(); bool wasInitialized() { return initialized_; } @@ -53,6 +51,11 @@ namespace corsika::qgsjetII { void doInteraction(TSecondaryView&); private: + int count_ = 0; + bool initialized_ = false; + QgsjetIIHadronType alternate_ = + QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles + corsika::default_prng_type& rng_ = corsika::RNGManager::getInstance().getRandomStream("qgsjet"); const int maxMassNumber_ = 208; diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp index 3f6d3e9936f1fd13fcb831359f8faee6768f2948..69cb4a93750b3389749c279cdc0b6f8f84765ad5 100644 --- a/corsika/modules/urqmd/UrQMD.hpp +++ b/corsika/modules/urqmd/UrQMD.hpp @@ -27,7 +27,7 @@ namespace corsika::urqmd { class UrQMD : public InteractionProcess<UrQMD> { public: - UrQMD(boost::filesystem::path const& path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat")); + UrQMD(boost::filesystem::path const path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat")); template <typename TParticle> GrammageType getInteractionLength(TParticle const&) const; @@ -46,7 +46,7 @@ namespace corsika::urqmd { private: static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int); - void readXSFile(boost::filesystem::path const&); + void readXSFile(boost::filesystem::path); // data members default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("urqmd"); diff --git a/corsika/modules/writers/ObservationPlaneWriterParquet.hpp b/corsika/modules/writers/ObservationPlaneWriterParquet.hpp new file mode 100644 index 0000000000000000000000000000000000000000..2a5089c26f129adb6a5281b8474c509096de3ab6 --- /dev/null +++ b/corsika/modules/writers/ObservationPlaneWriterParquet.hpp @@ -0,0 +1,59 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/BaseOutput.hpp> +#include <corsika/output/ParquetStreamer.hpp> +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> + +namespace corsika { + + class ObservationPlaneWriterParquet : public BaseOutput { + + ParquetStreamer output_; ///< The primary output file. + + public: + /** + * Construct an ObservationPlane. + * + * @param name The name of this output. + */ + ObservationPlaneWriterParquet(); + + /** + * Called at the start of each library. + */ + void startOfLibrary(boost::filesystem::path const& directory) final override; + + /** + * Called at the end of each shower. + */ + void endOfShower() final override; + + /** + * Called at the end of each library. + * + * This must also increment the run number since we override + * the default behaviour of BaseOutput. + */ + void endOfLibrary() final override; + + protected: + /** + * Write a particle to the file. + */ + void write(Code const& pid, units::si::HEPEnergyType const& energy, + units::si::LengthType const& x, units::si::LengthType const& y); + + }; // class ObservationPlaneWriterParquet + +} // namespace corsika + +#include <corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl> diff --git a/corsika/modules/writers/TrackWriterParquet.hpp b/corsika/modules/writers/TrackWriterParquet.hpp new file mode 100644 index 0000000000000000000000000000000000000000..7caf90558b2ae0fc5c1cc13ca38b7d009e70b335 --- /dev/null +++ b/corsika/modules/writers/TrackWriterParquet.hpp @@ -0,0 +1,62 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/BaseOutput.hpp> +#include <corsika/output/ParquetStreamer.hpp> +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/geometry/QuantityVector.hpp> + +namespace corsika { + + class TrackWriterParquet : public BaseOutput { + + public: + /** + * Construct a new writer. + * + * @param name The name of this output. + */ + TrackWriterParquet(); + + /** + * Called at the start of each library. + */ + void startOfLibrary(boost::filesystem::path const& directory) final override; + + /** + * Called at the end of each shower. + */ + void endOfShower() final override; + + /** + * Called at the end of each library. + * + * This must also increment the run number since we override + * the default behaviour of BaseOutput. + */ + void endOfLibrary() final override; + + protected: + /** + * Write a track to the file. + */ + void write(Code const& pid, units::si::HEPEnergyType const& energy, + QuantityVector<length_d> const& start, + QuantityVector<length_d> const& end); + + private: + ParquetStreamer output_; ///< The primary output file. + + }; // class TrackWriterParquet + +} // namespace corsika + +#include <corsika/detail/modules/writers/TrackWriterParquet.inl> diff --git a/corsika/output/BaseOutput.hpp b/corsika/output/BaseOutput.hpp new file mode 100644 index 0000000000000000000000000000000000000000..0024c5b7a918e4705fbbadd3bdcdd59dfcab0449 --- /dev/null +++ b/corsika/output/BaseOutput.hpp @@ -0,0 +1,75 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +#include <boost/filesystem.hpp> + +#include <yaml-cpp/yaml.h> + +namespace corsika { + + /** + * This is the base class for all outputs so that they + * can be stored in homogeneous containers. + */ + class BaseOutput { + + protected: + BaseOutput(); + + public: + /** + * Called at the start of each run. + */ + virtual void startOfLibrary(boost::filesystem::path const& directory) = 0; + + /** + * Called at the start of each event/shower. + */ + virtual void startOfShower() {} + + /** + * Called at the end of each event/shower. + */ + virtual void endOfShower() = 0; + + /** + * Called at the end of each run. + */ + virtual void endOfLibrary() = 0; + + /** + * Get the configuration of this output. + */ + virtual YAML::Node getConfig() const = 0; + + /** + * Get any summary information for the entire library. + */ + virtual YAML::Node getSummary() { return YAML::Node(); } + + /** + * Flag to indicate readiness. + */ + bool isInit() const { return is_init_; } + + protected: + /** + * Set init flag. + */ + void setInit(bool const v) { is_init_ = v; } + + int shower_{0}; ///< The current event number. + + private: + bool is_init_{false}; ///< flag to indicate readiness + }; + +} // namespace corsika + +#include <corsika/detail/output/BaseOutput.inl> diff --git a/corsika/output/DummyOutputManager.hpp b/corsika/output/DummyOutputManager.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fc5cf4647d5daa227fceceda5961bb051a840058 --- /dev/null +++ b/corsika/output/DummyOutputManager.hpp @@ -0,0 +1,70 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +namespace corsika { + + /*! + * An output manager that does nothing. + */ + class DummyOutputManager final { + + public: + /** + * Construct an OutputManager instance with a name in a given directory. + * + * @param name The name of this output collection. + * @param dir The directory where the output directory will be stored. + */ + DummyOutputManager(); + + /** + * Handle graceful closure of the outputs upon destruction. + */ + ~DummyOutputManager(); + + /** + * Register an existing output to this manager. + * + * @param name The unique name of this output. + * @param args... These are perfect forwarded to the + * constructor of the output. + */ + template <typename TOutput> + void add(std::string const& name, TOutput& output); + + /** + * Called at the start of each library. + * + * This iteratively calls startOfLibrary on each registered output. + */ + void startOfLibrary(); + + /** + * Called at the start of each event/shower. + * This iteratively calls startOfEvent on each registered output. + */ + void startOfShower(); + + /** + * Called at the end of each event/shower. + * This iteratively calls endOfEvent on each registered output. + */ + void endOfShower(); + + /** + * Called at the end of each library. + * This iteratively calls endOfLibrary on each registered output. + */ + void endOfLibrary(); + + }; // class DummyOutputManager + +} // namespace corsika + +#include <corsika/detail/output/DummyOutputManager.inl> diff --git a/corsika/output/NoOutput.hpp b/corsika/output/NoOutput.hpp new file mode 100644 index 0000000000000000000000000000000000000000..289c19ad720f9f03d621d3b50ac031c40cc794cf --- /dev/null +++ b/corsika/output/NoOutput.hpp @@ -0,0 +1,62 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> + +namespace corsika { + + /** + * This is the base class for all outputs so that they + * can be stored in homogeneous containers. + */ + class NoOutput : public BaseOutput { + + protected: + NoOutput() {} + + public: + /** + * Called at the start of each run. + */ + void startOfLibrary(boost::filesystem::path const&) final override {} + + /** + * Called at the start of each event/shower. + */ + void startOfShower() final override {} + + /** + * Called at the end of each event/shower. + */ + void endOfShower() final override {} + + /** + * Called at the end of each run. + */ + void endOfLibrary() final override {} + + /** + * Get the configuration of this output. + */ + YAML::Node getConfig() const { return YAML::Node(); }; + + /** + * Get any summary information for the entire library. + */ + YAML::Node getSummary() final override { return YAML::Node(); }; + + protected: + void write(Code const&, units::si::HEPEnergyType const&, units::si::LengthType const&, + units::si::LengthType const&) {} + }; + +} // namespace corsika + +#include <corsika/detail/output/BaseOutput.inl> diff --git a/corsika/output/OutputManager.hpp b/corsika/output/OutputManager.hpp new file mode 100644 index 0000000000000000000000000000000000000000..2582220520454d028b660ca06f5ac93ed8a18b73 --- /dev/null +++ b/corsika/output/OutputManager.hpp @@ -0,0 +1,118 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +#include <chrono> +#include <string> +#include <boost/filesystem.hpp> +#include <corsika/output/BaseOutput.hpp> +#include <corsika/framework/core/Logging.hpp> + +namespace corsika { + + /*! + * Manages CORSIKA 8 output streams. + */ + class OutputManager final { + + /** + * Indicates the current state of this manager. + */ + enum class OutputState { + NoInit, + LibraryReady, + ShowerInProgress, + LibraryFinished, + }; + + OutputState state_{OutputState::NoInit}; ///< The current state of this manager. + std::string const name_; ///< The name of this simulation file. + boost::filesystem::path const root_; ///< The top-level directory for the output. + int count_{0}; ///< The current ID of this shower. + std::chrono::time_point<std::chrono::system_clock> const start_time{ + std::chrono::system_clock::now()}; ///< The time the manager is created. + inline static auto logger{get_logger("output")}; ///< A custom logger. + /** + * The outputs that have been registered with us. + */ + std::map<std::string, std::reference_wrapper<BaseOutput>> outputs_; + + /** + * Write a YAML-node to a file. + */ + void writeNode(YAML::Node const& node, boost::filesystem::path const& path) const; + + /** + * Write the top-level config of this simulation. + */ + void writeTopLevelConfig() const; + + /** + * Initialize the "registered" output with a given name. + */ + void initOutput(std::string const& name) const; + + /** + * Write the top-level summary of this library. + */ + void writeTopLevelSummary() const; + + public: + /** + * Construct an OutputManager instance with a name in a given directory. + * + * @param name The name of this output collection. + * @param dir The directory where the output directory will be stored. + */ + OutputManager(std::string const& name, boost::filesystem::path const& dir); + + /** + * Handle graceful closure of the outputs upon destruction. + */ + ~OutputManager(); + + /** + * Register an existing output to this manager. + * + * @param name The unique name of this output. + * @param args... These are perfect forwarded to the + * constructor of the output. + */ + template <typename TOutput> + void add(std::string const& name, TOutput& output); + + /** + * Called at the start of each library. + * + * This iteratively calls startOfLibrary on each registered output. + */ + void startOfLibrary(); + + /** + * Called at the start of each event/shower. + * This iteratively calls startOfEvent on each registered output. + */ + void startOfShower(); + + /** + * Called at the end of each event/shower. + * This iteratively calls endOfEvent on each registered output. + */ + void endOfShower(); + + /** + * Called at the end of each library. + * This iteratively calls endOfLibrary on each registered output. + */ + void endOfLibrary(); + + }; // class OutputManager + +} // namespace corsika + +#include <corsika/detail/output/OutputManager.inl> diff --git a/corsika/output/ParquetStreamer.hpp b/corsika/output/ParquetStreamer.hpp new file mode 100644 index 0000000000000000000000000000000000000000..dcf1f547f72338b45459583a591150487de937c8 --- /dev/null +++ b/corsika/output/ParquetStreamer.hpp @@ -0,0 +1,83 @@ + +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <string> + +// clang-format-off, NOTE: the order of these includes is *important* +// you will get unhelpful compiler errors about unknown +// operator definitions if these are reordered +#include <parquet/stream_writer.h> +#include <parquet/arrow/schema.h> +#include <arrow/io/file.h> +// clang-format-on + +namespace corsika { + + /** + * This class automates the construction of simple tabular + * Parquet files using the parquet::StreamWriter. + */ + class ParquetStreamer { + + public: + /** + * ParquetStreamer's take no constructor arguments. + */ + ParquetStreamer(); + + /** + * Initialize the streamer to write to a given file. + */ + void initStreamer(std::string const& filepath); + + /** + * Add a field to this streamer. + */ + template <typename... TArgs> + void addField(TArgs&&... args); + + /** + * Enable compression for this streamer. + */ + void enableCompression(int const level = 3); + + /** + * Finalize the streamer construction. + */ + void buildStreamer(); + + /** + * Finish writing this stream. + */ + void closeStreamer(); + + /** + * Return a reference to the underlying writer. + */ + std::shared_ptr<parquet::StreamWriter> getWriter(); + + /** + * @return status of streamer + */ + bool isInit() const { return isInit_; } + + private: + bool isInit_ = false; ///< flag to handle state of writer + parquet::WriterProperties::Builder builder_; ///< The writer properties builder. + parquet::schema::NodeVector fields_; ///< The fields in this file. + std::shared_ptr<parquet::schema::GroupNode> schema_; ///< The schema for this file. + std::shared_ptr<arrow::io::FileOutputStream> outfile_; ///< The output file. + std::shared_ptr<parquet::StreamWriter> writer_; ///< The stream writer to 'outfile' + + }; // class ParquetStreamer +} // namespace corsika + +#include <corsika/detail/output/ParquetStreamer.inl> diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp index 34bd1377d99a5167b2ac406e506a4c8492173bc3..19e9a44c32aa401147d57585c89a7c1ab3ec6f4e 100644 --- a/corsika/stack/NuclearStackExtension.hpp +++ b/corsika/stack/NuclearStackExtension.hpp @@ -92,6 +92,10 @@ namespace corsika::nuclear_stack { * Overwrite normal getParticleMass function with nuclear version */ HEPMassType getMass() const; + /** + * Overwrite normal getParticleCharge function with nuclear version + */ + ElectricChargeType getCharge() const; /** * Overwirte normal getChargeNumber function with nuclear version **/ diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp index adaae996b63a5b69efa85dc8ad8956e253de44b9..e5a765781ae44f38b38da7d4135a6363c697f5d3 100644 --- a/corsika/stack/VectorStack.hpp +++ b/corsika/stack/VectorStack.hpp @@ -87,8 +87,14 @@ namespace corsika { return this->getMomentum() / this->getEnergy(); } + VelocityVector getVelocity() const { + return this->getMomentum() / this->getEnergy() * constants::c; + } + HEPMassType getMass() const { return get_mass(this->getPID()); } + ElectricChargeType getCharge() const { return get_charge(this->getPID()); } + int16_t getChargeNumber() const { return get_charge_number(this->getPID()); } ///@} }; diff --git a/do-clang-format.py b/do-clang-format.py index d6ea597535e69da32de80ff31962ba998523305d..fcfe674116832aea4ab94eb0c930b5d5c512fced 100755 --- a/do-clang-format.py +++ b/do-clang-format.py @@ -10,6 +10,7 @@ import os import sys import re +debug = False do_progress = False try: from progress.bar import ChargingBar @@ -77,6 +78,9 @@ else: filelist_clean.append(f) filelist = filelist_clean +if debug: + print ("filelist: ", filelist) + cmd = "clang-format" if "CLANG_FORMAT" in os.environ: cmd = os.environ["CLANG_FORMAT"] @@ -99,11 +103,18 @@ if do_progress: bar = ChargingBar('Processing', max=len(filelist)) if args.apply: + changed = [] for filename in filelist: if bar: bar.next() + a = open(filename, "rb").read() subp.check_call(cmd.split() + ["-i", filename]) + b = open(filename, "rb").read() + if a != b: + changed.append(filename) if bar: bar.finish() - + if debug: + print ("changed: ", changed) + else: # only print files which need formatting files_need_formatting = 0 diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 8cd42babe67261c0f17e52805e67bfceb8a09ab7..3240a632cb15c390454e99e0515c52f7945336e7 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -7,7 +7,7 @@ set (CMAKE_VERBOSE_MAKEFILE OFF) # this can be changed with `make VERBOSE=1` find_package (corsika CONFIG REQUIRED) # this is only for CORSIKA_REGISTER_EXAMPLE -include ("${CMAKE_CURRENT_SOURCE_DIR}/CMakeHelper.cmake") +include ("${CMAKE_CURRENT_SOURCE_DIR}/corsikaExamples.cmake") add_executable (helix_example helix_example.cpp) target_link_libraries (helix_example CORSIKA8::CORSIKA8) diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp index bca45481a3841f690e80ea5d936063a4ffe91acf..ebb202597bc8f7d09e88b815e89e5f55cf5a3225 100644 --- a/examples/boundary_example.cpp +++ b/examples/boundary_example.cpp @@ -14,6 +14,8 @@ #include <corsika/framework/utility/CorsikaFenv.hpp> #include <corsika/framework/core/Logging.hpp> +#include <corsika/output/OutputManager.hpp> + #include <corsika/setup/SetupEnvironment.hpp> #include <corsika/setup/SetupStack.hpp> #include <corsika/setup/SetupTrajectory.hpp> @@ -42,6 +44,7 @@ #include <iostream> #include <limits> #include <typeinfo> +#include <fstream> using namespace corsika; using namespace std; @@ -112,6 +115,8 @@ int main() { world->addChild(std::move(target)); universe.addChild(std::move(world)); + OutputManager output("boundary_outputs"); + // setup processes, decays and interactions setup::Tracking tracking; @@ -121,7 +126,9 @@ int main() { ParticleCut cut(50_GeV, true, true); - TrackWriter trackWriter("boundary_tracks.dat"); + TrackWriter trackWriter; + output.add("tracks", trackWriter); // register TrackWriter + MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat"); // assemble all processes into an ordered process list @@ -164,7 +171,7 @@ int main() { } // define air shower object, run simulation - Cascade EAS(env, tracking, sequence, stack); + Cascade EAS(env, tracking, sequence, output, stack); EAS.run(); @@ -174,4 +181,6 @@ int main() { (cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy()); CORSIKA_LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV, (Efinal / E0 - 1.) * 100); + + output.endOfLibrary(); } diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp index eb0938d9fe08e93ac3a3f9d972362d9c58a0ef48..f7758db0a3163521274ee5b0b33e21cc478c4405 100644 --- a/examples/cascade_example.cpp +++ b/examples/cascade_example.cpp @@ -11,10 +11,11 @@ #include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/random/RNGManager.hpp> #include <corsika/framework/geometry/Sphere.hpp> - #include <corsika/framework/utility/CorsikaFenv.hpp> #include <corsika/framework/core/Logging.hpp> +#include <corsika/output/OutputManager.hpp> + #include <corsika/media/Environment.hpp> #include <corsika/media/HomogeneousMedium.hpp> #include <corsika/media/NuclearComposition.hpp> @@ -54,8 +55,7 @@ using namespace std; // int main() { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); - logging::set_level(logging::level::trace); + logging::set_level(logging::level::info); std::cout << "cascade_example" << std::endl; @@ -107,6 +107,8 @@ int main() { rootCS, 0_m, 0_m, height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe + OutputManager output("cascade_outputs"); + ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env}; { @@ -140,7 +142,9 @@ int main() { // cascade with only HE model ==> HE cut ParticleCut cut(80_GeV, true, true); - TrackWriter trackWriter("tracks.dat"); + TrackWriter trackWriter; + output.add("tracks", trackWriter); // register TrackWriter + BetheBlochPDG eLoss{showerAxis}; // assemble all processes into an ordered process list @@ -148,7 +152,7 @@ int main() { make_sequence(stackInspect, sibyll, sibyllNuc, decay, eLoss, cut, trackWriter); // define air shower object, run simulation - Cascade EAS(env, tracking, sequence, stack); + Cascade EAS(env, tracking, sequence, output, stack); EAS.run(); @@ -162,4 +166,6 @@ int main() { cout << "total dEdX energy (GeV): " << eLoss.getTotal() / 1_GeV << endl << "relative difference (%): " << eLoss.getTotal() / E0 * 100 << endl; cut.reset(); + + output.endOfLibrary(); } diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp index 4b24f3b0148045a2db97546c37c6a2916e3c1d1e..a02c914e4ad1d554160d36745e5d00a2ee6d981e 100644 --- a/examples/cascade_proton_example.cpp +++ b/examples/cascade_proton_example.cpp @@ -11,10 +11,11 @@ #include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/random/RNGManager.hpp> #include <corsika/framework/geometry/Sphere.hpp> - #include <corsika/framework/utility/CorsikaFenv.hpp> #include <corsika/framework/core/Logging.hpp> +#include <corsika/output/OutputManager.hpp> + #include <corsika/media/Environment.hpp> #include <corsika/media/HomogeneousMedium.hpp> #include <corsika/media/NuclearComposition.hpp> @@ -66,6 +67,8 @@ int main() { // initialize random number sequence(s) RNGManager::getInstance().registerRandomStream("cascade"); + OutputManager output("cascade_proton_outputs"); + // setup environment, geometry using EnvType = setup::Environment; EnvType env; @@ -115,7 +118,7 @@ int main() { // setup processes, decays and interactions setup::Tracking tracking; - StackInspector<setup::Stack> stackInspect(1, true, E0); + StackInspector<setup::Stack> stackInspect(1000, true, E0); RNGManager::getInstance().registerRandomStream("sibyll"); RNGManager::getInstance().registerRandomStream("pythia"); @@ -130,7 +133,9 @@ int main() { // HadronicElasticModel::HadronicElasticInteraction // hadronicElastic(env); - TrackWriter trackWriter("tracks.dat"); + TrackWriter trackWriter; + output.add("tracks", trackWriter); // register TrackWriter + ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env}; BetheBlochPDG eLoss{showerAxis}; @@ -141,7 +146,7 @@ int main() { auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect); // define air shower object, run simulation - Cascade EAS(env, tracking, sequence, stack); + Cascade EAS(env, tracking, sequence, output, stack); EAS.run(); cout << "Result: E0=" << E0 / 1_GeV << endl; @@ -150,4 +155,6 @@ int main() { cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy(); cout << "total energy (GeV): " << Efinal / 1_GeV << endl << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl; + + output.endOfLibrary(); } diff --git a/examples/CMakeHelper.cmake b/examples/corsikaExamples.cmake similarity index 100% rename from examples/CMakeHelper.cmake rename to examples/corsikaExamples.cmake diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp index d4d24a8794cb54de46019bac71f7aad23e6657e0..606d72b41cf10445d62bb3d769e1465e90a2068b 100644 --- a/examples/em_shower.cpp +++ b/examples/em_shower.cpp @@ -18,6 +18,8 @@ #include <corsika/framework/process/InteractionCounter.hpp> #include <corsika/framework/core/Logging.hpp> +#include <corsika/output/OutputManager.hpp> + #include <corsika/media/Environment.hpp> #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> #include <corsika/media/NuclearComposition.hpp> @@ -66,7 +68,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; int main(int argc, char** argv) { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); if (argc != 2) { std::cerr << "usage: em_shower <energy/GeV>" << std::endl; @@ -137,6 +138,7 @@ int main(int argc, char** argv) { std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02 << std::endl; + OutputManager output("em_shower_outputs"); ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env}; // setup processes, decays and interactions @@ -146,7 +148,8 @@ int main(int argc, char** argv) { corsika::proposal::ContinuousProcess emContinuous(env); InteractionCounter emCascadeCounted(emCascade); - TrackWriter trackWriter("tracks.dat"); + TrackWriter trackWriter; + output.add("tracks", trackWriter); // register TrackWriter // long. profile; columns for photon, e+, e- still need to be added LongitudinalProfile longprof{showerAxis}; @@ -154,12 +157,13 @@ int main(int argc, char** argv) { Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}), "particles.dat"); + output.add("obsplane", observationLevel); auto sequence = make_sequence(emCascadeCounted, emContinuous, longprof, cut, observationLevel, trackWriter); // define air shower object, run simulation setup::Tracking tracking; - Cascade EAS(env, tracking, sequence, stack); + Cascade EAS(env, tracking, sequence, output, stack); // to fix the point of first interaction, uncomment the following two lines: // EAS.setNodes(); @@ -183,4 +187,6 @@ int main(int argc, char** argv) { save_hist(hists.labHist(), "inthist_lab_emShower.npz", true); save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true); longprof.save("longprof_emShower.txt"); + + output.endOfLibrary(); } diff --git a/examples/geometry_example.cpp b/examples/geometry_example.cpp index 7a25708a8d434b8b35e89a8cd23b1e4ad35e0f1f..596e6c1842011d56d8db27dc6e02fe81ab8bdbb9 100644 --- a/examples/geometry_example.cpp +++ b/examples/geometry_example.cpp @@ -21,7 +21,6 @@ using namespace corsika; int main() { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); CORSIKA_LOG_INFO("geometry_example"); diff --git a/examples/helix_example.cpp b/examples/helix_example.cpp index cf58531b4170aa5b74a3693033b71c700cde102e..e980931b74ace074bf6d0459cf69e6d158f32b91 100644 --- a/examples/helix_example.cpp +++ b/examples/helix_example.cpp @@ -21,7 +21,6 @@ using namespace corsika; int main() { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); CORSIKA_LOG_INFO("helix_example"); diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp index af379ca89bc2a14ccc706e653ca6e3fa25e7d49c..f1602f5ce24a7359d8e1de500369f70d32177b5f 100644 --- a/examples/hybrid_MC.cpp +++ b/examples/hybrid_MC.cpp @@ -25,6 +25,8 @@ #include <corsika/framework/core/Cascade.hpp> #include <corsika/framework/geometry/PhysicalGeometry.hpp> +#include <corsika/output/OutputManager.hpp> + #include <corsika/media/Environment.hpp> #include <corsika/media/FlatExponential.hpp> #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> @@ -86,7 +88,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; int main(int argc, char** argv) { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); CORSIKA_LOG_INFO("hybrid_MC"); @@ -171,6 +172,7 @@ int main(int argc, char** argv) { std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02 << std::endl; + OutputManager output("hybrid_MC_outputs"); ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env}; // setup processes, decays and interactions @@ -219,6 +221,7 @@ int main(int argc, char** argv) { Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}), "particles.dat"); + output.add("obsplane", observationLevel); corsika::urqmd::UrQMD urqmd_model; InteractionCounter urqmdCounted{urqmd_model}; @@ -240,7 +243,7 @@ int main(int argc, char** argv) { // define air shower object, run simulation setup::Tracking tracking; - Cascade EAS(env, tracking, sequence, stack); + Cascade EAS(env, tracking, sequence, output, stack); // to fix the point of first interaction, uncomment the following two lines: // EAS.SetNodes(); @@ -272,4 +275,6 @@ int main(int argc, char** argv) { std::ofstream finish("finished"); finish << "run completed without error" << std::endl; + + output.endOfLibrary(); } diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp index e5f6455d3f6782cef10d69c128e38534c50a2d8d..928197970e8c27f08b93d561f33aab1defea5796 100644 --- a/examples/vertical_EAS.cpp +++ b/examples/vertical_EAS.cpp @@ -25,6 +25,8 @@ #include <corsika/framework/core/Cascade.hpp> #include <corsika/framework/geometry/PhysicalGeometry.hpp> +#include <corsika/output/OutputManager.hpp> + #include <corsika/media/Environment.hpp> #include <corsika/media/FlatExponential.hpp> #include <corsika/media/HomogeneousMedium.hpp> @@ -97,7 +99,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; int main(int argc, char** argv) { - corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v"); logging::set_level(logging::level::info); CORSIKA_LOG_INFO("vertical_EAS"); @@ -134,6 +135,7 @@ int main(int argc, char** argv) { {{Code::Nitrogen, Code::Oxygen}, {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now + builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km); builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km); builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km); builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km); @@ -141,22 +143,6 @@ int main(int argc, char** argv) { builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean); builder.assemble(env); - CORSIKA_LOG_DEBUG( - "environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, " - "layer5={}", - fmt::ptr(env.getUniverse()->getContainingNode( - Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))), - fmt::ptr(env.getUniverse()->getContainingNode( - Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))), - fmt::ptr(env.getUniverse()->getContainingNode( - Point(rootCS, {constants::EarthRadius::Mean + 50_km, 0_m, 0_m}))), - fmt::ptr(env.getUniverse()->getContainingNode( - Point(rootCS, {constants::EarthRadius::Mean + 20_km, 0_m, 0_m}))), - fmt::ptr(env.getUniverse()->getContainingNode( - Point(rootCS, {constants::EarthRadius::Mean + 5_km, 0_m, 0_m}))), - fmt::ptr(env.getUniverse()->getContainingNode( - Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m})))); - // pre-setup particle stack unsigned short const A = std::stoi(std::string(argv[1])); Code beamCode; @@ -214,6 +200,9 @@ int main(int argc, char** argv) { ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env}; + // create the output manager that we then register outputs with + OutputManager output("vertical_EAS_outputs"); + // setup processes, decays and interactions // corsika::qgsjetII::Interaction qgsjet; @@ -280,8 +269,6 @@ int main(int argc, char** argv) { string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz"; string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz"; string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt"; - string const tracks_file = "tracks_" + to_string(i_shower) + ".dat"; - string const particles_file = "particles_" + to_string(i_shower) + ".dat"; std::cout << std::endl; std::cout << "Shower " << i_shower << "/" << number_showers << std::endl; @@ -289,6 +276,46 @@ int main(int argc, char** argv) { // setup particle stack, and add primary particle setup::Stack stack; stack.clear(); + unsigned short const A = std::stoi(std::string(argv[1])); + Code beamCode; + HEPEnergyType mass; + unsigned short Z = 0; + if (A > 0) { + beamCode = Code::Nucleus; + Z = std::stoi(std::string(argv[2])); + mass = get_nucleus_mass(A, Z); + } else { + int pdg = std::stoi(std::string(argv[2])); + beamCode = convert_from_PDG(PDGCode(pdg)); + mass = get_mass(beamCode); + } + HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3])); + double theta = 0.; + auto const thetaRad = theta / 180. * M_PI; + + auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { + return sqrt((Elab - m) * (Elab + m)); + }; + HEPMomentumType P0 = elab2plab(E0, mass); + auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) { + return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad)); + }; + + auto const [px, py, pz] = momentumComponents(thetaRad, P0); + auto plab = MomentumVector(rootCS, {px, py, pz}); + cout << "input particle: " << beamCode << endl; + cout << "input angles: theta=" << theta << endl; + cout << "input momentum: " << plab.getComponents() / 1_GeV + << ", norm = " << plab.getNorm() << endl; + + auto const observationHeight = 0_km + builder.getEarthRadius(); + auto const injectionHeight = 111.75_km + builder.getEarthRadius(); + auto const t = (injectionHeight - observationHeight) / cos(thetaRad); + Point const showerCore{rootCS, 0_m, 0_m, observationHeight}; + Point const injectionPos = + showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t; + + std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl; if (A > 1) { stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z)); @@ -308,9 +335,64 @@ int main(int argc, char** argv) { } } - TrackWriter trackWriter(tracks_file); - ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}), - particles_file); + // we make the axis much longer than the inj-core distance since the + // profile will go beyond the core, depending on zenith angle + std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5 + << std::endl; + + ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env, + false}; + + // setup processes, decays and interactions + + // corsika::qgsjetII::Interaction qgsjet; + corsika::sibyll::Interaction sibyll; + InteractionCounter sibyllCounted(sibyll); + + corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env); + InteractionCounter sibyllNucCounted(sibyllNuc); + + corsika::pythia8::Decay decayPythia; + + // use sibyll decay routine for decays of particles unknown to pythia + corsika::sibyll::Decay decaySibyll{{ + Code::N1440Plus, + Code::N1440MinusBar, + Code::N1440_0, + Code::N1440_0Bar, + Code::N1710Plus, + Code::N1710MinusBar, + Code::N1710_0, + Code::N1710_0Bar, + + Code::Pi1300Plus, + Code::Pi1300Minus, + Code::Pi1300_0, + + Code::KStar0_1430_0, + Code::KStar0_1430_0Bar, + Code::KStar0_1430_Plus, + Code::KStar0_1430_MinusBar, + }}; + + decaySibyll.printDecayConfig(); + + ParticleCut cut{60_GeV, true, true}; + // corsika::proposal::Interaction emCascade(env); + // corsika::proposal::ContinuousProcess emContinuous(env); + // InteractionCounter emCascadeCounted(emCascade); + BetheBlochPDG emContinuous(showerAxis); + + OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false); + TrackWriter trackWriter; + output.add("tracks", trackWriter); // register TrackWriter + + LongitudinalProfile longprof{showerAxis}; + + Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); + ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.})); + // register the observation plane with the output + output.add("obsplane", observationLevel); auto sequence = make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence, @@ -318,7 +400,7 @@ int main(int argc, char** argv) { // define air shower object, run simulation setup::Tracking tracking; - Cascade EAS(env, tracking, sequence, stack); + Cascade EAS(env, tracking, sequence, output, stack); // to fix the point of first interaction, uncomment the following two lines: // EAS.forceInteraction(); @@ -326,16 +408,16 @@ int main(int argc, char** argv) { EAS.run(); cut.showResults(); - emContinuous.showResults(); + // emContinuous.showResults(); observationLevel.showResults(); const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() + - cut.getEmEnergy() + emContinuous.getEnergyLost() + + cut.getEmEnergy() + // emContinuous.getEnergyLost() + observationLevel.getEnergyGround(); cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; observationLevel.reset(); cut.reset(); - emContinuous.reset(); + // emContinuous.reset(); auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() + urqmdCounted.getHistogram(); @@ -343,5 +425,7 @@ int main(int argc, char** argv) { save_hist(hists.labHist(), labHist_file, true); save_hist(hists.CMSHist(), cMSHist_file, true); longprof.save(longprof_file); + + output.endOfLibrary(); } } diff --git a/modules/qgsjetII/qgsjet-II-04.cpp b/modules/qgsjetII/qgsjet-II-04.cpp index eea2724a88c07f1b028aba0703173a1ebc1f0504..8ffd24df4bd7b860116b7792e16671562af051b6 100644 --- a/modules/qgsjetII/qgsjet-II-04.cpp +++ b/modules/qgsjetII/qgsjet-II-04.cpp @@ -2,7 +2,7 @@ #include <iostream> -datadir::datadir(const std::string& dir) { +datadir::datadir(std::string const& dir) { if (dir.length() > 130) { std::cerr << "QGSJetII error, will cut datadir \"" << dir << "\" to 130 characters: " << std::endl; @@ -13,11 +13,9 @@ datadir::datadir(const std::string& dir) { data[i + 1] = '\0'; } - /** @function qgran link to random number generation */ -double qgran_(int&) { return qgsjetII::rndm_interface(); } - +double qgran_(int&) { return qgsjetII::rndm_interface(); } diff --git a/python/.gitignore b/python/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..7b9e44cc96089e3d033d74eb46965d968bd28d52 --- /dev/null +++ b/python/.gitignore @@ -0,0 +1,147 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# ignore any generated output files +*.npz +*.dat + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +pip-wheel-metadata/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# static files generated from Django application using `collectstatic` +media +static diff --git a/python/Makefile b/python/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..ccc13257fe572b456857657929c41d8199ad6529 --- /dev/null +++ b/python/Makefile @@ -0,0 +1,30 @@ +## +# ## +# corsika +# +# @file + +# find python3 +PYTHON=`which python3` + +# our testing targets +.PHONY: tests flake black isort all + +all: mypy isort black flake tests + +tests: + ${PYTHON} -m pytest --cov=corsika tests + +flake: + ${PYTHON} -m flake8 corsika tests + +black: + ${PYTHON} -m black -t py37 corsika tests + +isort: + ${PYTHON} -m isort --atomic corsika tests + +mypy: + ${PYTHON} -m mypy corsika + +# end diff --git a/python/README.md b/python/README.md new file mode 100644 index 0000000000000000000000000000000000000000..7e35dfe389dde4d372b1e38301c41c42a45c6fa8 --- /dev/null +++ b/python/README.md @@ -0,0 +1,8 @@ +# CORSIKA 8 - Python Library + +To install this into your global environment using `pip` (not recommended), run + +``` shell +pip install --user git+https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/tree/master/python +``` + diff --git a/python/corsika/__init__.py b/python/corsika/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..078d688098095f06e2a14d2db4c203370993d578 --- /dev/null +++ b/python/corsika/__init__.py @@ -0,0 +1,17 @@ +""" + A Python interface to CORSIKA 8. + + (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" + +from . import io +from .io.library import Library + +# all imported objects +__all__ = ["io", "Library"] + +__version__: str = "8.0.0-alpha" diff --git a/python/corsika/io/__init__.py b/python/corsika/io/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..296572d6dee46a4b95a61b8ef18c549f120a1205 --- /dev/null +++ b/python/corsika/io/__init__.py @@ -0,0 +1,15 @@ +""" + The 'io' module provides for reading CORSIKA8 output files. + + (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" + +from .hist import read_hist +from .library import Library + +# all exported objects +__all__ = ["read_hist", "Library"] diff --git a/python/corsika/io/hist.py b/python/corsika/io/hist.py new file mode 100644 index 0000000000000000000000000000000000000000..615db0faec0ec63cff8fe00fc5ebe1d8ce08f2ff --- /dev/null +++ b/python/corsika/io/hist.py @@ -0,0 +1,70 @@ +""" + This file supports reading boost_histograms from CORSIKA8. + + (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" + +import boost_histogram as bh +import numpy as np + + +def read_hist(filename: str) -> bh.Histogram: + """ + Read a histogram produced with CORSIKA8's `SaveBoostHistogram()` function. + + Parameters + ---------- + filename: str + The filename of the .npy file containing the histogram. + + Returns + ------- + hist: bh.Histogram + An initialized bh.Histogram instance. + + Throws + ------ + ValueError: + If the histogram type is not supported. + """ + + # load the filenames + d = np.load(filename) + + # extract the axis and overflows information + axistypes = d["axistypes"].view("c") + overflow = d["overflow"] + underflow = d["underflow"] + + # this is where we store the axes that we extract from the file. + axes = [] + + # we now loop over the axes + for i, (at, has_overflow, has_underflow) in enumerate( + zip(axistypes, overflow, underflow) + ): + + # continuous histogram + if at == b"c": + axes.append( + bh.axis.Variable( + d[f"binedges_{i}"], overflow=has_overflow, underflow=has_underflow + ) + ) + # discrete histogram + elif at == b"d": + axes.append(bh.axis.IntCategory(d[f"bins_{i}"], growth=(not has_overflow))) + + # and unknown histogram type + else: + raise ValueError(f"'{at}' is not a valid C8 histogram axistype.") + + # create the histogram and fill it in + h = bh.Histogram(*axes) + h.view(flow=True)[:] = d["data"] + + return h diff --git a/python/corsika/io/library.py b/python/corsika/io/library.py new file mode 100644 index 0000000000000000000000000000000000000000..236fc97ac972a11c0fb366876bd92cae20d6ab18 --- /dev/null +++ b/python/corsika/io/library.py @@ -0,0 +1,231 @@ +""" + This file allows for reading/working with C8 libraries. + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import logging +import os +import os.path as op +import re +from typing import Any, Dict, Optional, List + +import yaml + +from . import outputs + + +class Library(object): + """ + Represents a library ("run") of showers produced by C8. + """ + + def __init__(self, path: str): + """ + + Parameters + ---------- + path: str + The path to the directory containing the library. + + + Raises + ------ + ValueError + If `path` does not contain a valid CORSIKA8 library. + """ + + # check that this is a valid library + if not self.__valid_library(path): + raise ValueError(f"'{path}' does not contain a valid CORSIKA8 library.") + + # store the top-level path + self.path = path + + # load the config and summary files + self.config = self.load_config(path) + self.summary = self.load_summary(path) + + # build the list of outputs + self.__outputs = self.__build_outputs(path) + + @property + def names(self) -> List[str]: + """ + Return the list of registered outputs. + """ + return list(self.__outputs.keys()) + + @property + def modules(self) -> Dict[str, str]: + """ + Return the list of registered outputs. + """ + pass + + def get(self, name: str) -> Optional[outputs.Output]: + """ + Return the output with a given name. + """ + if name in self.__outputs: + return self.__outputs[name] + else: + msg = f"Output with name '{name}' not available in this library." + logging.getLogger("corsika").warn(msg) + return None + + @staticmethod + def load_config(path: str) -> Dict[str, Any]: + """ + Load the top-level config from a given library path. + + + Parameters + ---------- + path: str + The path to the directory containing the library. + + Returns + ------- + dict: + The config as a python dictionary. + + Raises + ------ + FileNotFoundError + If the config file cannot be found + + """ + with open(op.join(path, "config.yaml"), "r") as f: + return yaml.load(f, Loader=yaml.Loader) + + @staticmethod + def load_summary(path: str) -> Dict[str, Any]: + """ + Load the top-level summary from a given library path. + + + Parameters + ---------- + path: str + The path to the directory containing the library. + + Returns + ------- + dict: + The summary as a python dictionary. + + Raises + ------ + FileNotFoundError + If the summary file cannot be found + + """ + with open(op.join(path, "summary.yaml"), "r") as f: + return yaml.load(f, Loader=yaml.Loader) + + @staticmethod + def __valid_library(path: str) -> bool: + """ + Check if the library pointed to by 'path' is a valid C8 library. + + Parameters + ---------- + path: str + The path to the directory containing the library. + + Returns + ------- + bool: + True if this is a valid C8 library. + + """ + + # check that the config file exists + if not op.exists(op.join(path, "config.yaml")): + return False + + # the config file exists, we load it + config = Library.load_config(path) + + # and check that the config's "writer" key is correct + return config["creator"] == "CORSIKA8" + + @staticmethod + def __build_outputs(path: str) -> Dict[str, outputs.Output]: + """ + Build the outputs contained in this library. + + This will print a warning message if a particular + output is invalid but will continue to load additional + outputs afterwards. + + Parameters + ---------- + path: str + The path to the directory containing this library. + + Returns + ------- + Dict[str, Output]: + A dictionary mapping names to initialized outputs. + """ + + # get a list of the subdirectories in the library + _, dirs, _ = next(os.walk(path)) + + # this is the dictionary where we store our components + components: Dict[str, Any] = {} + + # loop over the subdirectories + for subdir in dirs: + + # read the config file for this output + config = Library.load_config(op.join(path, subdir)) + + # the name keyword is our unique identifier + name = config.get("name") + + # get the "type" keyword to identify this component + out_type = config.get("type") + + # if `out_type` was None, this is an invalid output + if out_type is None or name is None: + msg = ( + f"'{subdir}' does not contain a valid config." + "Missing 'type' or 'name' keyword." + ) + logging.getLogger("corsika").warn(msg) + continue # skip to the next output, don't error + + # we now have a valid component type, get the corresponding + # type from the proccesses subdirectory + try: + + # instantiate the output and store it in our dict + component = getattr(outputs, out_type)(op.join(path, subdir)) + + # check if the read failed + if not component.is_good(): + msg = ( + f"'{name}' encountered an error while reading. " + "This process will be not be loaded." + ) + logging.getLogger("corsika").warn(msg) + else: + components[name] = component + + except AttributeError as e: + msg = ( + f"Unable to instantiate an instance of '{out_type}' " + f"for a process called '{name}'" + ) + logging.getLogger("corsika").warn(msg) + logging.getLogger("corsika").warn(e) + continue # skip to the next output, don't error + + # and we are done building - return the constructed outputs + return components diff --git a/python/corsika/io/outputs/__init__.py b/python/corsika/io/outputs/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..fa38acdde3c2dc8dd3cc351c16d27159cecd1efb --- /dev/null +++ b/python/corsika/io/outputs/__init__.py @@ -0,0 +1,14 @@ +""" + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" + +from .observation_plane import ObservationPlane +from .track_writer import TrackWriter +from .output import Output + +__all__ = ["Output", "ObservationPlane", "TrackWriter"] diff --git a/python/corsika/io/outputs/observation_plane.py b/python/corsika/io/outputs/observation_plane.py new file mode 100644 index 0000000000000000000000000000000000000000..271d145b4f076f1c7d747d6cc5e05b4a183d3655 --- /dev/null +++ b/python/corsika/io/outputs/observation_plane.py @@ -0,0 +1,87 @@ +""" + Read data written by ObservationPlane. + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import logging +import os.path as op +from typing import Any + +import pyarrow.parquet as pq + +from .output import Output + + +class ObservationPlane(Output): + """ + Read particle data from an ObservationPlane. + """ + + def __init__(self, path: str): + """ + Load the particle data into a parquet table. + + Parameters + ---------- + path: str + The path to the directory containing this output. + """ + super().__init__(path) + + # try and load our data + try: + self.__data = pq.read_table(op.join(path, "particles.parquet")) + except Exception as e: + logging.getLogger("corsika").warn( + f"An error occured loading an ObservationPlane: {e}" + ) + + def is_good(self) -> bool: + """ + Returns true if this output has been read successfully + and has the correct files/state/etc. + + Returns + ------- + bool: + True if this is a good output. + """ + return self.__data is not None + + def astype(self, dtype: str = "pandas", **kwargs: Any) -> Any: + """ + Load the particle data from this observation plane. + + All additional keyword arguments are passed to `parquet.read_table` + + Parameters + ---------- + dtype: str + The data format to return the data in (i.e. numpy, pandas, etc.) + + Returns + ------- + Any: + The return type of this method is determined by `dtype`. + """ + if dtype == "arrow": + return self.__data + elif dtype == "pandas": + return self.__data.to_pandas() + else: + raise ValueError( + ( + f"Unknown format '{dtype}' for ObservationPlane. " + "We currently only support ['arrow', 'pandas']." + ) + ) + + def __repr__(self) -> str: + """ + Return a string representation of this class. + """ + return f"ObservationPlane('{self.config['name']}')" diff --git a/python/corsika/io/outputs/output.py b/python/corsika/io/outputs/output.py new file mode 100644 index 0000000000000000000000000000000000000000..83a4c7928b786b9bad4301bec39b7e03aaf567f8 --- /dev/null +++ b/python/corsika/io/outputs/output.py @@ -0,0 +1,165 @@ +""" + This file defines the API for all output readers. + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import os.path as op +from abc import ABC, abstractmethod +from typing import Any, Dict + +import yaml + + +class Output(ABC): + """ + This class defines the abstract interface for all classes + that wish to provide reading support for CORSIKA8 outputs. + """ + + def __init__(self, path: str): + """ + __init__ must load the output files and check + that it is valid. + + Parameters + ---------- + path: str + The path to the directory containing this output. + """ + # load and store our path and config + self.path = path + self.__config = self.load_config(path) + self.__summary = self.load_summary(path) + + @abstractmethod + def is_good(self) -> bool: + """ + Returns true if this output has been read successfully + and has the correct files/state/etc. + + Returns + ------- + bool: + True if this is a good output. + """ + pass + + @abstractmethod + def astype(self, dtype: str, **kwargs: Any) -> Any: + """ + Return the data for this output in the data format given by 'dtype' + + Parameters + ---------- + dtype: str + The data format to return the data in (i.e. numpy, pandas, etc.) + *args: Any + Additional arguments can be accepted by the output. + **kwargs: Any + Additional keyword arguments can be accepted by the output. + + Returns + ------- + Any: + The return type of this method is determined by `dtype`. + """ + pass + + @property + def config(self) -> Dict[str, Any]: + """ + Return the config file for this output. + + Parameters + ---------- + + Returns + ------- + Dict[str, any] + The configuration file for this output. + """ + return self.__config + + @property + def summary(self) -> Dict[str, Any]: + """ + Return the summary file for this output. + + Parameters + ---------- + + Returns + ------- + Dict[str, any] + The summary file for this output. + """ + return self.__summary + + @property + def data(self) -> Any: + """ + Return the data in its default format. + + We try to use Pandas as the default format for most data. + + Parameters + ---------- + + Returns + ------- + Any: + The data in its default format. + """ + return self.astype() + + @staticmethod + def load_config(path: str) -> Dict[str, Any]: + """ + Load the top-level config from a given library path. + + Parameters + ---------- + path: str + The path to the directory containing the library. + + Returns + ------- + dict: + The config as a python dictionary. + + Raises + ------ + FileNotFoundError + If the config file cannot be found + + """ + with open(op.join(path, "config.yaml"), "r") as f: + return yaml.load(f, Loader=yaml.Loader) + + @staticmethod + def load_summary(path: str) -> Dict[str, Any]: + """ + Load the top-level summary from a given library path. + + Parameters + ---------- + path: str + The path to the directory containing the library. + + Returns + ------- + dict: + The summary as a python dictionary. + + Raises + ------ + FileNotFoundError + If the summary file cannot be found + + """ + with open(op.join(path, "summary.yaml"), "r") as f: + return yaml.load(f, Loader=yaml.Loader) diff --git a/python/corsika/io/outputs/track_writer.py b/python/corsika/io/outputs/track_writer.py new file mode 100644 index 0000000000000000000000000000000000000000..0859660d590182578e08ad4beaef5bb9a5c53498 --- /dev/null +++ b/python/corsika/io/outputs/track_writer.py @@ -0,0 +1,87 @@ +""" + Read data written by TrackWriter + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import logging +import os.path as op +from typing import Any + +import pyarrow.parquet as pq + +from .output import Output + + +class TrackWriter(Output): + """ + Read particle data from a TrackWriter + """ + + def __init__(self, path: str): + """ + Load the particle data into a parquet table. + + Parameters + ---------- + path: str + The path to the directory containing this output. + """ + super().__init__(path) + + # try and load our data + try: + self.__data = pq.read_table(op.join(path, "tracks.parquet")) + except Exception as e: + logging.getLogger("corsika").warn( + f"An error occured loading a TrackWriter: {e}" + ) + + def is_good(self) -> bool: + """ + Returns true if this output has been read successfully + and has the correct files/state/etc. + + Returns + ------- + bool: + True if this is a good output. + """ + return self.__data is not None + + def astype(self, dtype: str = "pandas", **kwargs: Any) -> Any: + """ + Load the particle data from this track writer. + + All additional keyword arguments are passed to `parquet.read_table` + + Parameters + ---------- + dtype: str + The data format to return the data in (i.e. numpy, pandas, etc.) + + Returns + ------- + Any: + The return type of this method is determined by `dtype`. + """ + if dtype == "arrow": + return self.__data + elif dtype == "pandas": + return self.__data.to_pandas() + else: + raise ValueError( + ( + f"Unknown format '{dtype}' for TrackWriter. " + "We currently only support ['arrow', 'pandas']." + ) + ) + + def __repr__(self) -> str: + """ + Return a string representation of this class. + """ + return f"TrackWriter('{self.config['name']}')" diff --git a/python/setup.cfg b/python/setup.cfg new file mode 100644 index 0000000000000000000000000000000000000000..f58d04eeb061f000b2d890aa26f542f70f3aa0c0 --- /dev/null +++ b/python/setup.cfg @@ -0,0 +1,69 @@ +[flake8] +# use a slightly longer line to be consistent with black +max-line-length = 88 + +# E231 is missing whitespace after a comma +# this contradicts the black formatting rules +# and therefore creates spurious errors +ignore = E231 + +# we set various directories we want to exclude +exclude = + # Don't bother checking in cache directories + __pycache__ + + +[isort] +# use parenthesis for multi-line imports +use_parentheses = true + + +[mypy] +# the primary Python version +python_version = 3.7 + +# allow returning Any +# this creates excessive errors when using libraries +# that don't have MyPy typing support +warn_return_any = False + +# don't allow untyped functions +disallow_untyped_defs = True + +# warn if any part of this config is mispelled +warn_unused_configs = True + +# warn for missing type information +warn_incomplete_stub = True + +# warn us if we don't return from a function explicitly +warn_no_return = True + +# use incremental typing to speed things up +incremental = True + +# show error contexts +show_error_context = True + +# and show the column numbers for errors +show_column_numbers = True + +# ignore missing types for setuptools +[mypy-setuptools.*] +ignore_missing_imports = True + +# ignore missing types for numpy +[mypy-numpy.*] +ignore_missing_imports = True + +# ignore missing types for matplotlib +[mypy-matplotlib.*] +ignore_missing_imports = True + +# ignore missing types for boost_histogram +[mypy-boost_histogram.*] +ignore_missing_imports = True + +# ignore missing types for pyarow +[mypy-pyarrow.*] +ignore_missing_imports = True diff --git a/python/setup.py b/python/setup.py new file mode 100644 index 0000000000000000000000000000000000000000..196c2ecbbccbc7e5caafee93be0af1c41746c94f --- /dev/null +++ b/python/setup.py @@ -0,0 +1,51 @@ +from os import path +from setuptools import setup + +# the stereo version +__version__ = "8.0.0-alpha" + +# get the absolute path of this project +here = path.abspath(path.dirname(__file__)) + +# Get the long description from the README file +with open(path.join(here, "README.md"), encoding="utf-8") as f: + long_description = f.read() + +# the standard setup info +setup( + name="corsika", + version=__version__, + description="A Python package for working with CORSIKA 8.", + long_description=long_description, + long_description_content_type="text/markdown", + url="https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika", + author="CORSIKA 8 Collaboration", + author_email="corsika-devel@lists.kit.edu", + classifiers=[ + "Development Status :: 3 - Alpha", + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Physics", + "Programming Language :: Python :: 3.6", + "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.8", + ], + keywords=["cosmic ray", "physics", "air shower", "simulation"], + packages=["corsika"], + python_requires=">=3.6*, <4", + install_requires=["numpy", "pyyaml", "pyarrow", "boost_histogram"], + extras_require={ + "test": [ + "pytest", + "black", + "mypy", + "isort", + "coverage", + "pytest-cov", + "flake8", + ], + "pandas": ["pandas"], + }, + scripts=[], + project_urls={"code": "https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika"}, + include_package_data=False, +) diff --git a/python/tests/__init__.py b/python/tests/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..c7d9b672b4c1bde7cc7725c619c0390a79e682c7 --- /dev/null +++ b/python/tests/__init__.py @@ -0,0 +1,56 @@ +""" + Tests for `corsika` + + (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" + +import os +import os.path as op + +__all__ = ["build_directory"] + + +def find_build_directory() -> str: + """ + Return the absolute path to the CORSIKA8 build directory. + + Returns + ------- + str + The absolute path to the build directory. + + Raises + ------ + RuntimeError + If the build directory cannot be found or it is empty. + """ + + # check if we are running on Gitlab + gitlab_build_dir = os.getenv("CI_BUILDS_DIR") + + # if we are running on Gitlab + if gitlab_build_dir is not None: + build_dir = op.abspath(gitlab_build_dir) + else: # otherwise, we are running locally + build_dir = op.abspath( + op.join( + op.dirname(__file__), + op.pardir, + op.pardir, + "build", + ) + ) + + # check that the build directory contains 'CMakeCache.txt' + if not op.exists(op.join(build_dir, "CMakeCache.txt")): + raise RuntimeError("Python tests cannot find C8 build directory.") + + return build_dir + + +# find the build_directory once +build_directory = find_build_directory() diff --git a/python/tests/io/__init__.py b/python/tests/io/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..0d08be944644953fd08405a4502456e76c1203b9 --- /dev/null +++ b/python/tests/io/__init__.py @@ -0,0 +1,9 @@ +""" + Tests for `corsika.io` + + (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" diff --git a/python/tests/io/test_hist.py b/python/tests/io/test_hist.py new file mode 100644 index 0000000000000000000000000000000000000000..7cf5b38fb7798c81680cf14d108e3b31dd3dcaf1 --- /dev/null +++ b/python/tests/io/test_hist.py @@ -0,0 +1,70 @@ +""" + Tests for `corsika.io.hist` + + (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import os +import os.path as op +import subprocess + +import corsika + +from .. import build_directory + +# the directory containing 'testSaveBoostHistogram' +bindir = op.join(build_directory, "Framework", "Utilities") + + +def generate_hist() -> str: + """ + Generate a test with `testSaveBoostHistogram`. + + Returns + ------- + str + The path to the generated histogram. + bool + If True, this file was regenerated. + """ + + # we construct the name of the bin + bin = op.join(bindir, "testSaveBoostHistogram") + + # check if a histogram already exists + if op.exists(op.join(bin, "hist.npz")): + return op.join(bin, "hist.npz"), False + else: + # run the program - this generates "hist.npz" in the CWD + subprocess.call(bin) + + return op.join(os.getcwd(), "hist.npz"), True + + +def test_corsika_io() -> None: + """ + Test I can corsika.io without a further import. + """ + corsika.io.read_hist + + +def test_corsika_read_hist() -> None: + """ + Check that I can read in the test histograms with `read_hist`. + """ + + # generate a test histogram + filename, delete = generate_hist() + + # and try and read in a histogram + h = corsika.io.read_hist(filename) + + # and delete the generated histogram + if delete: + os.remove(filename) + + # and check that it isn't empty + assert not h.empty() diff --git a/src/framework/core/CMakeLists.txt b/src/framework/core/CMakeLists.txt index 558b227f68801ec98c794c832ae5d8e654c1a703..f0edcdebd80fd61962c258b8e8ab43963030e5f6 100644 --- a/src/framework/core/CMakeLists.txt +++ b/src/framework/core/CMakeLists.txt @@ -5,7 +5,7 @@ file (MAKE_DIRECTORY ${output_dir}) add_custom_command ( OUTPUT ${output_dir}/GeneratedParticleProperties.inc - OUTPUT ${output_dir}/GeneratedParticleClasses.inc + ${output_dir}/GeneratedParticleClasses.inc ${output_dir}/particle_db.pkl COMMAND ${input_dir}/pdxml_reader.py ${input_dir}/ParticleData.xml ${input_dir}/NuclearData.xml @@ -22,7 +22,9 @@ add_custom_command ( add_custom_target (GenParticlesHeaders DEPENDS ${output_dir}/GeneratedParticleProperties.inc - ${output_dir}/GeneratedParticleClasses.inc) + ${output_dir}/GeneratedParticleClasses.inc + ${output_dir}/particle_db.pkl + ) add_dependencies (CORSIKA8 GenParticlesHeaders) install ( diff --git a/src/modules/qgsjetII/CMakeLists.txt b/src/modules/qgsjetII/CMakeLists.txt index 36ebe9b693b927de9acfca03be3dc016440c4050..9b9777356e511db01238c2a68389c5c476e738a7 100644 --- a/src/modules/qgsjetII/CMakeLists.txt +++ b/src/modules/qgsjetII/CMakeLists.txt @@ -10,7 +10,7 @@ add_custom_command ( ${input_dir}/qgsjet-II-04-codes.dat DEPENDS ${input_dir}/code_generator.py ${input_dir}/qgsjet-II-04-codes.dat - GenParticlesHeaders # for particle_db.pkl + ${PROJECT_BINARY_DIR}/corsika/framework/core/particle_db.pkl WORKING_DIRECTORY ${output_dir} COMMENT "Generate conversion tables for particle codes QGSJetII <-> CORSIKA" diff --git a/src/modules/sibyll/CMakeLists.txt b/src/modules/sibyll/CMakeLists.txt index b3d97a12349a5a5f2cbd9917404a8685894f0096..969d2f1d2e0e1bada1a4c247901ef381c121d368 100644 --- a/src/modules/sibyll/CMakeLists.txt +++ b/src/modules/sibyll/CMakeLists.txt @@ -10,7 +10,7 @@ add_custom_command ( ${input_dir}/sibyll_codes.dat DEPENDS ${input_dir}/code_generator.py ${input_dir}/sibyll_codes.dat - GenParticlesHeaders # for particle_db.pkl + ${PROJECT_BINARY_DIR}/corsika/framework/core/particle_db.pkl WORKING_DIRECTORY ${output_dir}/ COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA" diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3a12d7cc7401177052d3374be42c583073b3ae20..b8f1e709815fa0d41f9628482981cc4d879caec9 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -3,3 +3,4 @@ add_subdirectory (framework) add_subdirectory (media) add_subdirectory (stack) add_subdirectory (modules) +add_subdirectory (output) diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt index 35f61737b13dfedef24d6746cf92e56a26658ebe..829709779bd7ec19b8b9fb0037422a1f6ca105da 100644 --- a/tests/framework/CMakeLists.txt +++ b/tests/framework/CMakeLists.txt @@ -20,6 +20,7 @@ set (test_framework_sources testNullModel.cpp testHelix.cpp testInteractionCounter.cpp + testSolver.cpp ) CORSIKA_ADD_TEST (testFramework SOURCES ${test_framework_sources}) diff --git a/tests/framework/testCOMBoost.cpp b/tests/framework/testCOMBoost.cpp index caf6f08965bf4468294df23dee13e4f265dd4c47..d092a4b61ea46dfbcb734a4c8f31013550451204 100644 --- a/tests/framework/testCOMBoost.cpp +++ b/tests/framework/testCOMBoost.cpp @@ -40,7 +40,6 @@ auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) { TEST_CASE("rotation") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); // define projectile kinematics in lab frame HEPMassType const projectileMass = 1_GeV; @@ -172,7 +171,6 @@ TEST_CASE("rotation") { TEST_CASE("boosts") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); // define target kinematics in lab frame HEPMassType const targetMass = 1_GeV; @@ -329,7 +327,6 @@ TEST_CASE("boosts") { TEST_CASE("rest frame") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); HEPMassType const projectileMass = 1_GeV; HEPMomentumType const P0 = 1_TeV; diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp index b860e9955d3c3122c4295cdaa6f32875f9aa6021..44d277bae1bf03a7fdfe9fa12eef3ae2c2197f95 100644 --- a/tests/framework/testCascade.cpp +++ b/tests/framework/testCascade.cpp @@ -24,6 +24,8 @@ #include <corsika/media/HomogeneousMedium.hpp> #include <corsika/media/NuclearComposition.hpp> +#include <corsika/output/DummyOutputManager.hpp> + #include <SetupTestTrajectory.hpp> #include <catch2/catch.hpp> @@ -137,7 +139,6 @@ public: TEST_CASE("Cascade", "[Cascade]") { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); logging::set_level(logging::level::trace); HEPEnergyType E0 = 100_GeV; @@ -164,8 +165,10 @@ TEST_CASE("Cascade", "[Cascade]") { Point(rootCS, {0_m, 0_m, 10_km}), 0_ns)); DummyTracking tracking; - Cascade<DummyTracking, decltype(sequence), TestCascadeStack, TestCascadeStackView> EAS( - env, tracking, sequence, stack); + DummyOutputManager output; + Cascade<DummyTracking, decltype(sequence), DummyOutputManager, TestCascadeStack, + TestCascadeStackView> + EAS(env, tracking, sequence, output, stack); SECTION("full cascade") { EAS.run(); diff --git a/tests/framework/testClassTimer.cpp b/tests/framework/testClassTimer.cpp index d51a56b1a254146847f742bfc0283ee9f6e44be2..b88d2fc97a0f2c99e5100268e168b6134598ca63 100644 --- a/tests/framework/testClassTimer.cpp +++ b/tests/framework/testClassTimer.cpp @@ -106,7 +106,6 @@ public: TEST_CASE("ClassTimer", "[Timer]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("Measure runtime of a function without arguments") { diff --git a/tests/framework/testCombinedStack.cpp b/tests/framework/testCombinedStack.cpp index 7d01f3fcd3139507c41617a2fd276ddfedb46138..de6defc6ad7a0f94fcbcc289a012de38761112dd 100644 --- a/tests/framework/testCombinedStack.cpp +++ b/tests/framework/testCombinedStack.cpp @@ -89,7 +89,6 @@ using StackTest = CombinedStack<TestStackData, TestStackData2, CombinedTestInter TEST_CASE("Combined Stack", "[stack]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); // helper function for sum over stack data auto sum = [](const StackTest& stack) { @@ -218,6 +217,16 @@ TEST_CASE("Combined Stack", "[stack]") { CHECK(s.getEntries() == 0); CHECK(s.isEmpty()); } + + SECTION("exceptions") { + StackTest s; + auto p1 = s.addParticle(std::tuple{9.9}); + auto p2 = s.addParticle(std::tuple{9.9}); + ++p2; + CHECK_THROWS(s.copy(p1, p2)); + CHECK_THROWS(s.swap(p1, p2)); + CHECK(s.getSize() == 2); + } } //////////////////////////////////////////////////////////// @@ -288,7 +297,6 @@ using StackTest2 = CombinedStack<typename StackTest::stack_implementation_type, TEST_CASE("Combined Stack - multi", "[stack]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("create secondaries") { @@ -386,7 +394,6 @@ using Particle2 = typename StackTest2::particle_type; TEST_CASE("Combined Stack - secondary view") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("create secondaries via secondaryview") { diff --git a/tests/framework/testCorsikaFenv.cpp b/tests/framework/testCorsikaFenv.cpp index 2b6cdb58034d06c82cf336c2826cdb2e4f4248f6..030d33dfac78b92464b677b489f7f62b8fb6fa4c 100644 --- a/tests/framework/testCorsikaFenv.cpp +++ b/tests/framework/testCorsikaFenv.cpp @@ -22,7 +22,6 @@ static void handle_fpe(int /*signo*/) { gRESULT = 0; } TEST_CASE("CorsikaFenv", "[fenv]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("Enable all exceptions") { feenableexcept(FE_ALL_EXCEPT); } diff --git a/tests/framework/testFourVector.cpp b/tests/framework/testFourVector.cpp index 09b458295ef0b118684ed2707389a0a5f0fcfd9c..489d6261b17d6c68f85a01dae591992cbbba6a90 100644 --- a/tests/framework/testFourVector.cpp +++ b/tests/framework/testFourVector.cpp @@ -20,7 +20,6 @@ using namespace corsika; TEST_CASE("four vectors") { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); logging::set_level(logging::level::info); // this is just needed as a baseline diff --git a/tests/framework/testFunctionTimer.cpp b/tests/framework/testFunctionTimer.cpp index 46d97f710f5793a5f27bbc4ed255c49cbdad627b..1328ea7badbc86f13809197a06ca773593381137 100644 --- a/tests/framework/testFunctionTimer.cpp +++ b/tests/framework/testFunctionTimer.cpp @@ -32,7 +32,6 @@ public: TEST_CASE("FunctionTimer", "[Timer]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("Measure runtime of a free function") { diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp index e8847ee5e824242344147f46afeb37dfe0eb3708..809ae5e1391071b948a6307600cabb4bfc201136 100644 --- a/tests/framework/testGeometry.cpp +++ b/tests/framework/testGeometry.cpp @@ -14,6 +14,7 @@ #include <corsika/framework/geometry/Line.hpp> #include <corsika/framework/geometry/Helix.hpp> #include <corsika/framework/geometry/Point.hpp> +#include <corsika/framework/geometry/Path.hpp> #include <corsika/framework/geometry/RootCoordinateSystem.hpp> #include <corsika/framework/geometry/Sphere.hpp> #include <corsika/framework/geometry/StraightTrajectory.hpp> @@ -312,3 +313,78 @@ TEST_CASE("Geometry Trajectories") { .magnitude() == Approx(0).margin(absMargin)); } } + +TEST_CASE("Distance between points") { + // define a known CS + CoordinateSystemPtr root = get_root_CoordinateSystem(); + + // define known points + Point p1(root, {0_m, 0_m, 0_m}); + Point p2(root, {0_m, 0_m, 5_m}); + Point p3(root, {1_m, 0_m, 0_m}); + Point p4(root, {5_m, 0_m, 0_m}); + Point p5(root, {0_m, 4_m, 0_m}); + Point p6(root, {0_m, 5_m, 0_m}); + + // check distance() function + CHECK(distance(p1, p2) / 1_m == Approx(5)); + CHECK(distance(p3, p4) / 1_m == Approx(4)); + CHECK(distance(p5, p6) / 1_m == Approx(1)); +} + +TEST_CASE("Path") { + // define a known CS + CoordinateSystemPtr root = get_root_CoordinateSystem(); + + // define known points + Point p1(root, {0_m, 0_m, 0_m}); + Point p2(root, {0_m, 0_m, 1_m}); + Point p3(root, {0_m, 0_m, 2_m}); + Point p4(root, {0_m, 0_m, 3_m}); + Point p5(root, {0_m, 0_m, 4_m}); + // define paths + Path P1(p1); + Path P2({p1, p2}); + Path P3({p1, p2, p3}); + // define deque that include point(s) + std::deque<Point> l1 = {p1}; + std::deque<Point> l2 = {p1, p2}; + std::deque<Point> l3 = {p1, p2, p3}; + + // test the various path constructors + SECTION("Test Constructors") { + // check constructor for one point + CHECK(std::equal(P1.begin(), P1.end(), l1.begin(), + [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; })); + // check constructor for collection of points + CHECK(std::equal(P3.begin(), P3.end(), l3.begin(), + [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; })); + } + + // test the length and access methods + SECTION("Test getLength() and modifications to Path") { + P1.addToEnd(p2); + P2.removeFromEnd(); + // Check modifications to path + CHECK(std::equal(P1.begin(), P1.end(), l2.begin(), + [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; })); + CHECK(std::equal(P2.begin(), P2.end(), l1.begin(), + [](Point a, Point b) { return (a - b).getNorm() / 1_m < 1e-5; })); + // Check GetStart(), GetEnd(), GetPoint() + CHECK((P3.getEnd() - P3.getStart()).getNorm() / 1_m == Approx(2)); + CHECK((P1.getPoint(1) - p2).getNorm() / 1_m == Approx(0)); + // Check GetLength() + CHECK(P1.getLength() / 1_m == Approx(1)); + CHECK(P2.getLength() / 1_m == Approx(0)); + CHECK(P3.getLength() / 1_m == Approx(2)); + P2.removeFromEnd(); + CHECK(P2.getLength() / 1_m == Approx(0)); // Check the length of an empty path + P3.addToEnd(p4); + P3.addToEnd(p5); + CHECK(P3.getLength() / 1_m == Approx(4)); + P3.removeFromEnd(); + CHECK(P3.getLength() / 1_m == Approx(3)); // Check RemoveFromEnd() else case + // Check GetNSegments() + CHECK(P3.getNSegments() - 3 == Approx(0)); + } +} \ No newline at end of file diff --git a/tests/framework/testHelix.cpp b/tests/framework/testHelix.cpp index cbe997d203a02705001f9304bc767cfedce848db..042e4212f122bb538c86e9c78317fb86259b3ba8 100644 --- a/tests/framework/testHelix.cpp +++ b/tests/framework/testHelix.cpp @@ -22,7 +22,6 @@ double constexpr absMargin = 1.0e-8; TEST_CASE("Helix class") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); const CoordinateSystemPtr rootCS = get_root_CoordinateSystem(); Point r0(rootCS, {0_m, 0_m, 0_m}); diff --git a/tests/framework/testInteractionCounter.cpp b/tests/framework/testInteractionCounter.cpp index 556d4733b60b82bf295ff47463a24d3775c48553..4b4b7ff25dca0ba6a5cf44018c9efde66a150e82 100644 --- a/tests/framework/testInteractionCounter.cpp +++ b/tests/framework/testInteractionCounter.cpp @@ -43,7 +43,6 @@ struct DummyProcess { TEST_CASE("InteractionCounter", "[process]") { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); logging::set_level(logging::level::info); DummyProcess d; diff --git a/tests/framework/testLogging.cpp b/tests/framework/testLogging.cpp index b737d7dfa778b0cdf78b235479177e51534311b2..bf92cf8df00e8f12d64c0d5779f59b7976c4e09d 100644 --- a/tests/framework/testLogging.cpp +++ b/tests/framework/testLogging.cpp @@ -14,8 +14,6 @@ using namespace corsika; TEST_CASE("Logging", "[Logging]") { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); - SECTION("top level functions using default corsika logger") { logging::info("(1) This is an info message!"); logging::warn("(1) This is a warning message!"); diff --git a/tests/framework/testNullModel.cpp b/tests/framework/testNullModel.cpp index 5b1f351e13638a1ac2358bd5f0b9d2f295563350..872fc34dbbeb685db1806c74c93f4709443db545 100644 --- a/tests/framework/testNullModel.cpp +++ b/tests/framework/testNullModel.cpp @@ -21,7 +21,6 @@ using namespace corsika; TEST_CASE("NullModel", "[processes]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("interface") { [[maybe_unused]] NullModel nm; diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp index 9e87adbf40e46db7a4c5bc0bbee4fd870e87fcc2..c35f8bca97982e2d2388ab122cfad22269dd4960 100644 --- a/tests/framework/testParticles.cpp +++ b/tests/framework/testParticles.cpp @@ -17,7 +17,6 @@ using namespace corsika; TEST_CASE("ParticleProperties", "[Particles]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("Types") { CHECK(Electron::code == Code::Electron); diff --git a/tests/framework/testRandom.cpp b/tests/framework/testRandom.cpp index 5774575e6844754289cdcd2130313d3b03d95359..2d0c0720ae163f2777c80837d8179b9355aba044 100644 --- a/tests/framework/testRandom.cpp +++ b/tests/framework/testRandom.cpp @@ -55,7 +55,6 @@ SCENARIO("random-number streams can be registered and retrieved") { TEST_CASE("UniformRealDistribution") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); std::mt19937 rng; @@ -101,7 +100,6 @@ TEST_CASE("UniformRealDistribution") { TEST_CASE("ExponentialDistribution") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); std::mt19937 rng; diff --git a/tests/framework/testSaveBoostHistogram.cpp b/tests/framework/testSaveBoostHistogram.cpp index d5c16438a2ab6de1595393205d4eff4320a43d08..89182311a650e76968226e3dc2cf5a97c8b4e29a 100644 --- a/tests/framework/testSaveBoostHistogram.cpp +++ b/tests/framework/testSaveBoostHistogram.cpp @@ -16,7 +16,6 @@ namespace bh = boost::histogram; TEST_CASE("SaveHistogram") { - corsika::corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); corsika::logging::set_level(corsika::logging::level::info); std::mt19937 rng; diff --git a/tests/framework/testSecondaryView.cpp b/tests/framework/testSecondaryView.cpp index 87cfd91c0c74cb3b9c38cc8bda31af7a96e52c8c..fa113f88b6e111e9837d9ef96d28a432024974dd 100644 --- a/tests/framework/testSecondaryView.cpp +++ b/tests/framework/testSecondaryView.cpp @@ -44,7 +44,6 @@ using Particle = typename StackTest::particle_type; TEST_CASE("SecondaryStack", "[stack]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); // helper function for sum over stack data auto sum = [](const StackTest& stack) { diff --git a/tests/framework/testSolver.cpp b/tests/framework/testSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..52f1fa13fc29f59d7b90e4471b76bb51997b0613 --- /dev/null +++ b/tests/framework/testSolver.cpp @@ -0,0 +1,278 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <cmath> +#include <vector> +#include <complex> +#include <algorithm> + +#include <corsika/framework/core/Logging.hpp> +#include <corsika/framework/utility/QuadraticSolver.hpp> +#include <corsika/framework/utility/LinearSolver.hpp> +#include <corsika/framework/utility/CubicSolver.hpp> +#include <corsika/framework/utility/QuarticSolver.hpp> + +using namespace corsika; +using namespace std; + +double pol4(long double x, long double a, long double b, long double c, long double d, + long double e) { + return std::pow(x, 4) * a + std::pow(x, 3) * b + std::pow(x, 2) * c + x * d + e; +} + +void remove_duplicates(vector<double>& v, double const eps) { + std::sort(v.begin(), v.end()); + v.erase(std::unique(v.begin(), v.end(), + [eps](auto v1, auto v2) { + return (std::abs(v2) > eps ? std::abs(2 * (v1 - v2) / (v1 + v2)) < + eps // relative + : std::abs(v1 - v2) < eps); // absolute + }), + v.end()); +} + +TEST_CASE("Solver") { + + logging::set_level(logging::level::trace); + + double epsilon_check = 1e-3; // for catch2 asserts + + SECTION("linear") { + + std::vector<std::pair<double, double>> vals = {{13.33, 8338.3}, + {-333.99, -8.15}, + {-58633.3, 2343.54}, + {87632.87, -982.37}, + {7e8, -1e-7}}; + + for (auto v : vals) { + + { + double a = v.first; + double b = v.second; + + vector<double> s1 = solve_linear_real(a, b); + + CORSIKA_LOG_INFO("linear: a={}, b={}, N={}, s1[0]={}", a, b, s1.size(), s1[0]); + + double expected = -b / a; + + CHECK(s1.size() == 1); + CHECK((s1[0] == Approx(expected).epsilon(epsilon_check))); + + vector<complex<double>> s2 = solve_linear(a, b); + CHECK(s2.size() == 1); + CHECK((s2[0].real() == Approx(expected).epsilon(epsilon_check))); + } + } + + CHECK(solve_linear_real(0, 55.).size() == 0); + } + + SECTION("quadratic") { + + // tests of type: (x-z1)*(x-z2) = 0 --> x1=z1, x2=z2 and a=1, b=-z2-z1, + // c=z1*z2 + + std::vector<std::vector<double>> zeros = { + {13.33, 8338.3}, {-333.99, -8.15}, {-58633.3, 2343.54}, {87632.87, -982.37}, + {1.3e-6, 5e7}, {-4.2e-7, 65e6}, {0.1, -13e-5}, {7e8, -1e-7}}; + + for (unsigned int idegree = 0; idegree < 2; ++idegree) { + CORSIKA_LOG_INFO("------------- quadratic idegree={}", idegree); + for (auto z : zeros) { + + { + long double const z1 = z[0]; + long double const z2 = idegree <= 0 ? z1 : z[1]; + + long double a = 1; + long double b = -z1 - z2; + long double c = z1 * z2; + + std::vector<double> s1 = solve_quadratic_real(a, b, c); + remove_duplicates(s1, epsilon_check); + + CORSIKA_LOG_INFO("quadratic: z1={}, z2={}, N={}, s1=[{}]", z1, z2, s1.size(), + fmt::join(s1, ", ")); + + CHECK(s1.size() == idegree + 1); + for (auto value : s1) { + if (std::abs(value) < epsilon_check) { + CHECK(((value == Approx(z1).margin(epsilon_check)) || + (value == Approx(z2).margin(epsilon_check)))); + } else { + CHECK(((value == Approx(z1).epsilon(epsilon_check)) || + (value == Approx(z2).epsilon(epsilon_check)))); + } + } + /* + std::vector<complex<double>> s2 = solve_quadratic(a, b, c, epsilon); + CHECK(s2.size() == idegree + 1); + for (auto value : s2) { + if (std::abs(value) < epsilon) { + CHECK(((value == Approx(z1).margin(epsilon_check)) || + (value == Approx(z2).margin(epsilon_check)))); + } else { + CHECK(((value == Approx(z1).epsilon(epsilon_check)) || + (value == Approx(z2).epsilon(epsilon_check)))); + } + } + }*/ + } + } + } + } + + SECTION("cubic") { + + // tests of type: + // (x-z1) * (x-z2) * (x-z3) = 0 --> x1=z1, x2=z2, x3=z3 and + + // (x^2 - x*(z1+z2) + z1*z2) * (x-z3) = + // x^3 - x^2*(z1+z2) + x*z1*z2 - x^2 z3 + x*(z1+z2)*z3 - z1*z2*z3 = + // x^3 + x^2 (-z3-z1-z2) + x (z1*z2 + (z1+z2)*z3) - z1*z2*z3 + + epsilon_check = 1e-2; // for catch2 asserts + + std::vector<std::vector<double>> zeros = { + {13.33, 838.3, 44.}, {-333.99, -8.15, -33e8}, {-58633.3, 2343.54, -1e-5}, + {87632.87, -982.37, 0.}, {1.3e-4, 5e7, 6.6e9}, {-4.2e-7, 65e6, 433.3334}, + {23e-1, -13e-2, 10.}, {7e8, -1e-7, 1e8}}; + + for (unsigned int idegree = 0; idegree < 3; ++idegree) { + CORSIKA_LOG_INFO("------------- cubic idegree={}", idegree); + for (auto z : zeros) { + + { + long double const z1 = z[0]; + long double const z2 = idegree <= 0 ? z1 : z[1]; + long double const z3 = idegree <= 1 ? z1 : z[2]; + + long double const a = 1; + long double const b = -z1 - z2 - z3; + long double const c = z1 * z2 + (z1 + z2) * z3; + long double const d = -z1 * z2 * z3; + // + CORSIKA_LOG_INFO( + "cubic: z1={}, z2={}, z3={}, " + "a={}, b={}, c={}, d={}", + z1, z2, z3, a, b, c, d); + + vector<double> s1 = solve_cubic_real(a, b, c, d); + remove_duplicates(s1, epsilon_check); + + CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", ")); + + CHECK(s1.size() == idegree + 1); + for (double value : s1) { + CORSIKA_LOG_INFO("value={}, z1={} z2={} z3={} eps_check={}", value, z1, z2, + z3, epsilon_check); + if (std::abs(value) < epsilon_check) { + CHECK(((value == Approx(z1).margin(epsilon_check)) || + (value == Approx(z2).margin(epsilon_check)) || + (value == Approx(z3).margin(epsilon_check)))); + } else { + CHECK(((value == Approx(z1).epsilon(epsilon_check)) || + (value == Approx(z2).epsilon(epsilon_check)) || + (value == Approx(z3).epsilon(epsilon_check)))); + } + } + } + } + } + } // cubic + + SECTION("quartic") { + + epsilon_check = 1e-4; // for catch2 asserts + + // **clang-format-off** + // tests of type: + // (x-z1) * (x-z2) * (x-z3) * (x-z4) = 0 --> x1=z1, x2=z2, x3=z3, x4=z4 and + + // (x^2 - x (z1+z2) + z1*z2) * (x-z3) * + // (x-z4) = (x^3 - x^2 (z1+z2) + x z1*z2 - x^2 z3 + x*(z1+z2)*z3 - + // z1*z2*z3) * (x-z4) = (x^3 + x^2 (-z1-z2-z3) + x (z1*z2 + (z1+z2)*z3) - + // z1*z2*z3) * (x-z4) + // = + // + // x^4 + x^3 (-z1-z2-z3) + x^2 (z1*z2 + (z1+z2)*z3) - x z1*z2*z3 + // - x^3 z4 - x^2 (-z1-z2-z3)*z4 - x (z1*z2 + (z1+z2)*z3)*z4 + // + z1*z2*z3*z4 + // + // = x^4 + // + x^3 (-z1-z2-z3-z4) + // + x^2 (z1*z2 + (z1+z2)*z3 + (z1+z2+z3)*z4) + // - x (z1*z2*z3 + (z1*z2 + (z1+z2)*z3)*z4)) + // + z1*z2*z3*z4 + // **clang-format-on** + + std::vector<std::vector<double>> zeros = { + {13.33, 838.3, 44., 2.3}, {-3333.99, -8.15, -33e4, 8.8e3}, + {-58633.3, 2343.54, -1e-1, 0.}, {87632.87, -982.37, 10., 1e-2}, + {1.3e2, 5e5, 6.6e5, 1e3}, {-4.9, 65e2, 433.3334, 6e5}, + {23e-1, -13e-2, 10., 3.4e6}, {7e6, -1e-1, 1e5, 2.55e4}}; + + for (unsigned int idegree = 2; idegree < 4; + ++idegree) { // idegree==1 is not very stable !!!!!!! + CORSIKA_LOG_INFO("------------- quartic idegree={}", idegree); + for (auto z : zeros) { + + { + long double const z1 = z[0]; + long double const z2 = idegree <= 0 ? z1 : z[1]; + long double const z3 = idegree <= 1 ? z1 : z[2]; + long double const z4 = idegree <= 2 ? z1 : z[3]; + + long double const a = 1; + long double const b = -z1 - z2 - z3 - z4; + long double const c = z1 * z2 + (z1 + z2) * z3 + (z1 + z2 + z3) * z4; + long double const d = -z1 * z2 * z3 - (z1 * z2 + (z1 + z2) * z3) * z4; + long double const e = z1 * z2 * z3 * z4; + + // + CORSIKA_LOG_INFO( + "quartic: z1={}, z2={}, z3={}, z4={}, " + "a={}, b={}, c={}, d={}, e={}", + z1, z2, z3, z4, a, b, c, d, e); + + // make sure the tested cases are all ZERO (printout) + CORSIKA_LOG_INFO("quartic trace: {} {} {} {}", pol4(z1, a, b, c, d, e), + pol4(z2, a, b, c, d, e), pol4(z3, a, b, c, d, e), + pol4(z4, a, b, c, d, e)); + + vector<double> s1 = andre::solve_quartic_real(a, b, c, d, e); + // vector<double> s1 = solve_quartic_real(a, b, c, d, e, epsilon); + remove_duplicates(s1, epsilon_check); + + CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", ")); + + CHECK(s1.size() == idegree + 1); + for (double value : s1) { + CORSIKA_LOG_INFO("value={}, z1={} z2={} z3={} z4={} eps_check={}", value, z1, + z2, z3, z4, epsilon_check); + if (std::abs(value) < epsilon_check) { + CHECK(((value == Approx(z1).margin(epsilon_check)) || + (value == Approx(z2).margin(epsilon_check)) || + (value == Approx(z3).margin(epsilon_check)) || + (value == Approx(z4).margin(epsilon_check)))); + } else { + CHECK(((value == Approx(z1).epsilon(epsilon_check)) || + (value == Approx(z2).epsilon(epsilon_check)) || + (value == Approx(z3).epsilon(epsilon_check)) || + (value == Approx(z4).epsilon(epsilon_check)))); + } + } + } + } + } + } // quartic +} diff --git a/tests/framework/testStackInterface.cpp b/tests/framework/testStackInterface.cpp index b37dc495d7b03850ce135dd87942f8227531963d..4171647af5a547baeb55e6a2e3c3701eaeb77b29 100644 --- a/tests/framework/testStackInterface.cpp +++ b/tests/framework/testStackInterface.cpp @@ -27,7 +27,6 @@ typedef Stack<TestStackData, TestParticleInterface> StackTest; TEST_CASE("Stack", "[Stack]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); // helper function for sum over stack data auto sum = [](const StackTest& stack) { diff --git a/tests/framework/testUnits.cpp b/tests/framework/testUnits.cpp index 8d15544aca0d5c2dc1923af5853e6de7cafeb4f4..5c90aa9cd6ad65500a13a178d7dd6b0e52355006 100644 --- a/tests/framework/testUnits.cpp +++ b/tests/framework/testUnits.cpp @@ -19,7 +19,6 @@ using namespace corsika; TEST_CASE("PhysicalUnits", "[Units]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("Consistency") { CHECK(1_m / 1_m == Approx(1)); diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp index fb1e71fec7b914fcb449f8a2bd2d06dce952222e..1aee5973e8dca12f3f45b743e47658ef324e5600 100644 --- a/tests/media/testRefractiveIndex.cpp +++ b/tests/media/testRefractiveIndex.cpp @@ -16,6 +16,7 @@ #include <corsika/media/InhomogeneousMedium.hpp> #include <corsika/media/NuclearComposition.hpp> #include <corsika/media/UniformRefractiveIndex.hpp> +#include <corsika/media/ExponentialRefractiveIndex.hpp> #include <SetupTestTrajectory.hpp> @@ -89,3 +90,81 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") { CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); } + +TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") { + + logging::set_level(logging::level::info); + corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); + + // get a CS and a point + CoordinateSystemPtr const& gCS = get_root_CoordinateSystem(); + + Point const gOrigin(gCS, {0_m, 0_m, 0_m}); + + // setup our interface types + using IModelInterface = IRefractiveIndexModel<IMediumModel>; + using AtmModel = ExponentialRefractiveIndex<HomogeneousMedium<IModelInterface>>; + + // the constant density + const auto density{19.2_g / cube(1_cm)}; + + // the composition we use for the homogenous medium + NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, + std::vector<float>{1.f}); + + // a new refractive index + const double n0{2}; + const InverseLengthType lambda{56 / 1_m}; + + // create the atmospheric model and check refractive index + AtmModel medium(n0, lambda, density, protonComposition); + CHECK(n0 == medium.getRefractiveIndex(Point(gCS, 10_m, 14_m, 0_m))); + + // another refractive index + const double n0_{1}; + const InverseLengthType lambda_{1 / 1_km}; + + // create the atmospheric model and check refractive index + AtmModel medium_(n0_, lambda_, density, protonComposition); + CHECK(medium_.getRefractiveIndex(Point(gCS, 12_m, 56_m, 1_km)) == Approx(0.3678794)); + + // another refractive index + const double n0__{4}; + const InverseLengthType lambda__{2 / 1_m}; + + // create the atmospheric model and check refractive index + AtmModel medium__(n0__, lambda__, density, protonComposition); + CHECK(medium__.getRefractiveIndex(Point(gCS, 0_m, 0_m, 35_km)) == Approx(0)); + + // define our axis vector + Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1)); + + // check the density and nuclear composition + REQUIRE(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); + medium.getNuclearComposition(); + REQUIRE(density == medium_.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); + medium_.getNuclearComposition(); + REQUIRE(density == medium__.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); + medium__.getNuclearComposition(); + + // create a line of length 1 m + Line const line(gOrigin, Vector<SpeedType::dimension_type>( + gCS, {1_m / second, 0_m / second, 0_m / second})); + + // the end time of our line + auto const tEnd = 1_s; + + // and the associated trajectory + setup::Trajectory const track = + setup::testing::make_track<setup::Trajectory>(line, tEnd); + // // and the associated trajectory + // Trajectory<Line> const trajectory(line, tEnd); + + // and check the integrated grammage + REQUIRE((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); + REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); + REQUIRE((medium_.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); + REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); + REQUIRE((medium__.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); + REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); +} diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp index 19ab9976871312ca01f9ea45fc40a08cc4cc1bf7..63ca1fd573b3e7e2d0bbfff907d75ea609c7f6f2 100644 --- a/tests/modules/testCONEX.cpp +++ b/tests/modules/testCONEX.cpp @@ -39,7 +39,6 @@ using MExtraEnvirnoment = MediumPropertyModel<UniformMagneticField<T>>; TEST_CASE("CONEXSourceCut") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); RNGManager::getInstance().registerRandomStream("cascade"); RNGManager::getInstance().registerRandomStream("sibyll"); @@ -123,7 +122,6 @@ TEST_CASE("CONEXSourceCut") { TEST_CASE("ConexOutput", "[output validation]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); auto file = GENERATE(as<std::string>{}, "conex_fit", "conex_output"); diff --git a/tests/modules/testExecTime.cpp b/tests/modules/testExecTime.cpp index 68a85d81b2af7343ed323777daef39f55286a69e..4ebc03342914dc28e9137e7fd96dde010e8ac544 100644 --- a/tests/modules/testExecTime.cpp +++ b/tests/modules/testExecTime.cpp @@ -27,7 +27,6 @@ using namespace corsika::process::example_processors; TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); int tmp = 0; diff --git a/tests/modules/testNullModel.cpp b/tests/modules/testNullModel.cpp index a24433b080b4272955f01dfb2d116b43e0d5a2d6..6485b7df95a804e11301c0bc01095d705586fdea 100644 --- a/tests/modules/testNullModel.cpp +++ b/tests/modules/testNullModel.cpp @@ -14,7 +14,6 @@ using namespace corsika; TEST_CASE("NullModel", "[processes]") { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); logging::set_level(logging::level::trace); SECTION("interface") { diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp index 90a12288bd0b24c4c3021adaf66630fd1265d688..62983a2e08f49eef546e017b104cb104ead3de41 100644 --- a/tests/modules/testObservationPlane.cpp +++ b/tests/modules/testObservationPlane.cpp @@ -17,16 +17,17 @@ #include <corsika/framework/core/ParticleProperties.hpp> #include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/output/NoOutput.hpp> + #include <SetupTestEnvironment.hpp> #include <SetupTestStack.hpp> #include <SetupTestTrajectory.hpp> using namespace corsika; -TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") { +TEST_CASE("ObservationPlane", "interface") { logging::set_level(logging::level::trace); - corsika_logger->set_pattern("[%n:%^%-8l%$]: %v"); auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen); auto const& cs = *csPtr; @@ -54,23 +55,38 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") { SECTION("horizontal plane") { Plane const obsPlane(Point(cs, {10_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.})); - ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}), "particles.dat", - true); + ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.})); LengthType const length = obs.getMaxStepLength(particle, no_used_track); ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true); CHECK(length / 10_m == Approx(1).margin(1e-4)); CHECK(ret == ProcessReturn::ParticleAbsorbed); + + // particle does not reach plane: + { + setup::Trajectory no_hit_track = + setup::testing::make_track<setup::Trajectory>(line, 1_nm / constants::c); + LengthType const no_hit = obs.getMaxStepLength(particle, no_hit_track); + CHECK(no_hit == std::numeric_limits<double>::infinity() * 1_m); + } + + // particle past plane: + { + particle.setPosition({cs, {0_m, 0_m, -1_m}}); + setup::Trajectory no_hit_track = + setup::testing::make_track<setup::Trajectory>(line, 1_nm / constants::c); + LengthType const no_hit = obs.getMaxStepLength(particle, no_hit_track); + CHECK(no_hit == std::numeric_limits<double>::infinity() * 1_m); + } } SECTION("transparent plane") { Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.})); - ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 0., 1.}), "particles.dat", - false); + ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 0., 1.}), false); LengthType const length = obs.getMaxStepLength(particle, no_used_track); - ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true); + ProcessReturn const ret = obs.doContinuous(particle, no_used_track, false); CHECK(length / 1_m == Approx(1).margin(1e-4)); CHECK(ret == ProcessReturn::Ok); @@ -85,8 +101,7 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") { Plane const obsPlane(Point(cs, {10_m, 5_m, 5_m}), DirectionVector(cs, {1, 0.1, -0.05})); - ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}), "particles.dat", - true); + ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.})); LengthType const length = obs.getMaxStepLength(particle, no_used_track); ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true); @@ -94,4 +109,11 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") { CHECK(length / 10_m == Approx(1.1375).margin(1e-4)); CHECK(ret == ProcessReturn::ParticleAbsorbed); } + + SECTION("output") { + Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.})); + ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 0., 1.}), false); + auto const cfg = obs.getConfig(); + CHECK(cfg["type"]); + } } diff --git a/tests/modules/testOnShellCheck.cpp b/tests/modules/testOnShellCheck.cpp index ce1cb0b2d10175ef6e696177599265d9b8a70ef0..799049f29d94dd56d1cd1518bbf40029ead4a818 100644 --- a/tests/modules/testOnShellCheck.cpp +++ b/tests/modules/testOnShellCheck.cpp @@ -25,7 +25,6 @@ using namespace corsika; TEST_CASE("OnShellCheck", "[processes]") { logging::set_level(logging::level::debug); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); feenableexcept(FE_INVALID); using EnvType = setup::Environment; diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp index 62711c3f1da201a7ce5048520e1f32ba41d8d815..ab47527d80a254975e2585d6a8a356ebf7a1a420 100644 --- a/tests/modules/testParticleCut.cpp +++ b/tests/modules/testParticleCut.cpp @@ -26,7 +26,6 @@ using namespace corsika; TEST_CASE("ParticleCut", "processes") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] %v"); feenableexcept(FE_INVALID); using EnvType = setup::Environment; diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp index 8e178affc3f7c26c5a34df8caf2777386dba8d24..a4a93126a6c7456608f0cbd3f126e72910ae876d 100644 --- a/tests/modules/testPythia8.cpp +++ b/tests/modules/testPythia8.cpp @@ -23,7 +23,6 @@ using namespace corsika; TEST_CASE("Pythia", "[processes]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("linking pythia") { using namespace Pythia8; @@ -95,7 +94,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) { TEST_CASE("PythiaInterface", "[processes]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton); auto const& cs = *csPtr; diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp index 4a604fd8bbf32f440db2eb0e073aed684b412fe4..3d68eb0cd4d14cd07c09bfd9bbc8afc0a2116350 100644 --- a/tests/modules/testQGSJetII.cpp +++ b/tests/modules/testQGSJetII.cpp @@ -49,7 +49,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) { TEST_CASE("CORSIKA_DATA", "[processes]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("check CORSIKA_DATA") { @@ -68,7 +67,6 @@ TEST_CASE("CORSIKA_DATA", "[processes]") { TEST_CASE("QgsjetII", "[processes]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("Corsika -> QgsjetII") { CHECK(corsika::qgsjetII::convertToQgsjetII(PiMinus::code) == @@ -133,7 +131,6 @@ TEST_CASE("QgsjetII", "[processes]") { TEST_CASE("QgsjetIIInterface", "[processes]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen); [[maybe_unused]] auto const& env_dummy = env; @@ -184,7 +181,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { corsika::qgsjetII::Interaction model; model.doInteraction(view); // this also should produce some fragments - CHECK(view.getSize() == Approx(188).margin(2)); // this is not physics validation + CHECK(view.getSize() == Approx(228).margin(100)); // this is not physics validation int countFragments = 0; for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); } CHECK(countFragments == Approx(2).margin(1)); // this is not physics validation @@ -219,6 +216,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { auto [stackPtr, secViewPtr] = setup::testing::setup_stack( Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr); + [[maybe_unused]] setup::StackView& view = *(secViewPtr.get()); auto particle = stackPtr->first(); corsika::qgsjetII::Interaction model; GrammageType const length = model.getInteractionLength(particle); @@ -228,7 +226,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { auto [stackPtr, secViewPtr] = setup::testing::setup_stack( Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr); - setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] auto particle = stackPtr->first(); corsika::qgsjetII::Interaction model; model.doInteraction(view); CHECK(view.getSize() == Approx(18).margin(2)); // this is not physics validation @@ -237,28 +236,31 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { auto [stackPtr, secViewPtr] = setup::testing::setup_stack( Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr); - setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] auto particle = stackPtr->first(); corsika::qgsjetII::Interaction model; model.doInteraction(view); - CHECK(view.getSize() == Approx(7).margin(2)); // this is not physics validation + CHECK(view.getSize() == Approx(12).margin(8)); // this is not physics validation } { // Lambda is internally converted into neutron auto [stackPtr, secViewPtr] = setup::testing::setup_stack( Code::Lambda0, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr); - setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] auto particle = stackPtr->first(); corsika::qgsjetII::Interaction model; model.doInteraction(view); - CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation + CHECK(view.getSize() == Approx(15).margin(10)); // this is not physics validation } { // AntiLambda is internally converted into anti neutron auto [stackPtr, secViewPtr] = setup::testing::setup_stack( - Code::Lambda0Bar, 0, 0, 100_GeV, + Code::Lambda0Bar, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr); - setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] setup::StackView& view = *(secViewPtr.get()); + [[maybe_unused]] auto particle = stackPtr->first(); corsika::qgsjetII::Interaction model; model.doInteraction(view); - CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation + CHECK(view.getSize() == Approx(40).margin(20)); // this is not physics validation } } } diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp index dc237b4a379613dff02be6802ca7047776aa45d4..dc5249a7218e2e2dea96b7a23f791833ef4e8c9a 100644 --- a/tests/modules/testSibyll.cpp +++ b/tests/modules/testSibyll.cpp @@ -32,7 +32,6 @@ using namespace corsika::sibyll; TEST_CASE("Sibyll", "[processes]") { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); logging::set_level(logging::level::trace); SECTION("Sibyll -> Corsika") { @@ -99,7 +98,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) { TEST_CASE("SibyllInterface", "[processes]") { - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); logging::set_level(logging::level::trace); auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen); diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp index 00187ed9addbf6076641540dd494beccbd85758b..0b2a3ab7b9a489df955f00814ce371dfc4a026e9 100644 --- a/tests/modules/testStackInspector.cpp +++ b/tests/modules/testStackInspector.cpp @@ -24,7 +24,6 @@ using namespace corsika; TEST_CASE("StackInspector", "[processes]") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); auto const& rootCS = get_root_CoordinateSystem(); Point const origin(rootCS, {0_m, 0_m, 0_m}); diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp index 2040bc7fdd6326496d0b8eb11edb4e420c0472fc..73539692be7a8d2c4c791afbd5012d986760819a 100644 --- a/tests/modules/testTracking.cpp +++ b/tests/modules/testTracking.cpp @@ -20,32 +20,52 @@ #include <catch2/catch.hpp> +#include <boost/type_index.hpp> + using namespace corsika; +struct NonExistingDummyObject : public IVolume { + NonExistingDummyObject const& getVolume() const { return *this; } + bool contains(Point const&) const { return false; } +}; + template <typename T> int sgn(T val) { return (T(0) < val) - (val < T(0)); } /** - \file testTracking.cpp + @file testTracking.cpp - This is the unified and commond unit test for all Tracking algorithms: + This is the unified and common unit test for all Tracking algorithms: - tracking_leapfrog_curved::Tracking - tracking_leapfrog_straight::Tracking - tracking_line::Tracking + + The main part of tests are to inject particles at 10GeV momentum at + (-Rg,0,0) in +x direction into a sphere of radius Rg, where Rg is + the gyroradius (or 10m for neutral particles). Then it is checked + where the particles leaves the sphere for different charges + (-1,0,+1) and field strength (-50uT, 0T, +50uT). + + Each test is perfromed once, with the particles starting logically + outside of the Rg sphere (thus it first has to enter insides) and a + second time with the particle already logically inside the sphere. + + There is a second smaller, internal sphere at +z displacement. Only + neutral particles are allowed and expected to hit this. + + All those tests are parameterized, thus, they can be easily extended + or applied to new algorithms. + */ -TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking", - tracking_leapfrog_curved::Tracking, +TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking, tracking_leapfrog_straight::Tracking, tracking_line::Tracking) { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); - - logging::set_level(logging::level::trace); const HEPEnergyType P0 = 10_GeV; @@ -69,17 +89,19 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking", Bfield / 1_uT, outer)) { CORSIKA_LOG_DEBUG( - "********************\n TEST section PID={}, Bfield={} " - "uT, outer (?)={}", - PID, Bfield / 1_uT, outer); + "********************\n TEST algo={} section PID={}, " + "Bfield={} " + "uT, start_outside={}", + boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT, outer); const int chargeNumber = get_charge_number(PID); LengthType radius = 10_m; int deflect = 0; if (chargeNumber != 0 and Bfield != 0_T) { deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection - LengthType const gyroradius = - P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV); + LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) * + constants::c / (abs(get_charge(PID)) * abs(Bfield))); + CORSIKA_LOG_DEBUG("Rg={} deflect={}", gyroradius, deflect); radius = gyroradius; } @@ -91,6 +113,12 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking", TestType tracking; Point const center(cs, {0_m, 0_m, 0_m}); auto target = setup::Environment::createNode<Sphere>(center, radius); + // every particle should hit target_2 + auto target_2 = setup::Environment::createNode<Sphere>( + Point(cs, {-radius * 3 / 4, 0_m, 0_m}), radius * 0.2); + // only neutral particles hit_target_neutral + auto target_neutral = setup::Environment::createNode<Sphere>( + Point(cs, {radius / 2, 0_m, 0_m}), radius * 0.1); using MyHomogeneousModel = MediumPropertyModel< UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>; @@ -99,8 +127,18 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking", target->setModelProperties<MyHomogeneousModel>( Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m), NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.})); + target_neutral->setModelProperties<MyHomogeneousModel>( + Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m), + NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.})); + target_2->setModelProperties<MyHomogeneousModel>( + Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m), + NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.})); auto* targetPtr = target.get(); + auto* targetPtr_2 = target_2.get(); + auto* targetPtr_neutral = target_neutral.get(); worldPtr->addChild(std::move(target)); + targetPtr->addChild(std::move(target_2)); + targetPtr->addChild(std::move(target_neutral)); auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs); { [[maybe_unused]] auto& viewPtr_dum = viewPtr; } @@ -123,18 +161,29 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking", particle.setMomentum(traj.getDirection(1) * particle.getMomentum().getNorm()); if (outer) { // now we know we are in target volume, depending on "outer" - CHECK(traj.getLength(1) == 0_m); + CHECK(traj.getLength(1) / 1_m == Approx(0).margin(1e-3)); CHECK(nextVol == targetPtr); } // move forward, until we leave target volume - while (nextVol == targetPtr) { + bool hit_neutral = false; + bool hit_2nd = false; + while (nextVol != worldPtr) { + if (nextVol == targetPtr_neutral) hit_neutral = true; + if (nextVol == targetPtr_2) hit_2nd = true; const auto [traj2, nextVol2] = tracking.getTrack(particle); nextVol = nextVol2; particle.setNode(nextVol); particle.setPosition(traj2.getPosition(1)); particle.setMomentum(traj2.getDirection(1) * particle.getMomentum().getNorm()); + CORSIKA_LOG_TRACE("pos={}, p={}, |p|={} |v|={}, delta-l={}, delta-t={}", + particle.getPosition(), particle.getMomentum(), + particle.getMomentum().getNorm(), + particle.getVelocity().getNorm(), traj2.getLength(1), + traj2.getLength(1) / particle.getVelocity().getNorm()); } CHECK(nextVol == worldPtr); + CHECK(hit_2nd == true); + CHECK(hit_neutral == (deflect == 0 ? true : false)); Point pointCheck(cs, (deflect == 0 ? radius : 0_m), (deflect * radius), 0_m); @@ -147,3 +196,166 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking", Approx(0).margin(1e-3)); } } + +/** specifc test for curved leap-frog algorithm **/ + +TEST_CASE("TrackingLeapFrogCurved") { + + logging::set_level(logging::level::info); + + const HEPEnergyType P0 = 10_GeV; + + corsika::Code PID = Code::MuPlus; + + using MyHomogeneousModel = MediumPropertyModel< + UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>; + + SECTION("infinite sphere / universe") { + + auto [env, csPtr, worldPtr] = + corsika::setup::testing::setup_environment(Code::Oxygen, 100_uT); + { [[maybe_unused]] const auto& env_dummy = env; } + auto const& cs = *csPtr; + + tracking_leapfrog_curved::Tracking tracking; + Point const center(cs, {0_m, 0_m, 0_m}); + + auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs); + { [[maybe_unused]] auto& viewPtr_dum = viewPtr; } + auto particle = stack->first(); + // Note: momentum in X-direction + // magnetic field in X-direction + + particle.setPosition(Point(cs, 0_m, 0_m, 0_m)); + + Sphere sphere(center, 1_km * std::numeric_limits<double>::infinity()); + + auto intersections = tracking.intersect(particle, sphere); + // this must be a "linear trajectory" with no curvature + + CHECK(intersections.hasIntersections() == false); + } + + SECTION("momentum along field") { + + auto [env, csPtr, worldPtr] = + corsika::setup::testing::setup_environment(Code::Oxygen, 100_uT); + { [[maybe_unused]] const auto& env_dummy = env; } + auto const& cs = *csPtr; + + tracking_leapfrog_curved::Tracking tracking; + Point const center(cs, {0_m, 0_m, 0_m}); + auto target = setup::Environment::createNode<Sphere>(center, 10_km); + + MagneticFieldVector magneticfield(cs, 100_T, 0_T, 0_uT); + target->setModelProperties<MyHomogeneousModel>( + Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m), + NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.})); + auto* targetPtr = target.get(); + worldPtr->addChild(std::move(target)); + + auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs); + { [[maybe_unused]] auto& viewPtr_dum = viewPtr; } // prevent warning + auto particle = stack->first(); + // Note: momentum in X-direction + // magnetic field in X-direction + + particle.setNode(targetPtr); // set particle outside "target" volume + particle.setPosition(Point(cs, 0_m, 0_m, 0_m)); + + auto [traj, nextVol] = tracking.getTrack(particle); + { [[maybe_unused]] auto const& dummy = nextVol; } // prevent warning + + // this must be a "linear trajectory" with no curvature + CHECK(traj.getDirection(0).getComponents() == traj.getDirection(1).getComponents()); + } +} + +TEMPLATE_TEST_CASE("TrackingFail", "doesntwork", tracking_leapfrog_curved::Tracking, + tracking_leapfrog_straight::Tracking, tracking_line::Tracking) { + + logging::set_level(logging::level::info); + + const HEPEnergyType P0 = 10_GeV; + + auto [env, csPtr, worldPtr] = + corsika::setup::testing::setup_environment(Code::Oxygen, 1_mT); + { [[maybe_unused]] const auto& env_dummy = env; } + auto const& cs = *csPtr; + + TestType tracking; + Point const center(cs, {0_m, 0_m, 0_m}); + + auto [stack, viewPtr] = + setup::testing::setup_stack(Code::Proton, 0, 0, P0, worldPtr, cs); + { [[maybe_unused]] auto& viewPtr_dum = viewPtr; } + auto particle = stack->first(); + NonExistingDummyObject const dummy; + CHECK_THROWS(tracking.intersect(particle, dummy)); +} + +TEMPLATE_TEST_CASE("TrackingPlane", "plane", tracking_leapfrog_curved::Tracking, + tracking_leapfrog_straight::Tracking, tracking_line::Tracking) { + + logging::set_level(logging::level::trace); + + const HEPEnergyType P0 = 10_GeV; + + auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Gamma); + // for algorithms that know magnetic deflections choose: +-50uT, 0uT + // otherwise just 0uT + auto Bfield = GENERATE(filter( + []([[maybe_unused]] MagneticFluxType v) { + if constexpr (std::is_same_v<TestType, tracking_line::Tracking>) + return v == 0_uT; + else + return true; + }, + values<MagneticFluxType>({50_uT, 0_uT, -50_uT}))); + + SECTION(fmt::format("Tracking PID={}, Bfield={} uT", PID, Bfield / 1_uT)) { + + CORSIKA_LOG_DEBUG( + "********************\n TEST algo={} section PID={}, " + "Bfield={}uT ", + boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT); + + const int chargeNumber = get_charge_number(PID); + LengthType radius = 10_m; + int deflect = 0; + if (chargeNumber != 0 and Bfield != 0_T) { + deflect = 1; + LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) * + constants::c / (abs(get_charge(PID)) * abs(Bfield))); + CORSIKA_LOG_DEBUG("Rg={} deflect={}", gyroradius, deflect); + radius = gyroradius; + } + + auto [env, csPtr, worldPtr] = + corsika::setup::testing::setup_environment(Code::Oxygen, Bfield); + { [[maybe_unused]] const auto& env_dummy = env; } + auto const& cs = *csPtr; + + TestType tracking; + Point const center(cs, {0_m, 0_m, 0_m}); + + auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs); + { [[maybe_unused]] auto& viewPtr_dum = viewPtr; } + auto particle = stack->first(); + // Note: momentum in X-direction + // magnetic field in Z-direction + + particle.setPosition(Point(cs, -radius, 0_m, 0_m)); + + // put plane where we expect deflection by radius/2 + Plane const plane(Point(cs, radius * (1 - sqrt(3. / 4.)), 0_m, 0_m), + DirectionVector(cs, {-1., 0., 0.})); + Intersections const hit = tracking.intersect(particle, plane); + + CORSIKA_LOG_DEBUG("entry={} exit={}", hit.getEntry(), hit.getExit()); + + CHECK(hit.hasIntersections()); + CHECK(hit.getEntry() / 1_s == Approx(0.00275 * deflect).margin(0.0003)); + CHECK(hit.getExit() == 1_s * std::numeric_limits<double>::infinity()); + } +} diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp index 0f9f4a86243ec0f7f5b0c16e989b9459a5e0a293..f55aa032948c22026b170776d8946dc074c6e688 100644 --- a/tests/modules/testUrQMD.cpp +++ b/tests/modules/testUrQMD.cpp @@ -59,7 +59,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) { TEST_CASE("UrQMD") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); SECTION("conversion") { CHECK_THROWS(corsika::urqmd::convertFromUrQMD(106, 0)); @@ -99,6 +98,7 @@ TEST_CASE("UrQMD") { { [[maybe_unused]] auto const& env_dummy = env; } auto [stack, view] = setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs); + [[maybe_unused]] setup::StackView& viewRef = *(view.get()); CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) == Approx(105).margin(5)); } @@ -109,6 +109,7 @@ TEST_CASE("UrQMD") { { [[maybe_unused]] auto const& env_dummy = env; } auto [stack, view] = setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs); + [[maybe_unused]] setup::StackView& viewRef = *(view.get()); CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle())); } @@ -121,6 +122,7 @@ TEST_CASE("UrQMD") { auto [stackPtr, secViewPtr] = setup::testing::setup_stack( Code::Nucleus, A, Z, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr); + [[maybe_unused]] setup::StackView& viewRef = *(secViewPtr.get()); CHECK(stackPtr->getEntries() == 1); CHECK(secViewPtr->getEntries() == 0); @@ -172,6 +174,7 @@ TEST_CASE("UrQMD") { auto [stackPtr, secViewPtr] = setup::testing::setup_stack( Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr); + [[maybe_unused]] auto particle = stackPtr->first(); CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target } @@ -227,22 +230,4 @@ TEST_CASE("UrQMD") { CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() == Approx(0).margin(1e-2)); } - - SECTION("K0Short projectile") { - auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon); - [[maybe_unused]] auto const& env_dummy = env; // against warnings - [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings - - auto [stackPtr, secViewPtr] = setup::testing::setup_stack( - Code::K0Short, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, - *csPtr); - CHECK(stackPtr->getEntries() == 1); - CHECK(secViewPtr->getEntries() == 0); - - // must be assigned to variable, cannot be used as rvalue?! - auto projectile = secViewPtr->getProjectile(); - auto const projectileMomentum = projectile.getMomentum(); - - CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); - } } diff --git a/tests/output/CMakeLists.txt b/tests/output/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..98a43b4c218bf8d9e21bb7b9e56ac1275cbc9b0e --- /dev/null +++ b/tests/output/CMakeLists.txt @@ -0,0 +1,16 @@ +set (test_output_sources + TestMain.cpp + testOutputManager.cpp + testDummyOutputManager.cpp + testParquetStreamer.cpp + testWriterObservationPlane.cpp + testWriterTracks.cpp + ) + +CORSIKA_ADD_TEST (testOutput SOURCES ${test_output_sources}) + +target_compile_definitions ( + testOutput + PRIVATE + REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}" +) diff --git a/tests/output/TestMain.cpp b/tests/output/TestMain.cpp new file mode 100644 index 0000000000000000000000000000000000000000..51532584b8e03b35d79301543ac8f80b598ba544 --- /dev/null +++ b/tests/output/TestMain.cpp @@ -0,0 +1,11 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one + // cpp file +#include <catch2/catch.hpp> diff --git a/tests/output/testDummyOutputManager.cpp b/tests/output/testDummyOutputManager.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cc571832f601c1e9814ca25e1de9d4db152d1afc --- /dev/null +++ b/tests/output/testDummyOutputManager.cpp @@ -0,0 +1,32 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <corsika/output/DummyOutputManager.hpp> +#include <corsika/framework/core/Logging.hpp> + +using namespace corsika; + +class DummyOutput {}; + +TEST_CASE("DummyOutputManager", "interface") { + + logging::set_level(logging::level::info); + + // output manager performs nothing, no action, just interface + DummyOutputManager output; + + DummyOutput test; + output.add("test", test); + + output.startOfLibrary(); + output.startOfShower(); + output.endOfShower(); + output.endOfLibrary(); +} diff --git a/tests/output/testOutputManager.cpp b/tests/output/testOutputManager.cpp new file mode 100644 index 0000000000000000000000000000000000000000..58f90589c2206ed1e01607cd216a6b76c53a20c9 --- /dev/null +++ b/tests/output/testOutputManager.cpp @@ -0,0 +1,184 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <boost/filesystem.hpp> + +#include <corsika/framework/core/Logging.hpp> + +#include <corsika/output/OutputManager.hpp> +#include <corsika/output/NoOutput.hpp> + +using namespace corsika; + +struct DummyNoOutput : public NoOutput { + void check() { + NoOutput::startOfLibrary("./"); + NoOutput::startOfShower(); + NoOutput::endOfShower(); + NoOutput::endOfLibrary(); + NoOutput::getConfig(); + NoOutput::getSummary(); + } + void checkWrite() { NoOutput::write(Code::Unknown, 1_eV, 1_m, 1_m); } +}; + +struct DummyOutput : public BaseOutput { + + mutable bool isConfig_ = false; + mutable bool isSummary_ = false; + bool startLibrary_ = false; + bool startShower_ = false; + bool endLibrary_ = false; + bool endShower_ = false; + + void startOfLibrary(boost::filesystem::path const&) { startLibrary_ = true; } + + void startOfShower() { + BaseOutput::startOfShower(); + startShower_ = true; + } + + void endOfShower() { endShower_ = true; } + + void endOfLibrary() { endLibrary_ = true; } + + YAML::Node getConfig() const { + isConfig_ = true; + return YAML::Node(); + } + + YAML::Node getSummary() { + isSummary_ = true; + return BaseOutput::getSummary(); + } +}; + +TEST_CASE("OutputManager") { + + logging::set_level(logging::level::info); + + SECTION("standard") { + + // preparation + if (boost::filesystem::exists("./out_test")) { + boost::filesystem::remove_all("./out_test"); + } + + // output manager performs nothing, no action, just interface + OutputManager output("check", "./out_test"); + + CHECK(boost::filesystem::is_directory("./out_test/check")); + + DummyOutput test; + output.add("test", test); + + CHECK_THROWS(output.add( + "test", + test)); // should emit warning which cannot be catched, but no action or failure + + CHECK(test.isConfig_); + test.isConfig_ = false; + + output.startOfLibrary(); + CHECK(test.startLibrary_); + test.startLibrary_ = false; + + output.startOfShower(); + CHECK(test.startShower_); + test.startShower_ = false; + + output.endOfShower(); + CHECK(test.endShower_); + test.endShower_ = false; + + output.endOfLibrary(); + CHECK(test.endLibrary_); + CHECK(test.isSummary_); + test.isSummary_ = false; + test.endLibrary_ = false; + + CHECK(boost::filesystem::exists("./out_test/check/test/summary.yaml")); + } + + SECTION("auto-write") { + + // preparation + if (boost::filesystem::exists("./out_test")) { + boost::filesystem::remove_all("./out_test"); + } + + // output manager performs nothing, no action, just interface + OutputManager* output = new OutputManager("check", "./out_test"); + + CHECK(boost::filesystem::is_directory("./out_test/check")); + + DummyOutput test; + output->add("test", test); + output->startOfLibrary(); + output->startOfShower(); + + // check support for closing automatically + delete output; + output = 0; + + CHECK(boost::filesystem::exists("./out_test/check/test/summary.yaml")); + } + + SECTION("failures") { + + logging::set_level(logging::level::info); + + // preparation + if (boost::filesystem::exists("./out_test")) { + boost::filesystem::remove_all("./out_test"); + } + + // output manager performs nothing, no action, just interface + OutputManager output("check", "./out_test"); + CHECK_THROWS(new OutputManager("check", "./out_test")); + + // CHECK_THROWS(output.startOfShower()); + // CHECK_THROWS(output.endOfShower()); + CHECK_THROWS(output.endOfLibrary()); + + output.startOfLibrary(); + + CHECK_THROWS(output.startOfLibrary()); + // CHECK_THROWS(output.endOfShower()); + // CHECK_THROWS(output.endOfLibrary()); + + output.startOfShower(); + + CHECK_THROWS(output.startOfLibrary()); + // CHECK_THROWS(output.startOfShower()); + // CHECK_THROWS(output.endOfLibrary()); + + output.endOfShower(); + + CHECK_THROWS(output.startOfLibrary()); + // CHECK_THROWS(output.startOfShower()); + // CHECK_THROWS(output.endOfShower()); + + output.endOfLibrary(); + + // CHECK_THROWS(output.endOfShower()); + // CHECK_THROWS(output.startOfShower()); + // CHECK_THROWS(output.endOfLibrary()); + } + + SECTION("NoOutput") { + // this is one of the classes where testing is a bit useless, but we can at least make + // sure the interface exists. + DummyNoOutput nothing; + + nothing.check(); + nothing.checkWrite(); + } +} \ No newline at end of file diff --git a/tests/output/testParquetStreamer.cpp b/tests/output/testParquetStreamer.cpp new file mode 100644 index 0000000000000000000000000000000000000000..387e5cade02d9e9cf1f1c1532cc1594bffb65d4a --- /dev/null +++ b/tests/output/testParquetStreamer.cpp @@ -0,0 +1,56 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <boost/filesystem.hpp> + +#include <corsika/output/ParquetStreamer.hpp> +#include <corsika/framework/core/Logging.hpp> + +using namespace corsika; + +TEST_CASE("ParquetStreamer") { + + logging::set_level(logging::level::info); + + SECTION("standard") { + + // preparation + if (boost::filesystem::exists("./parquet_test.parquet")) { + boost::filesystem::remove_all("./parquet_test.parquet"); + } + + ParquetStreamer test; + CHECK_FALSE(test.isInit()); + CHECK_THROWS(test.getWriter()); + test.initStreamer("./parquet_test.parquet"); + + test.addField("testint", parquet::Repetition::REQUIRED, parquet::Type::INT32, + parquet::ConvertedType::INT_32); + test.addField("testfloat", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + + test.enableCompression(1); + + test.buildStreamer(); + CHECK(test.isInit()); + + int testint = 1; + double testfloat = 2.0; + + std::shared_ptr<parquet::StreamWriter> writer = test.getWriter(); + (*writer) << static_cast<int>(testint) << static_cast<int>(testint) + << static_cast<float>(testfloat) << parquet::EndRow; + + test.closeStreamer(); + CHECK_THROWS(test.getWriter()); + CHECK_FALSE(test.isInit()); + CHECK(boost::filesystem::exists("./parquet_test.parquet")); + } +} \ No newline at end of file diff --git a/tests/output/testWriterObservationPlane.cpp b/tests/output/testWriterObservationPlane.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4b780a2b4bbcde53d531f4a6a113d4b2a7087f3a --- /dev/null +++ b/tests/output/testWriterObservationPlane.cpp @@ -0,0 +1,50 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <boost/filesystem.hpp> + +#include <corsika/modules/writers/ObservationPlaneWriterParquet.hpp> + +#include <corsika/framework/core/Logging.hpp> +#include <corsika/framework/geometry/QuantityVector.hpp> + +using namespace corsika; + +struct TestWriterPlane : public ObservationPlaneWriterParquet { + + YAML::Node getConfig() const { return YAML::Node(); } + + void checkWrite() { + ObservationPlaneWriterParquet::write(Code::Unknown, 1_eV, 2_m, 3_m); + } +}; + +TEST_CASE("ObservationPlaneWriterParquet") { + + logging::set_level(logging::level::info); + + SECTION("standard") { + + // preparation + if (boost::filesystem::exists("./output_dir")) { + boost::filesystem::remove_all("./output_dir"); + } + boost::filesystem::create_directory("./output_dir"); + + TestWriterPlane test; + test.startOfLibrary("./output_dir"); + test.startOfShower(); + test.checkWrite(); + test.endOfShower(); + test.endOfLibrary(); + + CHECK(boost::filesystem::exists("./output_dir/particles.parquet")); + } +} diff --git a/tests/output/testWriterTracks.cpp b/tests/output/testWriterTracks.cpp new file mode 100644 index 0000000000000000000000000000000000000000..38cf0cf5985a956aaef03b589b93a1b6b9f9b260 --- /dev/null +++ b/tests/output/testWriterTracks.cpp @@ -0,0 +1,49 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <boost/filesystem.hpp> + +#include <corsika/modules/writers/TrackWriterParquet.hpp> + +#include <corsika/framework/core/Logging.hpp> + +using namespace corsika; + +struct TestWriterTrack : public TrackWriterParquet { + + YAML::Node getConfig() const { return YAML::Node(); } + + void checkWrite() { + TrackWriterParquet::write(Code::Unknown, 1_eV, {2_m, 3_m, 4_m}, {5_m, 6_m, 7_m}); + } +}; + +TEST_CASE("TrackWriterParquet") { + + logging::set_level(logging::level::info); + + SECTION("standard") { + + // preparation + if (boost::filesystem::exists("./output_dir_tracks")) { + boost::filesystem::remove_all("./output_dir_tracks"); + } + boost::filesystem::create_directory("./output_dir_tracks"); + + TestWriterTrack test; + test.startOfLibrary("./output_dir_tracks"); + test.startOfShower(); + test.checkWrite(); + test.endOfShower(); + test.endOfLibrary(); + + CHECK(boost::filesystem::exists("./output_dir_tracks/tracks.parquet")); + } +}