diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 7ae939f157b2250f7803ba14098957c6f4eae8bd..eea04e67c9192d1de66fd7ee3af456706abfdda5 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -480,7 +480,7 @@ coverage:
- ctest -j4
- cmake --build . --target coverage
- tar czf coverage-report.tar.gz coverage-report
- coverage: '/^.*functions\.+:\s(.*\%)\s/'
+ coverage: '/^.*lines\.+:\s(.*\%)\s/'
rules:
- if: '$CI_MERGE_REQUEST_LABELS =~ /Ready for code review/' # run on merge requests, if label 'Ready for code review' is set
- if: $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a4d6cd6c9d2a360cc5561085f010f12c9e5a8faf..a2b05773a97363db09577e32df92c91b64c5f272 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -165,7 +165,11 @@ if (CMAKE_BUILD_TYPE STREQUAL Coverage)
add_custom_command (
OUTPUT coverage.info
COMMAND ${LCOV_BIN_DIR}/lcov -q --remove raw-coverage.info "*/usr/*" "/usr/*" --output-file coverage2.info
- COMMAND ${LCOV_BIN_DIR}/lcov --remove coverage2.info "*/externals/*" "*/tests/*" "*/sibyll2.3d.cpp" "*/.conan/*" "*/include/Pythia8/*" "${CMAKE_SOURCE_DIR}/modules/*" "${CMAKE_BINARY_DIR}/modules/*" --output-file coverage.info
+ COMMAND ${LCOV_BIN_DIR}/lcov --remove coverage2.info
+ "*/externals/*" "*/tests/*" "*/sibyll2.3d.cpp" "*/.conan/*"
+ "*/include/Pythia8/*" "*/install/*" "${CMAKE_SOURCE_DIR}/modules/*"
+ "${CMAKE_BINARY_DIR}/modules/*"
+ --output-file coverage.info
COMMAND ${CMAKE_COMMAND} -E remove coverage2.info
DEPENDS raw-coverage.info
)
@@ -238,8 +242,8 @@ target_link_libraries (
CONAN_PKG::yaml-cpp
CONAN_PKG::arrow
cnpy # for SaveBoostHistogram
- stdc++fs # experimental::filesystem
)
+
# "src" is needed, since some headers (namely GeneratedParticleProperties.inc ec.) are produced by python script from e.g. ParticleData.xml
add_subdirectory (src)
add_subdirectory (documentation)
diff --git a/README.rst b/README.rst
new file mode 100644
index 0000000000000000000000000000000000000000..896a4e571c54f19ab09088c602d960f3a416ef00
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,201 @@
+CORSIKA 8 Framework for Particle Cascades in Astroparticle Physics
+===================================================================
+
+The purpose of CORSIKA is to simulate any particle cascades in
+astroparticle physics or astrophysical context. A lot of emphasis is
+put on modularity, flexibility, completeness, validation and
+correctness. To boost computational efficiency different techniques
+are provided, like thinning or cascade equations. The aim is that
+CORSIKA remains the most comprehensive framework for simulating
+particle cascades with stochastic and continuous processes.
+
+The software makes extensive use of static design patterns and
+compiler optimization. Thus, the most fundamental configuration
+decision of the user must be performed at compile time. At run time
+model parameters can still be changed.
+
+CORSIKA 8 is by default released under the GPLv3 license. See `license
+file <https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/LICENSE>`_
+which is part of every release and the source code.
+
+If you use, or want to refer to, CORSIKA 8 please cite `"Towards a Next
+Generation of CORSIKA: A Framework for the Simulation of Particle
+Cascades in Astroparticle Physics", Comput.Softw.Big Sci. 3 (2019)
+2 <https://doi.org/10.1007/s41781-018-0013-0>`_. We kindly ask (and
+require) any relevant improvement or addition to be offered or
+contributed to the main CORSIKA 8 repository for the benefit of the
+whole community.
+
+When you plan to contribute to CORSIKA 8 check the guidelines outlined here:
+`coding
+guidelines <https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/CONTRIBUTING.md>`_. Code
+that fails the review by the CORSIKA author group must be improved
+before it can be merged in the official code base. After your code has
+been accepted and merged, you become a contributor of the CORSIKA 8
+project (code author).
+
+IMPORTANT: Before you contribute, you need to read and agree to the
+`collaboration agreement
+<https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/COLLABORATION_AGREEMENT.md>`_. The
+agreement can be discussed, and eventually improved.
+
+We also want to point you to the `MCnet guidelines
+<https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/MCNET_GUIDELINES>`_,
+which are very useful also for us.
+
+
+Get in contact
+--------------
+ * Connect to https://gitlab.ikp.kit.edu register yourself and join the "Air Shower Physics" group. Write to me (ralf.ulrich@kit.edu) only in case there are problems with that.
+ * Connect to corsika-devel@lists.kit.edu (self-register at
+ https://www.lists.kit.edu/sympa/subscribe/corsika-devel) to get in
+ touch with the project.
+
+ * Register on the corsika slack channel.
+
+
+Installation
+------------
+CORSIKA 8 is tested regularly at least on gcc7.3.0 and clang-8.0.0.
+
+Prerequisites
+~~~~~~~~~~~~~
+You will also need:
+ * Python 3 (supported versions are Python >= 3.6), with pip
+ * conan (via pip)
+ * cmake
+ * git
+ * g++, gfortran, binutils, make
+
+On a bare Ubuntu 20.04, just add:
+::
+
+ sudo apt-get install python3 python3-pip cmake g++ gfortran git doxygen graphviz
+
+
+On a bare CentOS 7 install python3, pip3 (pip from python3) and cmake3. Any of the devtools 7, 8, 9 should work (at least).
+Also initialize devtools, before building CORSIKA 8:
+::
+
+ source /opt/rh/devtoolset-9/enable
+
+
+CORSIKA 8 uses the `conan <https://conan.io/>`_ package manager to
+manage our dependencies. If you do not have Conan installed, it can be
+installed with:
+::
+
+ pip install --user conan
+
+
+Compiling
+~~~~~~~~~
+
+Once Conan is installed, follow these steps to download and install CORSIKA 8:
+::
+
+ git clone --recursive git@gitlab.ikp.kit.edu:AirShowerPhysics/corsika.git
+ mkdir corsika-build
+ cd corsika-build
+ cmake ../corsika -DCMAKE_INSTALL_PREFIX=../corsika-install
+ make -j8
+ make install
+
+
+Installation (using docker containers)
+--------------------------------------
+
+There are docker containers prepared that bring all the environment and packages you need to run CORSIKA. See `docker hub <https://hub.docker.com/repository/docker/corsika/devel>`_ for a complete overview.
+
+Prerequisites
+~~~~~~~~~~~~~
+
+You only need docker, e.g. on Ubuntu: :code:`sudo apt-get install docker` and of course root access.
+
+Compiling
+---------
+
+Follow these steps to download and install CORSIKA 8, master development version
+::
+
+ git clone --recursive https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika.git
+ sudo docker run -v $PWD:/corsika -it corsika/devel:clang-8 /bin/bash
+ mkdir build
+ cd build
+ cmake ../corsika -DCMAKE_INSTALL_PREFIX=../corsika-install
+ make -j8
+ make install
+
+Runing Unit Tests
+-----------------
+
+Note, before you run *any* executbale you must also define the
+:code:`CORSIKA_DATA` environment variable to point to the location where you
+cloned corsika :code:`modules/data`, thus typically
+::
+
+ export CORSIKA_DATA=$PWD/../corsika/modules/data
+
+To run the Unit Tests, just type :code:`ctest` in your build area.
+
+
+Running examples
+----------------
+
+To see how a relatively simple hadron cascade develops,
+see :code:`examples/cascade_example.cpp` for a starting point.
+
+To run the cascade_example, or any other CORSIKA 8 application, you
+must first compile it wrt. to the CORSIKA 8 header-only framework. This
+can be done best by copying
+e.g. :code:`corsika-install/share/corsika/examples/` to your working place
+(e.g. :code:`corsika-work`).
+
+Next, you need to define the environment variable :code:`corsika_DIR` to point to, either,
+your build, or your install area. Thus, e.g.
+::
+
+ export corsika_DIR=<dir where you installed CORSIKA 8 to, or where you buld it">
+
+Then compile your example/application with
+::
+
+ cd corsika-work
+ cmake .
+ make
+ bin/cascade_example
+
+Visualize output (needs gnuplot installed):
+::
+
+ bash $corsika_DIR/share/corsika/tools/plot_tracks.sh tracks.dat
+ firefox tracks.dat.gif
+
+(note, if you use the corsika_DIR to point to the build area: the script :code:`plot_tracks.sh` is
+not copied to the build area, it is only part of the source tree at :code:`tools/plot_tracks.sh`)
+
+Or also consider the :code:`vertical_EAS` example in the same directory,
+which can be configured with command line options.
+Run :code:`bin/vertical_EAS` to get basic help.
+
+
+Generating doxygen documentation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+To generate the documentation, you need doxygen and graphviz. If you work with
+the docker corsika/devel containers this is already included.
+Otherwise, e.g. on Ubuntu 18.04, do:
+::
+
+ sudo apt-get install doxygen graphviz
+
+Switch to the corsika build directory and do
+::
+
+ make doxygen
+ make install
+
+open with firefox:
+::
+
+ firefox ../corsika-install/share/corsika/doc/html/index.html
diff --git a/cmake/FindSphinx.cmake b/cmake/FindSphinx.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..9a6705c4840e57d3c7940289be5e498ad0fcce15
--- /dev/null
+++ b/cmake/FindSphinx.cmake
@@ -0,0 +1,11 @@
+# Look for an executable called sphinx-build
+find_program (SPHINX_EXECUTABLE
+ NAMES sphinx-build
+ DOC "Path to sphinx-build executable")
+
+include (FindPackageHandleStandardArgs)
+
+#Handle standard arguments to find_package like REQUIRED and QUIET
+find_package_handle_standard_args (Sphinx
+ "Failed to find sphinx-build executable"
+ SPHINX_EXECUTABLE)
diff --git a/cmake/corsikaConfig.cmake.in b/cmake/corsikaConfig.cmake.in
index 81a499c04bb861e6a660d90bee533b461b30f694..546d5141d2c467cf7293bddf4cabef6fa28d481f 100644
--- a/cmake/corsikaConfig.cmake.in
+++ b/cmake/corsikaConfig.cmake.in
@@ -25,6 +25,12 @@ if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
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/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index 2ca24ee5a4a9c080276731e9652f7b03202a82ec..a932e4de8b596217bde32b1bb672e69866f05dc6 100644
--- a/corsika/detail/framework/geometry/QuantityVector.inl
+++ b/corsika/detail/framework/geometry/QuantityVector.inl
@@ -128,7 +128,7 @@ namespace corsika {
}
template <typename TDimension>
- inline auto& QuantityVector<TDimension>::operator-() const {
+ inline auto QuantityVector<TDimension>::operator-() const {
return QuantityVector<TDimension>(-eigenVector_);
}
diff --git a/corsika/detail/framework/geometry/Vector.inl b/corsika/detail/framework/geometry/Vector.inl
index 352eb8b3095ddec2221c1af41dc8f71b499746e2..93907b944206f64364c254473739aa324f9c6aeb 100644
--- a/corsika/detail/framework/geometry/Vector.inl
+++ b/corsika/detail/framework/geometry/Vector.inl
@@ -194,7 +194,7 @@ namespace corsika {
}
template <typename TDimension>
- inline auto& Vector<TDimension>::operator-() const {
+ inline auto Vector<TDimension>::operator-() const {
return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
-BaseVector<TDimension>::getQuantityVector());
}
diff --git a/corsika/detail/framework/process/BoundaryCrossingProcess.hpp b/corsika/detail/framework/process/BoundaryCrossingProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1efd36454579d35b31d9f97421c42ca42a31674c
--- /dev/null
+++ b/corsika/detail/framework/process/BoundaryCrossingProcess.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/framework/process/ProcessTraits.hpp>
+
+namespace corsika {
+
+ /**
+ traits test for BoundaryCrossingProcess::doBoundaryCrossing method
+ */
+ template <class TProcess, typename TReturn, typename TParticle>
+ struct has_method_doBoundaryCrossing
+ : public detail::has_method_signature<TReturn, TParticle&,
+ typename TParticle::node_type const&,
+ typename TParticle::node_type const&> {
+
+ ///! method signature
+ using detail::has_method_signature<
+ TReturn, TParticle&, typename TParticle::node_type const&,
+ typename TParticle::node_type const&>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::template doBoundaryCrossing<TParticle>)) test(
+ std::nullptr_t);
+
+ //! non templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::doBoundaryCrossing)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! @file BoundaryCrossingProcess.hpp
+ //! value traits type
+
+ template <class TProcess, typename TReturn, typename TParticle>
+ bool constexpr has_method_doBoundaryCrossing_v =
+ has_method_doBoundaryCrossing<TProcess, TReturn, TParticle>::value;
+
+} // namespace corsika
diff --git a/corsika/detail/framework/process/ContinuousProcess.hpp b/corsika/detail/framework/process/ContinuousProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..451f45a68be3c592e3ac4d68e58f78cab7f682d3
--- /dev/null
+++ b/corsika/detail/framework/process/ContinuousProcess.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 <corsika/framework/process/ProcessTraits.hpp>
+
+namespace corsika {
+
+ /**
+ traits test for ContinuousProcess::doContinuous method
+ */
+ template <class TProcess, typename TReturn, typename TArg1, typename TArg2>
+ struct has_method_doContinuous
+ : public detail::has_method_signature<TReturn, TArg1, TArg2, bool> {
+
+ //! method signature
+ using detail::has_method_signature<TReturn, TArg1, TArg2, bool>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::template doContinuous<TArg1, TArg2>)) test(
+ std::nullptr_t);
+
+ //! non templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::doContinuous)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! @file ContinuousProcess.hpp
+ //! value traits type
+ template <class TProcess, typename TReturn, typename TArg1, typename TArg2>
+ bool constexpr has_method_doContinuous_v =
+ has_method_doContinuous<TProcess, TReturn, TArg1, TArg2>::value;
+
+ /**
+ traits test for ContinuousProcess::getMaxStepLength method
+ */
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ struct has_method_getMaxStepLength
+ : public detail::has_method_signature<TReturn, TArgs...> {
+
+ //! method signature
+ using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated option
+ template <class T>
+ static decltype(testSignature(&T::template getMaxStepLength<TArgs...>)) test(
+ std::nullptr_t);
+
+ //! non templated option
+ template <class T>
+ static decltype(testSignature(&T::getMaxStepLength)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! value traits type
+ template <class TProcess, typename TReturn, typename... TArgs>
+ bool constexpr has_method_getMaxStepLength_v =
+ has_method_getMaxStepLength<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika
diff --git a/corsika/detail/framework/process/DecayProcess.hpp b/corsika/detail/framework/process/DecayProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b5664418a4d6db1692be00294c82d7e224782dbc
--- /dev/null
+++ b/corsika/detail/framework/process/DecayProcess.hpp
@@ -0,0 +1,94 @@
+/*
+ * (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/process/ProcessTraits.hpp>
+
+namespace corsika {
+
+ /**
+ traits test for DecayProcess::doDecay method
+ */
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ struct has_method_doDecay : public detail::has_method_signature<TReturn, TArgs...> {
+
+ //! type of process to be studied
+ typedef std::decay_t<TProcess> process_type;
+
+ ///! method signature
+ using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! signature of templated method
+ template <class T>
+ static decltype(testSignature(&T::template doDecay<TArgs...>)) test(std::nullptr_t);
+
+ //! signature of non-templated method
+ template <class T>
+ static decltype(testSignature(&T::doDecay)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<process_type>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! @file DecayProcess.hpp
+ //! value traits type
+ template <class TProcess, typename TReturn, typename... TArgs>
+ bool constexpr has_method_doDecay_v =
+ has_method_doDecay<TProcess, TReturn, TArgs...>::value;
+
+ /**
+ traits test for DecayProcess::getLifetime method
+ */
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ struct has_method_getLifetime : public detail::has_method_signature<TReturn, TArgs...> {
+
+ //! the moethod signature
+ using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! signature of templated method
+ template <class T>
+ static decltype(testSignature(&T::template getLifetime<TArgs...>)) test(
+ std::nullptr_t);
+
+ //! signature of non-templated method
+ template <class T>
+ static decltype(testSignature(&T::getLifetime)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! value traits type
+ template <class TProcess, typename TReturn, typename... TArgs>
+ bool constexpr has_method_getLifetime_v =
+ has_method_getLifetime<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika
diff --git a/corsika/detail/framework/process/InteractionProcess.hpp b/corsika/detail/framework/process/InteractionProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8fb157a688c8c4828ea5c0b3eff410af9ad9aa76
--- /dev/null
+++ b/corsika/detail/framework/process/InteractionProcess.hpp
@@ -0,0 +1,95 @@
+/*
+ * (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/process/ProcessTraits.hpp>
+
+namespace corsika {
+
+ /**
+ traits test for InteractionProcess::doInteraction method
+ */
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ struct has_method_doInteract : public detail::has_method_signature<TReturn, TArgs...> {
+
+ ///! method signature
+ using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! signature of templated method
+ template <class T>
+ static decltype(testSignature(&T::template doInteraction<TArgs...>)) test(
+ std::nullptr_t);
+
+ //! signature of non-templated method
+ template <class T>
+ static decltype(testSignature(&T::doInteraction)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! @file BoundaryCrossingProcess.hpp
+ //! value traits type
+ template <class TProcess, typename TReturn, typename... TArgs>
+ bool constexpr has_method_doInteract_v =
+ has_method_doInteract<TProcess, TReturn, TArgs...>::value;
+
+ /**
+ traits test for InteractionProcess::getInteractionLength method
+ */
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ struct has_method_getInteractionLength
+ : public detail::has_method_signature<TReturn, TArgs...> {
+
+ ///! method signature
+ using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::template getInteractionLength<TArgs...>)) test(
+ std::nullptr_t);
+
+ //! non templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::getInteractionLength)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! @file BoundaryCrossingProcess.hpp
+ //! value traits type
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ bool constexpr has_method_getInteractionLength_v =
+ has_method_getInteractionLength<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index b93331d05079962f5e710ac8de952935eb289517..6c9bd789289a6fa24204bfe10df34e5413b5b6fa 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -26,6 +26,23 @@
namespace corsika {
+ template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
+ int IndexProcess2>
+ inline ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::ProcessSequence(TProcess1 in_A, TProcess2 in_B)
+ : A_(in_A)
+ , B_(in_B) {
+
+ // make sure only BaseProcess types TProcess1/2 are passed
+
+ static_assert(is_process_v<TProcess1>,
+ "can only use process derived from BaseProcess in "
+ "ProcessSequence, for Process 1");
+ static_assert(is_process_v<TProcess2>,
+ "can only use process derived from BaseProcess in "
+ "ProcessSequence, for Process 2");
+ }
+
template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
int IndexProcess2>
template <typename TParticle>
@@ -34,18 +51,41 @@ namespace corsika {
IndexProcess2>::doBoundaryCrossing(TParticle& particle,
typename TParticle::node_type const& from,
typename TParticle::node_type const& to) {
+
ProcessReturn ret = ProcessReturn::Ok;
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
- process1_type> ||
- t1ProcSeq) {
- ret |= A_.doBoundaryCrossing(particle, from, to);
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (is_boundary_process_v<process1_type> ||
+ process1_type::is_process_sequence) {
+
+ // interface checking on TProcess1
+ static_assert(
+ has_method_doBoundaryCrossing_v<TProcess1, ProcessReturn, TParticle>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+ "required for "
+ "BoundaryCrossingProcess<TDerived>. ");
+
+ ret |= A_.doBoundaryCrossing(particle, from, to);
+ }
}
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>,
- process2_type> ||
- t2ProcSeq) {
- ret |= B_.doBoundaryCrossing(particle, from, to);
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr (is_boundary_process_v<process2_type> ||
+ process2_type::is_process_sequence) {
+
+ // interface checking on TProcess2
+ static_assert(
+ has_method_doBoundaryCrossing_v<TProcess2, ProcessReturn, TParticle>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+ "required for "
+ "BoundaryCrossingProcess<TDerived>. ");
+
+ ret |= B_.doBoundaryCrossing(particle, from, to);
+ }
}
return ret;
@@ -59,18 +99,50 @@ namespace corsika {
doContinuous(TParticle& particle, TTrack& vT,
[[maybe_unused]] ContinuousProcessIndex const limitId) {
ProcessReturn ret = ProcessReturn::Ok;
- if constexpr (t1ProcSeq) {
- ret |= A_.doContinuous(particle, vT, limitId);
- } else if constexpr (is_continuous_process_v<process1_type>) {
- ret |=
- A_.doContinuous(particle, vT, limitId == ContinuousProcessIndex(IndexProcess1));
+
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (process1_type::is_process_sequence) {
+
+ ret |= A_.doContinuous(particle, vT, limitId);
+ } else if constexpr (is_continuous_process_v<process1_type>) {
+
+ // interface checking on TProcess1
+ static_assert(
+ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle&, TTrack&> ||
+ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle&,
+ TTrack const&> ||
+ has_method_doContinuous_v<TProcess1, ProcessReturn, TParticle const&,
+ TTrack const&>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doContinuous(TParticle [const]&,TTrack [const]&,bool)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ ret |= A_.doContinuous(particle, vT,
+ limitId == ContinuousProcessIndex(IndexProcess1));
+ }
}
- if constexpr (t2ProcSeq) {
- ret |= B_.doContinuous(particle, vT, limitId);
- } else if constexpr (is_continuous_process_v<process2_type>) {
- ret |=
- B_.doContinuous(particle, vT, limitId == ContinuousProcessIndex(IndexProcess2));
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr (process2_type::is_process_sequence) {
+ ret |= B_.doContinuous(particle, vT, limitId);
+ } else if constexpr (is_continuous_process_v<process2_type>) {
+
+ // interface checking on TProcess2
+ static_assert(
+ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle&, TTrack&> ||
+ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle&,
+ TTrack const&> ||
+ has_method_doContinuous_v<TProcess2, ProcessReturn, TParticle const&,
+ TTrack const&>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doContinuous(TParticle [const]&,TTrack [const]&,bool)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ ret |= B_.doContinuous(particle, vT,
+ limitId == ContinuousProcessIndex(IndexProcess2));
+ }
}
return ret;
@@ -79,15 +151,42 @@ namespace corsika {
template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
int IndexProcess2>
template <typename TSecondaries>
- inline void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
- IndexProcess2>::doSecondaries(TSecondaries& vS) {
- if constexpr (std::is_base_of_v<SecondariesProcess<process1_type>, process1_type> ||
- t1ProcSeq) {
- A_.doSecondaries(vS);
+ inline void
+ ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
+ IndexProcess2>::doSecondaries([[maybe_unused]] TSecondaries& vS) {
+
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (is_secondaries_process_v<process1_type> ||
+ process1_type::is_process_sequence) {
+
+ // interface checking on TProcess1
+ static_assert(
+ has_method_doSecondaries_v<TProcess1, void, TSecondaries&> ||
+ has_method_doSecondaries_v<TProcess1, void, TSecondaries const&>,
+ "TDerived has no method with correct signature \"void "
+ "doSecondaries(TStackView [const]&)\" required for "
+ "SecondariesProcessProcess<TDerived>. ");
+
+ A_.doSecondaries(vS);
+ }
}
- if constexpr (std::is_base_of_v<SecondariesProcess<process2_type>, process2_type> ||
- t2ProcSeq) {
- B_.doSecondaries(vS);
+
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr (is_secondaries_process_v<process2_type> ||
+ process2_type::is_process_sequence) {
+
+ // interface checking on TProcess2
+ static_assert(
+ has_method_doSecondaries_v<TProcess2, void, TSecondaries&> ||
+ has_method_doSecondaries_v<TProcess2, void, TSecondaries const&>,
+ "TDerived has no method with correct signature \"void "
+ "doSecondaries(TStackView [const]&)\" required for "
+ "SecondariesProcessProcess<TDerived>. ");
+
+ B_.doSecondaries(vS);
+ }
}
}
@@ -96,13 +195,21 @@ namespace corsika {
inline bool ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
IndexProcess2>::checkStep() {
bool ret = false;
- if constexpr (std::is_base_of_v<StackProcess<process1_type>, process1_type> ||
- (t1ProcSeq && !t1SwitchProcSeq)) {
- ret |= A_.checkStep();
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (is_stack_process_v<process1_type> ||
+ (process1_type::is_process_sequence &&
+ !process1_type::is_switch_process_sequence)) {
+ ret |= A_.checkStep();
+ }
}
- if constexpr (std::is_base_of_v<StackProcess<process2_type>, process2_type> ||
- (t2ProcSeq && !t2SwitchProcSeq)) {
- ret |= B_.checkStep();
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr (is_stack_process_v<process2_type> ||
+ (process2_type::is_process_sequence &&
+ !process2_type::is_switch_process_sequence)) {
+ ret |= B_.checkStep();
+ }
}
return ret;
}
@@ -112,13 +219,37 @@ namespace corsika {
template <typename TStack>
inline void ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
IndexProcess2>::doStack(TStack& stack) {
- if constexpr (std::is_base_of_v<StackProcess<process1_type>, process1_type> ||
- (t1ProcSeq && !t1SwitchProcSeq)) {
- if (A_.checkStep()) { A_.doStack(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> ||
+ (process1_type::is_process_sequence &&
+ !process1_type::is_switch_process_sequence)) {
+
+ // interface checking on TProcess1
+ static_assert(has_method_doStack_v<TProcess1, void, TStack&> ||
+ has_method_doStack_v<TProcess1, void, TStack const&>,
+ "TDerived has no method with correct signature \"void "
+ "doStack(TStack [const]&)\" required for "
+ "StackProcess<TDerived>. ");
+
+ if (A_.checkStep()) { A_.doStack(stack); }
+ }
}
- if constexpr (std::is_base_of_v<StackProcess<process2_type>, process2_type> ||
- (t2ProcSeq && !t2SwitchProcSeq)) {
- if (B_.checkStep()) { B_.doStack(stack); }
+ 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> ||
+ (process2_type::is_process_sequence &&
+ !process2_type::is_switch_process_sequence)) {
+
+ // interface checking on TProcess1
+ static_assert(has_method_doStack_v<TProcess2, void, TStack&> ||
+ has_method_doStack_v<TProcess2, void, TStack const&>,
+ "TDerived has no method with correct signature \"void "
+ "doStack(TStack [const]&)\" required for "
+ "StackProcess<TDerived>. ");
+
+ if (B_.checkStep()) { B_.doStack(stack); }
+ }
}
}
@@ -132,22 +263,44 @@ namespace corsika {
ContinuousProcessStepLength max_length(std::numeric_limits<double>::infinity() *
meter);
- if constexpr (t1ProcSeq) {
- ContinuousProcessStepLength const step = A_.getMaxStepLength(particle, vTrack);
- max_length = std::min(max_length, step);
- } else if constexpr (is_continuous_process_v<process1_type>) {
- ContinuousProcessStepLength const step(A_.getMaxStepLength(particle, vTrack),
- ContinuousProcessIndex(IndexProcess1));
- max_length = std::min(max_length, step);
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (process1_type::is_process_sequence) {
+ ContinuousProcessStepLength const step = A_.getMaxStepLength(particle, vTrack);
+ max_length = std::min(max_length, step);
+ } else if constexpr (is_continuous_process_v<process1_type>) {
+
+ // interface checking on TProcess1
+ static_assert(has_method_getMaxStepLength_v<TProcess1, LengthType,
+ TParticle const&, TTrack const&>,
+ "TDerived has no method with correct signature \"LengthType "
+ "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ ContinuousProcessStepLength const step(A_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess1));
+ max_length = std::min(max_length, step);
+ }
}
- // return ContinuousProcessStepLength(std::min(max_length, len), IndexProcess1);
- if constexpr (t2ProcSeq) {
- ContinuousProcessStepLength const step = B_.getMaxStepLength(particle, vTrack);
- max_length = std::min(max_length, step);
- } else if constexpr (is_continuous_process_v<process2_type>) {
- ContinuousProcessStepLength const step(B_.getMaxStepLength(particle, vTrack),
- ContinuousProcessIndex(IndexProcess2));
- max_length = std::min(max_length, step);
+
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr (process2_type::is_process_sequence) {
+ ContinuousProcessStepLength const step = B_.getMaxStepLength(particle, vTrack);
+ max_length = std::min(max_length, step);
+ } else if constexpr (is_continuous_process_v<process2_type>) {
+
+ // interface checking on TProcess2
+ static_assert(has_method_getMaxStepLength_v<TProcess2, LengthType,
+ TParticle const&, TTrack const&>,
+ "TDerived has no method with correct signature \"LengthType "
+ "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ ContinuousProcessStepLength const step(B_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess2));
+ max_length = std::min(max_length, step);
+ }
}
return max_length;
}
@@ -161,13 +314,19 @@ namespace corsika {
InverseGrammageType tot = 0 * meter * meter / gram; // default value
- if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> ||
- t1ProcSeq) {
- tot += A_.getInverseInteractionLength(particle);
+ 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> ||
+ process1_type::is_process_sequence) {
+ tot += A_.getInverseInteractionLength(particle);
+ }
}
- if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> ||
- t2ProcSeq) {
- tot += B_.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> ||
+ process2_type::is_process_sequence) {
+ tot += B_.getInverseInteractionLength(particle);
+ }
}
return tot;
}
@@ -181,39 +340,62 @@ namespace corsika {
[[maybe_unused]] InverseGrammageType lambda_inv_select,
[[maybe_unused]] InverseGrammageType lambda_inv_sum) {
- // TODO: add check for lambda_inv_select>lambda_inv_tot
-
- if constexpr (t1ProcSeq) {
- // if A is a process sequence --> check inside
- ProcessReturn const ret =
- 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>) {
- // 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
- if (lambda_inv_select <= lambda_inv_sum) {
- A_.doInteraction(view);
- return ProcessReturn::Interacted;
- }
- } // end branch A
-
- if constexpr (t2ProcSeq) {
- // 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>) {
- // if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
- // soon as SecondaryView::parent() is migrated!
- // check if we should execute THIS process and then EXIT
- if (lambda_inv_select <= lambda_inv_sum) {
- B_.doInteraction(view);
- return ProcessReturn::Interacted;
- }
- } // end branch B_
+ // TODO: add check for lambda_inv_select > lambda_inv_tot
+
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (process1_type::is_process_sequence) {
+ // if A is a process sequence --> check inside
+ ProcessReturn const ret =
+ 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>) {
+ // 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
+ if (lambda_inv_select <= lambda_inv_sum) {
+
+ // interface checking on TProcess1
+ static_assert(has_method_doInteract_v<TProcess1, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doInteraction(TSecondaryView&)\" required for "
+ "InteractionProcess<TDerived>. ");
+
+ A_.template doInteraction(view);
+ return ProcessReturn::Interacted;
+ }
+ } // end branch A
+ }
+
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+
+ 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>) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
+ // soon as SecondaryView::parent() is migrated!
+ // check if we should execute THIS process and then EXIT
+ if (lambda_inv_select <= lambda_inv_sum) {
+
+ // interface checking on TProcess1
+ static_assert(has_method_doInteract_v<TProcess2, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doInteraction(TSecondaryView&)\" required for "
+ "InteractionProcess<TDerived>. ");
+
+ B_.doInteraction(view);
+ return ProcessReturn::Interacted;
+ }
+ } // end branch B_
+ }
return ProcessReturn::Ok;
}
@@ -226,13 +408,19 @@ namespace corsika {
InverseTimeType tot = 0 / second; // default value
- if constexpr (std::is_base_of_v<DecayProcess<process1_type>, process1_type> ||
- t1ProcSeq) {
- tot += A_.getInverseLifetime(particle);
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (is_decay_process_v<process1_type> ||
+ process1_type::is_process_sequence) {
+ tot += A_.getInverseLifetime(particle);
+ }
}
- if constexpr (std::is_base_of_v<DecayProcess<process2_type>, process2_type> ||
- t2ProcSeq) {
- tot += B_.getInverseLifetime(particle);
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ if constexpr (is_decay_process_v<process2_type> ||
+ process2_type::is_process_sequence) {
+ tot += B_.getInverseLifetime(particle);
+ }
}
return tot;
}
@@ -249,34 +437,54 @@ namespace corsika {
// TODO: add check for decay_inv_select>decay_inv_tot
- if constexpr (t1ProcSeq) {
- // if A_ is a process sequence --> check inside
- 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>) {
- // 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
- if (decay_inv_select <= decay_inv_sum) { // more pedagogical: rndm_select <
- // decay_inv_sum / decay_inv_tot
- A_.doDecay(view);
- return ProcessReturn::Decayed;
- }
- } // end branch A_
-
- if constexpr (t2ProcSeq) {
- // 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>) {
- // 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
- if (decay_inv_select <= decay_inv_sum) {
- B_.doDecay(view);
- return ProcessReturn::Decayed;
- }
- } // end branch B_
+ if constexpr (is_process_v<process1_type>) { // to protect from further compiler
+ // errors if process1_type is invalid
+ if constexpr (process1_type::is_process_sequence) {
+ // if A_ is a process sequence --> check inside
+ 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 (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
+ if (decay_inv_select <= decay_inv_sum) { // more pedagogical: rndm_select <
+ // decay_inv_sum / decay_inv_tot
+ // interface checking on TProcess1
+ static_assert(has_method_doDecay_v<TProcess1, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doDecay(TSecondaryView&)\" required for "
+ "DecayProcess<TDerived>. ");
+
+ A_.doDecay(view);
+ return ProcessReturn::Decayed;
+ }
+ } // end branch A_
+ }
+
+ if constexpr (is_process_v<process2_type>) { // to protect from further compiler
+ // errors if process2_type is invalid
+ 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 (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
+ if (decay_inv_select <= decay_inv_sum) {
+
+ // interface checking on TProcess1
+ static_assert(has_method_doDecay_v<TProcess2, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doDecay(TSecondaryView&)\" required for "
+ "DecayProcess<TDerived>. ");
+
+ B_.doDecay(view);
+ return ProcessReturn::Decayed;
+ }
+ } // end branch B_
+ }
+
return ProcessReturn::Ok;
}
@@ -285,15 +493,12 @@ namespace corsika {
**/
namespace detail {
// need helper alias to achieve this:
- template <typename TProcess1, typename TProcess2,
- typename = typename std::enable_if_t<
- contains_stack_process_v<TProcess1> ||
- std::is_base_of_v<StackProcess<typename std::decay_t<TProcess1>>,
- typename std::decay_t<TProcess1>> ||
- contains_stack_process_v<TProcess2> ||
- std::is_base_of_v<StackProcess<typename std::decay_t<TProcess2>>,
- typename std::decay_t<TProcess2>>,
- int>>
+ template <
+ typename TProcess1, typename TProcess2,
+ typename = typename std::enable_if_t<
+ contains_stack_process_v<TProcess1> || is_stack_process_v<TProcess1> ||
+ contains_stack_process_v<TProcess2> || is_stack_process_v<TProcess2>,
+ int>>
using enable_if_stack = ProcessSequence<TProcess1, TProcess2>;
} // namespace detail
diff --git a/corsika/detail/framework/process/SecondariesProcess.hpp b/corsika/detail/framework/process/SecondariesProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0310bf3dfa2b73826c10c2ad25bf4eb2d4ab5a21
--- /dev/null
+++ b/corsika/detail/framework/process/SecondariesProcess.hpp
@@ -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
+
+#include <corsika/framework/process/ProcessTraits.hpp>
+
+namespace corsika {
+
+ /**
+ traits test for SecondariesProcess::doSecondaries method
+ */
+ template <class TProcess, typename TReturn, typename... TArg>
+ struct has_method_doSecondaries
+ : public detail::has_method_signature<TReturn, TArg...> {
+
+ //! method signature
+ using detail::has_method_signature<TReturn, TArg...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::template doSecondaries<TArg...>)) test(
+ std::nullptr_t);
+
+ //! non templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::doSecondaries)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! @file DecayProcess.hpp
+ //! value traits type
+ template <class TProcess, typename TReturn, typename... TArg>
+ bool constexpr has_method_doSecondaries_v =
+ has_method_doSecondaries<TProcess, TReturn, TArg...>::value;
+
+} // namespace corsika
diff --git a/corsika/detail/framework/process/StackProcess.hpp b/corsika/detail/framework/process/StackProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..68f350adf1bd1db6cc2be74b329e6a4824b2df07
--- /dev/null
+++ b/corsika/detail/framework/process/StackProcess.hpp
@@ -0,0 +1,57 @@
+/*
+ * (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/process/ProcessTraits.hpp>
+
+namespace corsika {
+
+ /**
+ traits test for StackProcess::doStack method
+ */
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ struct has_method_doStack : public detail::has_method_signature<TReturn, TArgs...> {
+
+ //! method signature
+ using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::template doStack<TArgs...>)) test(std::nullptr_t);
+
+ //! templated-template parameter option
+ template <template <typename> typename T>
+ static decltype(testSignature(&T<TArgs...>::template doStack)) test(std::nullptr_t);
+
+ //! non-templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::doStack)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! @file StackProcess.hpp
+ //! value traits type
+ template <class TProcess, typename TReturn, typename... TArgs>
+ bool constexpr has_method_doStack_v =
+ has_method_doStack<TProcess, TReturn, TArgs...>::value;
+
+} // namespace corsika
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index c4f8ede12f859e7c063f91cabbcc01af8b99cbc4..2f87ade6d329ab249f5a61f21b79c72d07d315fd 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -27,113 +27,183 @@
namespace corsika {
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
- template <typename TParticle, typename TVTNType>
- inline ProcessReturn
- SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
- IndexProcess2>::doBoundaryCrossing(TParticle& particle,
- TVTNType const& from,
- TVTNType const& to) {
- switch (select_(particle)) {
- case SwitchResult::First: {
+ template <typename TParticle>
+ inline ProcessReturn SwitchProcessSequence<
+ TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+ IndexProcess2>::doBoundaryCrossing(TParticle& particle,
+ 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> ||
+ process1_type::is_process_sequence) {
+
+ // interface checking on TSequence
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
- process1_type> ||
- t1ProcSeq) {
- return A_.doBoundaryCrossing(particle, from, to);
+ process1_type>) {
+
+ static_assert(
+ has_method_doBoundaryCrossing_v<TSequence, ProcessReturn, TParticle&>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+ "required for "
+ "BoundaryCrossingProcess<TDerived>. ");
}
- break;
+
+ return A_.doBoundaryCrossing(particle, from, to);
}
- case SwitchResult::Second: {
+ } else {
+
+ if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>,
+ process2_type> ||
+ process2_type::is_process_sequence) {
+
+ // interface checking on USequence
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>,
- process2_type> ||
- t2ProcSeq) {
- return B_.doBoundaryCrossing(particle, from, to);
+ process2_type>) {
+
+ static_assert(
+ has_method_doBoundaryCrossing_v<USequence, ProcessReturn, TParticle>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doBoundaryCrossing(TParticle&, VolumeNode const&, VolumeNode const&)\" "
+ "required for "
+ "BoundaryCrossingProcess<TDerived>. ");
}
- break;
+
+ return B_.doBoundaryCrossing(particle, from, to);
}
}
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TTrack>
inline ProcessReturn SwitchProcessSequence<
- TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ TCondition, TSequence, USequence, IndexStart, IndexProcess1,
IndexProcess2>::doContinuous(TParticle& particle, TTrack& vT,
ContinuousProcessIndex const idLimit) {
- switch (select_(particle)) {
- case SwitchResult::First: {
- if constexpr (t1ProcSeq) { return A_.doContinuous(particle, vT, idLimit); }
- if constexpr (is_continuous_process_v<process1_type>) {
- return A_.doContinuous(particle, vT,
- idLimit == ContinuousProcessIndex(IndexProcess1));
- }
- break;
+ if (select_(particle)) {
+ if constexpr (process1_type::is_process_sequence) {
+ return A_.doContinuous(particle, vT, idLimit);
}
- case SwitchResult::Second: {
- if constexpr (t2ProcSeq) { return B_.doContinuous(particle, vT, idLimit); }
- if constexpr (is_continuous_process_v<process2_type>) {
- return B_.doContinuous(particle, vT,
- idLimit == ContinuousProcessIndex(IndexProcess2));
- }
- break;
+ if constexpr (is_continuous_process_v<process1_type>) {
+
+ static_assert(
+ has_method_doContinuous_v<TSequence, ProcessReturn, TParticle&, TTrack&> ||
+ has_method_doContinuous_v<TSequence, ProcessReturn, TParticle&,
+ TTrack const&> ||
+ has_method_doContinuous_v<TSequence, ProcessReturn, TParticle const&,
+ TTrack const&>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doContinuous(TParticle[const]&,TTrack[const]&,bool)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ return A_.doContinuous(particle, vT,
+ idLimit == ContinuousProcessIndex(IndexProcess1));
+ }
+ } else {
+ if constexpr (process2_type::is_process_sequence) {
+ return B_.doContinuous(particle, vT, idLimit);
+ }
+ if constexpr (is_continuous_process_v<process2_type>) {
+
+ // interface checking on USequence
+ static_assert(
+ has_method_doContinuous_v<USequence, ProcessReturn, TParticle&, TTrack&> ||
+ has_method_doContinuous_v<USequence, ProcessReturn, TParticle&,
+ TTrack const&> ||
+ has_method_doContinuous_v<USequence, ProcessReturn, TParticle const&,
+ TTrack const&>,
+ "TDerived has no method with correct signature \"ProcessReturn "
+ "doContinuous(TParticle [const]&,TTrack[const]&,bool)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ return B_.doContinuous(particle, vT,
+ idLimit == ContinuousProcessIndex(IndexProcess2));
}
}
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TSecondaries>
inline void
- SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
IndexProcess2>::doSecondaries(TSecondaries& vS) {
const auto& particle = vS.parent();
- switch (select_(particle)) {
- case SwitchResult::First: {
- if constexpr (std::is_base_of_v<SecondariesProcess<process1_type>,
- process1_type> ||
- t1ProcSeq) {
- A_.doSecondaries(vS);
- }
- break;
+ if (select_(particle)) {
+ if constexpr (std::is_base_of_v<SecondariesProcess<process1_type>, process1_type> ||
+ process1_type::is_process_sequence) {
+
+ // interface checking on TSequence
+ static_assert(
+ has_method_doSecondaries_v<TSequence, void, TSecondaries&> ||
+ has_method_doSecondaries_v<TSequence, void, TSecondaries const&>,
+ "TDerived has no method with correct signature \"void "
+ "doSecondaries(TStackView [const]&)\" required for "
+ "SecondariesProcessProcess<TDerived>. ");
+
+ A_.doSecondaries(vS);
}
- case SwitchResult::Second: {
- if constexpr (std::is_base_of_v<SecondariesProcess<process2_type>,
- process2_type> ||
- t2ProcSeq) {
- B_.doSecondaries(vS);
- }
- break;
+ } else {
+ if constexpr (std::is_base_of_v<SecondariesProcess<process2_type>, process2_type> ||
+ process2_type::is_process_sequence) {
+
+ // interface checking on USequence
+ static_assert(
+ has_method_doSecondaries_v<USequence, void, TSecondaries&> ||
+ has_method_doSecondaries_v<USequence, void, TSecondaries const&>,
+ "TDerived has no method with correct signature \"void "
+ "doSecondaries(TStackView [const]&)\" required for "
+ "SecondariesProcessProcess<TDerived>. ");
+
+ B_.doSecondaries(vS);
}
}
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TParticle, typename TTrack>
inline ContinuousProcessStepLength
- SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
IndexProcess2>::getMaxStepLength(TParticle& particle,
TTrack& vTrack) {
- switch (select_(particle)) {
- case SwitchResult::First: {
- if constexpr (t1ProcSeq) { return A_.getMaxStepLength(particle, vTrack); }
- if constexpr (is_continuous_process_v<process1_type>) {
- return ContinuousProcessStepLength(A_.getMaxStepLength(particle, vTrack),
- ContinuousProcessIndex(IndexProcess1));
- }
- break;
+ if (select_(particle)) {
+ if constexpr (process1_type::is_process_sequence) {
+ return A_.getMaxStepLength(particle, vTrack);
}
- case SwitchResult::Second: {
- if constexpr (t2ProcSeq) { return B_.getMaxStepLength(particle, vTrack); }
- if constexpr (is_continuous_process_v<process2_type>) {
- return ContinuousProcessStepLength(B_.getMaxStepLength(particle, vTrack),
- ContinuousProcessIndex(IndexProcess2));
- }
- break;
+ if constexpr (is_continuous_process_v<process1_type>) {
+
+ // interface checking on TSequence
+ static_assert(has_method_getMaxStepLength_v<TSequence, LengthType,
+ TParticle const&, TTrack const&>,
+ "TDerived has no method with correct signature \"LengthType "
+ "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ return ContinuousProcessStepLength(A_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess1));
+ }
+ } else {
+ if constexpr (process2_type::is_process_sequence) {
+ return B_.getMaxStepLength(particle, vTrack);
+ }
+ if constexpr (is_continuous_process_v<process2_type>) {
+
+ // interface checking on USequence
+ static_assert(has_method_getMaxStepLength_v<USequence, LengthType,
+ TParticle const&, TTrack const&>,
+ "TDerived has no method with correct signature \"LengthType "
+ "getMaxStepLength(TParticle const&, TTrack const&)\" required for "
+ "ContinuousProcess<TDerived>. ");
+
+ return ContinuousProcessStepLength(B_.getMaxStepLength(particle, vTrack),
+ ContinuousProcessIndex(IndexProcess2));
}
}
@@ -141,182 +211,169 @@ namespace corsika {
return ContinuousProcessStepLength(std::numeric_limits<double>::infinity() * meter);
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TParticle>
inline InverseGrammageType SwitchProcessSequence<
- TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ TCondition, TSequence, USequence, IndexStart, IndexProcess1,
IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
- switch (select_(particle)) {
- case SwitchResult::First: {
- if constexpr (std::is_base_of_v<InteractionProcess<process1_type>,
- process1_type> ||
- t1ProcSeq) {
- return A_.getInverseInteractionLength(particle);
- }
- break;
+ if (select_(particle)) {
+ if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> ||
+ process1_type::is_process_sequence) {
+ return A_.getInverseInteractionLength(particle);
}
- case SwitchResult::Second: {
- if constexpr (std::is_base_of_v<InteractionProcess<process2_type>,
- process2_type> ||
- t2ProcSeq) {
- return B_.getInverseInteractionLength(particle);
- }
- break;
+
+ } else {
+
+ if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> ||
+ process2_type::is_process_sequence) {
+ return B_.getInverseInteractionLength(particle);
}
}
return 0 * meter * meter / gram; // default value
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TSecondaryView>
- inline ProcessReturn SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart,
+ inline ProcessReturn SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
IndexProcess1, IndexProcess2>::
selectInteraction(TSecondaryView& view,
[[maybe_unused]] InverseGrammageType lambda_inv_select,
[[maybe_unused]] InverseGrammageType lambda_inv_sum) {
- switch (select_(view.parent())) {
- case SwitchResult::First: {
- if constexpr (t1ProcSeq) {
- // if A_ is a process sequence --> check inside
- ProcessReturn const ret =
- 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>) {
- // 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
- if (lambda_inv_select < lambda_inv_sum) {
- A_.doInteraction(view);
- return ProcessReturn::Interacted;
- }
- } // end branch A_
- break;
- }
+ if (select_(view.parent())) {
+ if constexpr (process1_type::is_process_sequence) {
+ // if A_ is a process sequence --> check inside
+ ProcessReturn const ret =
+ 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>) {
+ // 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
+ if (lambda_inv_select < lambda_inv_sum) {
- case SwitchResult::Second: {
-
- if constexpr (t2ProcSeq) {
- // 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>) {
- // 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
- if (lambda_inv_select < lambda_inv_sum) {
- B_.doInteraction(view);
- return ProcessReturn::Interacted;
- }
- } // end branch B_
- break;
- }
+ // interface checking on TSequence
+ static_assert(has_method_doInteract_v<TSequence, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doInteraction(TSecondaryView&)\" required for "
+ "InteractionProcess<TDerived>. ");
+
+ A_.doInteraction(view);
+ return ProcessReturn::Interacted;
+ }
+ } // end branch A_
+
+ } else {
+
+ 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>) {
+ // 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
+ if (lambda_inv_select < lambda_inv_sum) {
+
+ // interface checking on TSequence
+ static_assert(has_method_doInteract_v<USequence, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doInteraction(TSecondaryView&)\" required for "
+ "InteractionProcess<TDerived>. ");
+
+ B_.doInteraction(view);
+ return ProcessReturn::Interacted;
+ }
+ } // end branch B_
}
return ProcessReturn::Ok;
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TParticle>
inline InverseTimeType
- SwitchProcessSequence<TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart, IndexProcess1,
IndexProcess2>::getInverseLifetime(TParticle&& particle) {
- switch (select_(particle)) {
- case SwitchResult::First: {
- if constexpr (std::is_base_of_v<DecayProcess<process1_type>, process1_type> ||
- t1ProcSeq) {
- return A_.getInverseLifetime(particle);
- }
- break;
+ if (select_(particle)) {
+ if constexpr (std::is_base_of_v<DecayProcess<process1_type>, process1_type> ||
+ process1_type::is_process_sequence) {
+ return A_.getInverseLifetime(particle);
}
- case SwitchResult::Second: {
- if constexpr (std::is_base_of_v<DecayProcess<process2_type>, process2_type> ||
- t2ProcSeq) {
- return B_.getInverseLifetime(particle);
- }
- break;
+ } else {
+
+ if constexpr (std::is_base_of_v<DecayProcess<process2_type>, process2_type> ||
+ process2_type::is_process_sequence) {
+ return B_.getInverseLifetime(particle);
}
}
return 0 / second; // default value
}
- template <typename TProcess1, typename TProcess2, typename TSelect, int IndexStart,
+ template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
// select decay process
template <typename TSecondaryView>
inline ProcessReturn SwitchProcessSequence<
- TProcess1, TProcess2, TSelect, IndexStart, IndexProcess1,
+ TCondition, TSequence, USequence, IndexStart, IndexProcess1,
IndexProcess2>::selectDecay(TSecondaryView& view,
[[maybe_unused]] InverseTimeType decay_inv_select,
[[maybe_unused]] InverseTimeType decay_inv_sum) {
- switch (select_(view.parent())) {
- case SwitchResult::First: {
- if constexpr (t1ProcSeq) {
- // if A_ is a process sequence --> check inside
- 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>) {
- // 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
- if (decay_inv_select < decay_inv_sum) {
- // more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot
- A_.doDecay(view);
- return ProcessReturn::Decayed;
- }
- } // end branch A_
- break;
- }
+ if (select_(view.parent())) {
+ if constexpr (process1_type::is_process_sequence) {
+ // if A_ is a process sequence --> check inside
+ 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>) {
+ // 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
+ if (decay_inv_select < decay_inv_sum) {
+ // more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot
- case SwitchResult::Second: {
-
- if constexpr (t2ProcSeq) {
- // 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>) {
- // 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
- if (decay_inv_select < decay_inv_sum) {
- B_.doDecay(view);
- return ProcessReturn::Decayed;
- }
- } // end branch B_
- break;
- }
+ // interface checking on TSequence
+ static_assert(has_method_doDecay_v<TSequence, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doDecay(TSecondaryView&)\" required for "
+ "DecayProcess<TDerived>. ");
+
+ A_.doDecay(view);
+ return ProcessReturn::Decayed;
+ }
+ } // end branch A_
+
+ } else {
+
+ 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>) {
+ // 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
+ if (decay_inv_select < decay_inv_sum) {
+
+ // interface checking on TSequence
+ static_assert(has_method_doDecay_v<USequence, void, TSecondaryView&>,
+ "TDerived has no method with correct signature \"void "
+ "doDecay(TSecondaryView&)\" required for "
+ "DecayProcess<TDerived>. ");
+
+ B_.doDecay(view);
+ return ProcessReturn::Decayed;
+ }
+ } // end branch B_
}
return ProcessReturn::Ok;
}
- /*
- /// traits marker to identify objectas ProcessSequence
- template <typename TProcess1, typename TProcess2, typename TSelect>
- struct is_process_sequence<ProcessSequence<typename std::decay_t<TProcess1>,
- typename std::decay_t<TProcess2>,
- typename std::decay_t<TSelect>>>
- : std::true_type {};
- */
- /// traits marker to identify objectas ProcessSequence
- template <typename TProcess1, typename TProcess2, typename TSelect>
- struct is_process_sequence<SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
- : std::true_type {};
-
- template <typename TProcess1, typename TProcess2, typename TSelect>
- struct is_process_sequence<SwitchProcessSequence<TProcess1, TProcess2, TSelect>&>
- : std::true_type {};
-
- /// traits marker to identify objectas SwitchProcessSequence
- template <typename TProcess1, typename TProcess2, typename TSelect>
- struct is_switch_process_sequence<SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
- : std::true_type {};
-
} // namespace corsika
diff --git a/corsika/detail/framework/random/RNGManager.inl b/corsika/detail/framework/random/RNGManager.inl
index 575ab43b73e48efe733e2bde64a99503bce44199..1d0da48313f622e3c9b567127ed2104b15fc3029 100644
--- a/corsika/detail/framework/random/RNGManager.inl
+++ b/corsika/detail/framework/random/RNGManager.inl
@@ -8,6 +8,10 @@
#pragma once
+#include <iterator>
+#include <random>
+#include <sstream>
+
namespace corsika {
inline void RNGManager::registerRandomStream(string_type const& pStreamName) {
@@ -48,13 +52,16 @@ namespace corsika {
inline void RNGManager::seedAll(void) {
std::random_device rd;
- std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()};
- for (auto& entry : rngs_) {
- std::vector<std::uint32_t> seeds(1);
- sseq.generate(seeds.begin(), seeds.end());
- std::uint32_t seed = seeds[0];
- CORSIKA_LOG_TRACE("Random seed stream {} seed {}", entry.first, seed);
- entry.second.seed(seed);
+
+ for (auto& [streamName, rng] : rngs_) {
+ std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()}; // 6 really random values
+
+ // for logging collect sseq input values in string
+ std::stringstream ss;
+ sseq.param(std::ostream_iterator<int>{ss, " "});
+ CORSIKA_LOG_DEBUG("Random seed stream {} seed {}", streamName, ss.str());
+
+ rng.seed(sseq);
}
}
diff --git a/corsika/detail/framework/utility/CorsikaData.inl b/corsika/detail/framework/utility/CorsikaData.inl
index e2ea1e99bc3742783020312f725419fbd0d63846..fc638b81a7f2e1a7c66908835d06182d99e69610 100644
--- a/corsika/detail/framework/utility/CorsikaData.inl
+++ b/corsika/detail/framework/utility/CorsikaData.inl
@@ -7,14 +7,15 @@
*/
#pragma once
+#include <boost/filesystem/path.hpp>
+
#include <cstdlib>
#include <stdexcept>
#include <string>
-inline std::string corsika::corsika_data(std::string const& key) {
+inline boost::filesystem::path corsika::corsika_data(boost::filesystem::path const& key) {
if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
- auto const path = std::string(p) + "/" + key;
- return path;
+ return boost::filesystem::path(p) / key;
} else {
throw std::runtime_error("CORSIKA_DATA not set");
}
diff --git a/corsika/detail/media/NuclearComposition.inl b/corsika/detail/media/NuclearComposition.inl
index 85582510e5cb0e115c7e97133ee215de1193b2b2..4995799d58ba35e30c09448946d67c50f8e29ade 100644
--- a/corsika/detail/media/NuclearComposition.inl
+++ b/corsika/detail/media/NuclearComposition.inl
@@ -101,6 +101,10 @@ namespace corsika {
// must be updated, too!
inline size_t NuclearComposition::getHash() const { return hash_; }
+ inline bool NuclearComposition::operator==(NuclearComposition const& v) const {
+ return v.hash_ == hash_;
+ }
+
inline void NuclearComposition::updateHash() {
std::vector<std::size_t> hashes;
for (float ifrac : this->getFractions()) hashes.push_back(std::hash<float>{}(ifrac));
diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl
index b5b975818f585dc19a70afbf106360b33e5b6308..f77e763ac8a493706e44673250f2446a7d69a0e4 100644
--- a/corsika/detail/modules/TrackWriter.inl
+++ b/corsika/detail/modules/TrackWriter.inl
@@ -20,8 +20,8 @@ namespace corsika {
template <typename TOutput>
template <typename TParticle, typename TTrack>
- ProcessReturn TrackWriter<TOutput>::doContinuous(const TParticle& vP, const TTrack& vT,
- bool const) {
+ inline ProcessReturn TrackWriter<TOutput>::doContinuous(TParticle const& vP,
+ TTrack const& vT, bool const) {
auto const start = vT.getPosition(0).getCoordinates();
auto const end = vT.getPosition(1).getCoordinates();
@@ -34,8 +34,8 @@ namespace corsika {
template <typename TOutput>
template <typename TParticle, typename TTrack>
- inline LengthType TrackWriter<TOutput>::getMaxStepLength(const TParticle&,
- const TTrack&) {
+ inline LengthType TrackWriter<TOutput>::getMaxStepLength(TParticle const&,
+ TTrack const&) {
return meter * std::numeric_limits<double>::infinity();
}
diff --git a/corsika/detail/modules/conex/CONEXhybrid.inl b/corsika/detail/modules/conex/CONEXhybrid.inl
index 5b350013f27d75eeaf7ed2ed8f1ea5c0bed6cb61..6ab443e49bb7920d1e6a544ec84941236c8a1192 100644
--- a/corsika/detail/modules/conex/CONEXhybrid.inl
+++ b/corsika/detail/modules/conex/CONEXhybrid.inl
@@ -54,7 +54,8 @@ namespace corsika {
return b.normalized();
})}
- , y_sf_{showerAxis_.getDirection().cross(x_sf_)} {
+ , y_sf_{showerAxis_.getDirection().cross(x_sf_)}
+ , energy_em_(0_GeV) {
CORSIKA_LOG_DEBUG("x_sf (conexObservationCS): {}",
x_sf_.getComponents(conexObservationCS_));
@@ -180,6 +181,7 @@ namespace corsika {
double const E = energy / 1_GeV;
double const m = mass / 1_GeV;
+ energy_em_ += energy;
CORSIKA_LOG_DEBUG("CONEXhybrid: removing {} {:5e} GeV", egs_pid, energy);
@@ -284,4 +286,8 @@ namespace corsika {
fitout.close();
}
+ inline HEPEnergyType CONEXhybrid::getEnergyEM() const { return energy_em_; }
+
+ inline void CONEXhybrid::reset() { energy_em_ = 0_GeV; }
+
} // namespace corsika
diff --git a/corsika/detail/modules/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
index d58ac21d851ff8b19be258e09e9a63a500784b1a..eb8cbd2e5d6c6b1bc3b52b9999be634ffe079934 100644
--- a/corsika/detail/modules/pythia8/Decay.inl
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -55,7 +55,7 @@ namespace corsika::pythia8 {
Pythia8::Pythia::readString("Check:event = 1");
Pythia8::Pythia::readString("ProcessLevel:all = off");
- Pythia8::Pythia::readString("ProcessLevel:resonanceDecays = off");
+ Pythia8::Pythia::readString("ProcessLevel:resonanceDecays = on");
// making sure
setStable(Code::Pi0);
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
index a97265abc3a3aa442e9d373d579a37fe3da66412..217fe30b72c3fa27d1ad3794a4a6b1ec8852ffd5 100644
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ b/corsika/detail/modules/qgsjetII/Interaction.inl
@@ -66,15 +66,18 @@ namespace corsika::qgsjetII {
iTarget = targetA;
if (iTarget > maxMassNumber_ || iTarget <= 0) {
std::ostringstream txt;
- txt << "QgsjetII target outside range. iTarget=" << iTarget;
+ txt << "QgsjetII target outside range. Atarget=" << iTarget;
throw std::runtime_error(txt.str().c_str());
}
}
int iProjectile = 1;
if (is_nucleus(beamId)) {
iProjectile = Abeam;
- if (iProjectile > maxMassNumber_ || iProjectile <= 0)
- throw std::runtime_error("QgsjetII target outside range. ");
+ if (iProjectile > maxMassNumber_ || iProjectile <= 0) {
+ std::ostringstream txt;
+ txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
+ throw std::runtime_error(txt.str().c_str());
+ }
}
CORSIKA_LOG_DEBUG(
@@ -236,12 +239,16 @@ namespace corsika::qgsjetII {
if (is_nucleus(targetCode)) { // nucleus
targetMassNumber = get_nucleus_A(targetCode);
if (targetMassNumber > maxMassNumber_)
- throw std::runtime_error("QgsjetII target mass outside range.");
+ throw std::runtime_error(
+ "QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
} else {
- if (targetCode != Proton::code)
- throw std::runtime_error("QgsjetII Taget not possible.");
+ if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
+ throw std::runtime_error(
+ "QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
+ // here
}
- CORSIKA_LOG_DEBUG("Interaction: target qgsjetII code/A: ", targetMassNumber);
+ CORSIKA_LOG_DEBUG("Interaction: target qgsjetII code/A: {}", targetMassNumber);
int projectileMassNumber = 1; // "1" means "hadron"
QgsjetIIHadronType qgsjet_hadron_type =
@@ -249,7 +256,9 @@ namespace corsika::qgsjetII {
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
projectileMassNumber = projectile.getNuclearA();
if (projectileMassNumber > maxMassNumber_)
- throw std::runtime_error("QgsjetII projectile mass outside range.");
+ throw std::runtime_error(
+ "QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1);
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index a1abaa6abf994b7f07a165b4f42599c3b4aaa89f..ee29642e00bbd408c355a5c2c645abef312564ec 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -15,15 +15,112 @@
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <boost/filesystem.hpp>
+#include <boost/multi_array.hpp>
+
#include <algorithm>
+#include <cassert>
#include <functional>
#include <iostream>
+#include <fstream>
+#include <sstream>
#include <urqmd.hpp>
namespace corsika::urqmd {
- inline UrQMD::UrQMD() { ::urqmd::iniurqmdc8_(); }
+ inline UrQMD::UrQMD(boost::filesystem::path const& xs_file) {
+ readXSFile(xs_file);
+ ::urqmd::iniurqmdc8_();
+ }
+
+ inline CrossSectionType UrQMD::getTabulatedCrossSection(Code projectileCode,
+ Code targetCode,
+ HEPEnergyType labEnergy) const {
+ // translated to C++ from CORSIKA 7 subroutine cxtot_u
+
+ auto const kinEnergy = labEnergy - get_mass(projectileCode);
+
+ assert(kinEnergy >= HEPEnergyType::zero());
+
+ double const logKinEnergy = std::log10(kinEnergy * (1 / 1_GeV));
+ double const ye = std::max(10 * logKinEnergy + 10.5, 1.);
+ int const je = std::min(int(ye), int(xs_interp_support_table_.shape()[2] - 2));
+ std::array<double, 3> w;
+ w[2 - 1] = ye - je;
+ w[3 - 1] = w[2 - 1] * (w[2 - 1] - 1.) * .5;
+ w[1 - 1] = 1 - w[2 - 1] + w[3 - 1];
+ w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
+
+ int projectileIndex;
+ switch (projectileCode) {
+ case Code::Proton:
+ projectileIndex = 0;
+ break;
+ case Code::AntiProton:
+ projectileIndex = 1;
+ break;
+ case Code::Neutron:
+ projectileIndex = 2;
+ break;
+ case Code::AntiNeutron:
+ projectileIndex = 3;
+ break;
+ case Code::PiPlus:
+ projectileIndex = 4;
+ break;
+ case Code::PiMinus:
+ projectileIndex = 5;
+ break;
+ case Code::KPlus:
+ projectileIndex = 6;
+ break;
+ case Code::KMinus:
+ projectileIndex = 7;
+ break;
+ case Code::K0Short:
+ case Code::K0Long:
+ /* since K0Short and K0Long are treated the same, we can also add K0 and K0Bar
+ * to the list. This is a deviation from CORSIKA 7. */
+ case Code::K0:
+ case Code::K0Bar:
+ projectileIndex = 8;
+ break;
+ default: { // LCOV_EXCL_START since this can never happen due to canInteract
+ CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
+ return CrossSectionType::zero(); // LCOV_EXCL_STOP
+ }
+ }
+
+ int targetIndex;
+ switch (targetCode) {
+ case Code::Nitrogen:
+ targetIndex = 0;
+ break;
+ case Code::Oxygen:
+ targetIndex = 1;
+ break;
+ case Code::Argon:
+ targetIndex = 2;
+ break;
+ default:
+ std::stringstream ss;
+ ss << "UrQMD cross-section not tabluated for target " << targetCode;
+ throw std::runtime_error(ss.str().data());
+ }
+
+ auto result = CrossSectionType::zero();
+ for (int i = 0; i < 3; ++i) {
+ result +=
+ xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i];
+ }
+
+ CORSIKA_LOG_TRACE(
+ "UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
+ get_name(projectileCode), get_name(targetCode), labEnergy / 1_GeV, result);
+
+ return result;
+ }
inline CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
HEPEnergyType vLabEnergy,
@@ -63,20 +160,19 @@ namespace corsika::urqmd {
template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile
- inline CrossSectionType UrQMD::getCrossSection(TParticle const& vProjectile,
- Code vTargetCode) const {
- // TODO: return 0 for non-hadrons?
-
- auto const projectileCode = vProjectile.getPID();
- auto const projectileEnergyLab = vProjectile.getEnergy();
-
- if (projectileCode == Code::K0Long) {
- return 0.5 * (getCrossSection(Code::K0, vTargetCode, projectileEnergyLab) +
- getCrossSection(Code::K0Bar, vTargetCode, projectileEnergyLab));
+ inline CrossSectionType UrQMD::getCrossSection(TParticle const& projectile,
+ Code targetCode) const {
+ auto const projectileCode = projectile.getPID();
+
+ if (projectileCode == Code::Nucleus) {
+ /*
+ * unfortunately unavoidable at the moment until we have tools to get the actual
+ * inealstic cross-section from UrQMD
+ */
+ return CrossSectionType::zero();
}
- int const Ap = (projectileCode == Code::Nucleus) ? vProjectile.getNuclearA() : 1;
- return getCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
+ return getTabulatedCrossSection(projectileCode, targetCode, projectile.getEnergy());
}
inline bool UrQMD::canInteract(Code vCode) const {
@@ -85,10 +181,9 @@ namespace corsika::urqmd {
// only the usual long-lived species as input.
// TODO: Charmed mesons should be added to the list, too
- static Code const validProjectileCodes[] = {
- Code::Nucleus, Code::Proton, Code::AntiProton, Code::Neutron,
- Code::AntiNeutron, Code::PiPlus, Code::PiMinus, Code::KPlus,
- Code::KMinus, Code::K0, Code::K0Bar, Code::K0Long};
+ static std::array constexpr validProjectileCodes{
+ Code::Proton, Code::AntiProton, Code::Neutron, Code::AntiNeutron, Code::PiPlus,
+ Code::PiMinus, Code::KPlus, Code::KMinus, Code::K0Short, Code::K0Long};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
@@ -301,4 +396,41 @@ namespace corsika::urqmd {
return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
}
+ inline void UrQMD::readXSFile(boost::filesystem::path const& filename) {
+ boost::filesystem::ifstream file(filename, std::ios::in);
+
+ if (!file.is_open()) {
+ throw std::runtime_error(
+ filename.native() +
+ " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
+ }
+
+ std::string line;
+
+ std::getline(file, line);
+ std::stringstream ss(line);
+
+ char dummy;
+ int nTargets, nProjectiles, nSupports;
+ ss >> dummy >> nTargets >> nProjectiles >> nSupports;
+
+ decltype(xs_interp_support_table_)::extent_gen extents;
+ xs_interp_support_table_.resize(extents[nProjectiles][nTargets][nSupports]);
+
+ for (int i = 0; i < nTargets; ++i) {
+ for (int j = 0; j < nProjectiles; ++j) {
+ for (int k = 0; k < nSupports; ++k) {
+ std::getline(file, line);
+ std::stringstream s(line);
+ double energy, sigma;
+ s >> energy >> sigma;
+ xs_interp_support_table_[j][i][k] = sigma * 1_mb;
+ }
+
+ std::getline(file, line);
+ std::getline(file, line);
+ }
+ }
+ }
+
} // namespace corsika::urqmd
diff --git a/corsika/detail/setup/SetupStack.hpp b/corsika/detail/setup/SetupStack.hpp
index d31e940f6371826772f0c29a4915663a00397f28..e58ff5d9c7a7e88450147d8761e22585c1de8512 100644
--- a/corsika/detail/setup/SetupStack.hpp
+++ b/corsika/detail/setup/SetupStack.hpp
@@ -11,6 +11,7 @@
#include <corsika/framework/stack/CombinedStack.hpp>
#include <corsika/stack/GeometryNodeStackExtension.hpp>
#include <corsika/stack/NuclearStackExtension.hpp>
+#include <corsika/stack/WeightStackExtension.hpp>
#include <corsika/stack/history/HistorySecondaryProducer.hpp>
#include <corsika/stack/history/HistoryStackExtension.hpp>
diff --git a/corsika/detail/stack/WeightStackExtension.inl b/corsika/detail/stack/WeightStackExtension.inl
new file mode 100644
index 0000000000000000000000000000000000000000..719dc81a036c5d079317dc1d0bd10a549a3b0289
--- /dev/null
+++ b/corsika/detail/stack/WeightStackExtension.inl
@@ -0,0 +1,92 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <tuple>
+#include <utility>
+#include <vector>
+
+namespace corsika::weights {
+
+ // default version for particle-creation from input data
+ template <typename TParentStack>
+ inline void WeightDataInterface<TParentStack>::setParticleData(
+ std::tuple<double> const v) {
+ setWeight(std::get<0>(v));
+ }
+
+ template <typename TParentStack>
+ inline void WeightDataInterface<TParentStack>::setParticleData(
+ WeightDataInterface const& parent, std::tuple<double> const) {
+ setWeight(parent.getWeight()); // copy Weight from parent particle!
+ }
+
+ template <typename TParentStack>
+ inline void WeightDataInterface<TParentStack>::setParticleData() {
+ setWeight(1);
+ } // default weight
+
+ template <typename TParentStack>
+ inline void WeightDataInterface<TParentStack>::setParticleData(
+ WeightDataInterface const& parent) {
+ setWeight(parent.getWeight()); // copy Weight from parent particle!
+ }
+
+ template <typename TParentStack>
+ inline std::string WeightDataInterface<TParentStack>::asString() const {
+ return fmt::format("weight={}", getWeight());
+ }
+
+ template <typename TParentStack>
+ inline void WeightDataInterface<TParentStack>::setWeight(double const v) {
+ super_type::getStackData().setWeight(super_type::getIndex(), v);
+ }
+
+ template <typename TParentStack>
+ inline double WeightDataInterface<TParentStack>::getWeight() const {
+ return super_type::getStackData().getWeight(super_type::getIndex());
+ }
+
+ // definition of stack-data object to store geometry information
+
+ // these functions are needed for the Stack interface
+ inline void WeightData::clear() { weight_vector_.clear(); }
+
+ inline unsigned int WeightData::getSize() const { return weight_vector_.size(); }
+
+ inline unsigned int WeightData::getCapacity() const { return weight_vector_.size(); }
+
+ inline void WeightData::copy(int const i1, int const i2) {
+ weight_vector_[i2] = weight_vector_[i1];
+ }
+
+ inline void WeightData::swap(int const i1, int const i2) {
+ std::swap(weight_vector_[i1], weight_vector_[i2]);
+ }
+
+ // custom data access function
+ inline void WeightData::setWeight(int const i, double const v) {
+ weight_vector_[i] = v;
+ }
+
+ inline double WeightData::getWeight(int const i) const { return weight_vector_[i]; }
+
+ // these functions are also needed by the Stack interface
+ inline void WeightData::incrementSize() {
+ weight_vector_.push_back(1);
+ } // default weight
+
+ inline void WeightData::decrementSize() {
+ if (weight_vector_.size() > 0) { weight_vector_.pop_back(); }
+ }
+
+} // namespace corsika::weights
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index 99f0c97bda5de1f19615d7c9945d8b468052ae1f..bcab9ba975820b61060c5061aaea8dfbf31d3899 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -91,8 +91,10 @@ namespace corsika {
, output_(out)
, stack_(stack) {
CORSIKA_LOG_INFO(c8_ascii_);
+ CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(),
+ TTracking::getVersion());
if constexpr (TStackView::has_event) {
- CORSIKA_LOG_INFO(" - With full cascade HISTORY.");
+ CORSIKA_LOG_INFO("Stack - with full cascade HISTORY.");
}
}
//! \}
diff --git a/corsika/framework/core/Logging.hpp b/corsika/framework/core/Logging.hpp
index 2daea27c596be86c039f31f7a0b608e28d47ff4e..67f059c2fe6d1dc5d4b228e551f4ebbbcd1acf9e 100644
--- a/corsika/framework/core/Logging.hpp
+++ b/corsika/framework/core/Logging.hpp
@@ -7,7 +7,7 @@
*/
/**
- * @File Logging.hpp
+ * @file Logging.hpp
*
* CORSIKA8 logging utilities.
*
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index e133d68679c27bc597050a46d84c5bba2bfed70e..d0d17d264d0dbe8e7a7338ecce743482d0e700c1 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -10,11 +10,6 @@
@file ParticleProperties.hpp
Interface to particle properties
-
- The properties of all particles are saved in static and flat
- arrays. There is a enum corsika::Code to identify each
- particles, and each individual particles has its own static class,
- which can be used to retrieve its physical properties.
*/
#pragma once
@@ -29,21 +24,46 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
-/**
- * \file ParticleProperties.hpp
- *
- * The properties of all elementary particles are accessible here. The data
- * are taken from the Pythia ParticleData.xml file.
- *
- */
-
namespace corsika {
+
/**
- * @enum Code
- * The Code enum is the actual place to define CORSIKA 8 particle codes.
+ @defgroup Particles Particle Properties
+
+ The properties of all particles are saved in static and flat
+ arrays. There is a enum corsika::Code to identify each
+ particles, and each individual particles has its own static class,
+ which can be used to retrieve its physical properties.
+
+ The properties of all elementary particles are accessible here. The data
+ are taken from the Pythia ParticleData.xml file.
+
+ Particle data can be accessed via global function in namespace corsika, or via
+ static classes for each particle type. These classes all have the interface (example
+ for the class corsika::Electron):
+
+ @code{.cpp}
+ static constexpr Code code{Code::Electron};
+ static constexpr Code anti_code{Code::Positron};
+ static constexpr HEPMassType mass{corsika::get_mass(code)};
+ static constexpr ElectricChargeType charge{corsika::get_charge(code)};
+ static constexpr int charge_number{corsika::get_charge_number(code)};
+ static constexpr std::string_view name{corsika::get_name(code)};
+ static constexpr bool is_nucleus{corsika::is_nucleus(code)};
+ @endcode
+
+ The names, relations and properties of all particles known to CORSIKA 8 are listed
+ below.
+
+ @addtogroup Particles
+ @{
*/
+
+ /** The Code enum is the actual place to define CORSIKA 8 particle codes. */
enum class Code : int16_t;
+
+ /** Specifically for PDG ids */
enum class PDGCode : int32_t;
+
using CodeIntType = std::underlying_type<Code>::type;
using PDGCodeType = std::underlying_type<PDGCode>::type;
@@ -92,6 +112,9 @@ namespace corsika {
//! the output stream operator for human-readable particle codes
std::ostream& operator<<(std::ostream&, corsika::Code);
+
+ /** @}*/
+
} // namespace corsika
// data arrays, etc., as generated automatically
diff --git a/corsika/framework/geometry/QuantityVector.hpp b/corsika/framework/geometry/QuantityVector.hpp
index d9199dc57490bea7c0204c6256eddc7d4f80ef8c..ac7afc51e96e6a0cc4a4a59b8102c3accbb0c476 100644
--- a/corsika/framework/geometry/QuantityVector.hpp
+++ b/corsika/framework/geometry/QuantityVector.hpp
@@ -89,7 +89,7 @@ namespace corsika {
auto& operator-=(QuantityVector<TDimension> const& pQVec);
- auto& operator-() const;
+ auto operator-() const;
auto normalized() const;
diff --git a/corsika/framework/geometry/Vector.hpp b/corsika/framework/geometry/Vector.hpp
index fec720e4df43f33acb09aa536b6d6b482087cf5d..2a4b34a196a9363fc87926daeb9d5f5fcdb86076 100644
--- a/corsika/framework/geometry/Vector.hpp
+++ b/corsika/framework/geometry/Vector.hpp
@@ -131,7 +131,7 @@ namespace corsika {
auto& operator-=(Vector<TDimension> const& pVec);
- auto& operator-() const;
+ auto operator-() const;
auto normalized() const;
diff --git a/corsika/framework/process/BaseProcess.hpp b/corsika/framework/process/BaseProcess.hpp
index b8c96a54faa504a84c79272ce23f67534b6950b8..23e8667b20c9b92c9200135fcfe1b107c6e6c386 100644
--- a/corsika/framework/process/BaseProcess.hpp
+++ b/corsika/framework/process/BaseProcess.hpp
@@ -12,26 +12,25 @@
#include <type_traits>
-namespace corsika {
+//! @file BaseProcess.hpp
- class TDerived; // fwd decl
+namespace corsika {
/**
+ @ingroup Processes
+ @{
+
Each process in C8 must derive from BaseProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
- are of type BaseProcess<T>
+ are of type BaseProcess
- \todo rename BaseProcess into just Process
- \todo rename _BaseProcess, or find better alternative in FIXME
- ./Processes/AnalyticProcessors/ExecTime.h, see e.g. how this is done in
- ProcessSequence.hpp/make_sequence
+ @todo rename BaseProcess into just Process
*/
- class _BaseProcess {};
template <typename TDerived>
- struct BaseProcess : _BaseProcess {
+ struct BaseProcess {
protected:
friend TDerived;
@@ -39,20 +38,27 @@ namespace corsika {
// derived classes to be created, not
// BaseProcess itself
+ /** @name getRef Return reference to underlying type
+ @{
+ */
TDerived& ref() { return static_cast<TDerived&>(*this); }
const TDerived& ref() const { return static_cast<const TDerived&>(*this); }
+ //! @}
public:
- //! Default number of processes ist just one, obviously
+ static bool const is_process_sequence = false;
+ static bool const is_switch_process_sequence = false;
+
+ //! Default number of processes is just one, obviously
static unsigned int constexpr getNumberOfProcesses() { return 1; }
- // Base processor type for use in other template classes
+ //! Base processor type for use in other template classes
using process_type = TDerived;
};
/**
- * ProcessTraits specialization
- **/
+ is_process traits specialization to indicate inheritance from BaseProcess
+ */
template <typename TProcess>
struct is_process<
TProcess,
@@ -60,11 +66,17 @@ namespace corsika {
typename std::decay_t<TProcess>>>>
: std::true_type {};
+ /**
+ count_processes traits specialization to increase process count by one.
+ */
template <typename TProcess, int N>
- struct count_processes<TProcess, N,
- typename std::enable_if_t<is_process_v<TProcess> &&
- !is_process_sequence_v<TProcess>>> {
+ struct count_processes<
+ TProcess, N,
+ typename std::enable_if_t<is_process_v<std::decay_t<TProcess>> &&
+ !std::decay_t<TProcess>::is_process_sequence>> {
static unsigned int constexpr count = N + 1;
};
+ //! @}
+
} // namespace corsika
diff --git a/corsika/framework/process/BoundaryCrossingProcess.hpp b/corsika/framework/process/BoundaryCrossingProcess.hpp
index 1e9a4bdee4eb3de6de009e7a1274bddb3e967d2a..553b617bc136e1c128ff16cbde017186cf210f59 100644
--- a/corsika/framework/process/BoundaryCrossingProcess.hpp
+++ b/corsika/framework/process/BoundaryCrossingProcess.hpp
@@ -13,33 +13,49 @@
#include <type_traits>
-namespace corsika {
+#include <corsika/detail/framework/process/BoundaryCrossingProcess.hpp> // for extra traits, method/interface checking
- /*
- struct passepartout {
- template <typename T>
- operator T&();
+namespace corsika {
- template <typename T>
- operator T &&();
- };*/
+ /**
+ @ingroup Processes
+ @{
- template <typename TDerived>
- class BoundaryCrossingProcess : public BaseProcess<TDerived> {
+ Processes acting on the particles traversion from one volume into
+ another volume.
- /* static_assert(std::is_invocable_v<decltype(&TDerived<>::doBoundaryCrossing),
- TDerived&, passepartout>, "BoundaryCrossingProcess needs
- doBoundaryCrossing(TParticle, " "TParticle::node_type, TParticle::node_type)");*/
+ Create a new BoundaryCrossingProcess, e.g. for XYModel, via
+ @code{.cpp}
+ class XYModel : public BoundaryCrossingProcess<XYModel> {};
+ @endcode
- public:
- /**
- * This method is called when a particle crosses the boundary between the nodes
- * \p from and \p to.
- */
- template <typename TParticle>
- ProcessReturn doBoundaryCrossing(TParticle&,
+ and provide the necessary interface method:
+ @code{.cpp}
+ template <typename TParticle>
+ ProcessReturn XYModel::doBoundaryCrossing(TParticle& Particle,
typename TParticle::node_type const& from,
typename TParticle::node_type const& to);
+ @endcode
+
+ where Particle is the object to read particle data from a
+ Stack. The volume the particle is originating from is `from`, the
+ volume where it goes to is `to`.
+ */
+
+ template <typename TDerived>
+ class BoundaryCrossingProcess : public BaseProcess<TDerived> {
+ public:
};
+ /**
+ * ProcessTraits specialization to flag BoundaryProcess objects
+ **/
+ template <typename TProcess>
+ struct is_boundary_process<TProcess,
+ std::enable_if_t<std::is_base_of_v<
+ BoundaryCrossingProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>> : std::true_type {};
+
+ /** @} */
+
} // namespace corsika
diff --git a/corsika/framework/process/ContinuousProcess.hpp b/corsika/framework/process/ContinuousProcess.hpp
index 51d73d332fbe00ea3c5845608a0c582f224d5b4c..c549bc714da8e4c6b655afbb8690568a2f4c4e24 100644
--- a/corsika/framework/process/ContinuousProcess.hpp
+++ b/corsika/framework/process/ContinuousProcess.hpp
@@ -13,46 +13,62 @@
#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/detail/framework/process/ContinuousProcess.hpp> // for extra traits, method/interface checking
+
namespace corsika {
/**
+ @ingroup Processes
+ @{
+
Processes with continuous effects along a particle Trajectory
- The structural base type of a process object in a
- ProcessSequence. Both, the ProcessSequence and all its elements
- are of type ContinuousProcess<T>
+ 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.
*/
template <typename TDerived>
class ContinuousProcess : public BaseProcess<TDerived> {
- private:
- protected:
public:
- // here starts the interface part
- /**
- * Applies the effects of this ContinuousProcess on a Particle on a Track.
- *
- * Note, the stepLimit is a flag, if this particular process was responsible for the
- * track-length limit. This can be used by the process to trigger activity.
- *
- * \todo -> enforce TDerived to implement doContinuous...
- **/
- template <typename TParticle, typename TTrack>
- ProcessReturn doContinuous(TParticle&, TTrack const&, bool const stepLimit) const;
-
- /**
- * Calculates/returns a possible step length limitation of this continuousprocess.
- *
- *
- * \todo -> enforce TDerived to implement getMaxStepLength...
- **/
- template <typename TParticle, typename TTrack>
- LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
};
/**
- * ProcessTraits specialization
+ * ProcessTraits specialization to flag ContinuousProcess objects
**/
template <typename TProcess>
struct is_continuous_process<
@@ -61,10 +77,6 @@ namespace corsika {
typename std::decay_t<TProcess>>>>
: std::true_type {};
- template <typename TProcess, int N>
- struct count_continuous<TProcess, N,
- typename std::enable_if_t<is_continuous_process_v<TProcess>>> {
- enum { count = N + 1 };
- };
+ /** @} */
} // namespace corsika
diff --git a/corsika/framework/process/ContinuousProcessIndex.hpp b/corsika/framework/process/ContinuousProcessIndex.hpp
index 9d6e84d0f39baaa9a1e0ad7dbf369aba6a967fb3..e36bf2812a33c6f03e4058033585f846a871b897 100644
--- a/corsika/framework/process/ContinuousProcessIndex.hpp
+++ b/corsika/framework/process/ContinuousProcessIndex.hpp
@@ -11,9 +11,11 @@
namespace corsika {
/**
- * To index individual processes (continuous processes) inside a
- * ProcessSequence.
- *
+ @ingroup Processes
+
+ To index individual processes (continuous processes) inside a
+ ProcessSequence.
+
**/
class ContinuousProcessIndex {
diff --git a/corsika/framework/process/ContinuousProcessStepLength.hpp b/corsika/framework/process/ContinuousProcessStepLength.hpp
index 32b921857aaf6762c584ac5bf68b6885c019e038..575be9063d49a32b5853320dc9d25bac40a7b9d8 100644
--- a/corsika/framework/process/ContinuousProcessStepLength.hpp
+++ b/corsika/framework/process/ContinuousProcessStepLength.hpp
@@ -14,9 +14,11 @@
namespace corsika {
/**
- * 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 {
diff --git a/corsika/framework/process/DecayProcess.hpp b/corsika/framework/process/DecayProcess.hpp
index 64298dc3442176149f6cd7747219f7d2a9bac482..dc1e9a3ed3984570787b2a6ce36f3d88a20d90b7 100644
--- a/corsika/framework/process/DecayProcess.hpp
+++ b/corsika/framework/process/DecayProcess.hpp
@@ -11,15 +11,40 @@
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/detail/framework/process/DecayProcess.hpp> // for extra traits, method/interface checking
+
namespace corsika {
/**
+ @ingroup Processes
+ @{
+
Process decribing the decay of particles
- The structural base type of a process object in a
- ProcessSequence. Both, the ProcessSequence and all its elements
- are of type DecayProcess<T>
+ Create a new DecayProcess, e.g. for XYModel, via
+ @code
+ class XYModel : public DecayProcess<XYModel> {};
+ @endcode
+
+ and provide the two necessary interface methods
+ @code
+ template <typename TSecondaryView>
+ void XYModel::doDecay(TSecondaryView&);
+ template <typename TParticle>
+ TimeType getLifetime(TParticle const&)
+ @endcode
+
+ Where, of course, SecondaryView and Particle are the valid
+ classes to access particles 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.
+
+ SecondaryView allows to retrieve the properties of the projectile
+ particles, AND to store new particles (secondaries) which then
+ subsequently can be processes by SecondariesProcess. This is how
+ the output of decays can be studied right away.
*/
template <typename TDerived>
@@ -27,18 +52,29 @@ namespace corsika {
public:
using BaseProcess<TDerived>::ref;
- /// here starts the interface-definition part
- // -> enforce TDerived to implement DoDecay...
template <typename TParticle>
- void doDecay(TParticle&);
+ InverseTimeType getInverseLifetime(TParticle const& particle) {
- template <typename TParticle>
- TimeType getLifetime(TParticle const&);
+ // interface checking on TProcess1
+ static_assert(has_method_getLifetime_v<TDerived, TimeType, TParticle const&>,
+ "TDerived has no method with correct signature \"GrammageType "
+ "getInteractionLength(TParticle const&)\" required for "
+ "InteractionProcess<TDerived>. ");
- template <typename TParticle>
- InverseTimeType getInverseLifetime(TParticle const& particle) {
return 1. / ref().getLifetime(particle);
}
};
+ /**
+ * ProcessTraits specialization to flag DecayProcess objects
+ **/
+ template <typename TProcess>
+ struct is_decay_process<
+ TProcess,
+ std::enable_if_t<std::is_base_of_v<DecayProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>>
+ : std::true_type {};
+
+ /** @} */
+
} // namespace corsika
diff --git a/corsika/framework/process/InteractionCounter.hpp b/corsika/framework/process/InteractionCounter.hpp
index 143fa430e79fea07937c67117646e79919f61f06..e38722049eaa0b00c62a4c231000821a4960120f 100644
--- a/corsika/framework/process/InteractionCounter.hpp
+++ b/corsika/framework/process/InteractionCounter.hpp
@@ -14,9 +14,19 @@
namespace corsika {
/*!
+ @ingroup Processes
+ @{
+
* Wrapper around an InteractionProcess that fills histograms of the number
- * of calls to DoInteraction() binned in projectile energy (both in
+ * of calls to `doInteraction()` binned in projectile energy (both in
* lab and center-of-mass frame) and species
+ *
+ * Use by wrapping a normal InteractionProcess
+ * @code{.cpp}
+ * InteractionProcess collision1;
+ * InteractionClounter<collision1> counted_collision1;
+ * @endcode
+ *
*/
template <class TCountedProcess>
class InteractionCounter
@@ -25,12 +35,17 @@ namespace corsika {
public:
InteractionCounter(TCountedProcess& process);
+ //! wrapper around internall process doInteraction
template <typename TSecondaryView>
void doInteraction(TSecondaryView& view);
+ ///! returns internal process getInteractionLength
template <typename TParticle>
GrammageType getInteractionLength(TParticle const& particle) const;
+ /** returns the filles histograms
+ @return InteractionHistogram, which contains the histogram data
+ */
InteractionHistogram const& getHistogram() const;
private:
@@ -38,6 +53,8 @@ namespace corsika {
InteractionHistogram histogram_;
};
+ //! @}
+
} // namespace corsika
#include <corsika/detail/framework/process/InteractionCounter.inl>
diff --git a/corsika/framework/process/InteractionHistogram.hpp b/corsika/framework/process/InteractionHistogram.hpp
index ab0fee8fd299a36d3f150db1e8d0f1213a615315..78796f609e3a83316a52d818cf045425af7c507a 100644
--- a/corsika/framework/process/InteractionHistogram.hpp
+++ b/corsika/framework/process/InteractionHistogram.hpp
@@ -23,11 +23,27 @@
namespace corsika {
+ /** @ingroup Processes
+ @{
+
+ Class that creates and stores histograms of collisions
+ @f$dN/dE_{lab}@f$, @f$dN/d\sqrt{s}@f$ which is used by class
+ InteractionCounter
+
+ Histograms are of type boost::histogram
+
+ */
+
class InteractionHistogram {
static double constexpr lower_edge_cms = 1e3, upper_edge_cms = 1e17; // eV sqrt s
static double constexpr lower_edge_lab = 1e3, upper_edge_lab = 1e21; // eV lab
static unsigned int constexpr num_bins_lab = 18 * 10, num_bins_cms = 14 * 10;
+ /**
+ hist_type is a boost::histogram with two axes
+ - a growing PDG id axis
+ - a fixed logarithmic energy axis as configured by the user
+ */
using hist_type =
decltype(detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab));
@@ -36,17 +52,28 @@ namespace corsika {
public:
InteractionHistogram();
- //! fill both CMS and lab histograms at the same time
+ /**
+ fill both CMS and lab histograms at the same time
+ @param projectile_id corsika::Code of particle
+ @param lab_energy Energy in lab. frame
+ @param mass_target Mass of target particle
+ @param A if projectile_id is corsika::Nucleus : Mass of nucleus
+ @param Z if projectile_id is corsika::Nucleus : Charge of nucleus
+ */
void fill(Code projectile_id, HEPEnergyType lab_energy, HEPEnergyType mass_target,
int A = 0, int Z = 0);
+ //! return histogram in c.m.s. frame
hist_type const& CMSHist() const { return inthist_cms_; }
+ /// return histogram in laboratory frame
hist_type const& labHist() const { return inthist_lab_; }
InteractionHistogram& operator+=(InteractionHistogram const& other);
InteractionHistogram operator+(InteractionHistogram other) const;
};
+ /** @} */
+
} // namespace corsika
#include <corsika/detail/framework/process/InteractionHistogram.inl> // for implementation
diff --git a/corsika/framework/process/InteractionProcess.hpp b/corsika/framework/process/InteractionProcess.hpp
index f80480926c74b5f908e3fe7b4c85643a30cfeadd..6965457a7fe1a6c495e42b62943ebe87b3771702 100644
--- a/corsika/framework/process/InteractionProcess.hpp
+++ b/corsika/framework/process/InteractionProcess.hpp
@@ -11,34 +11,73 @@
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/detail/framework/process/InteractionProcess.hpp> // for extra traits, method/interface checking
+
namespace corsika {
/**
+ @ingroup Processes
+ @{
+
Process describing the interaction of particles
- The structural base type of a process object in a
- ProcessSequence. Both, the ProcessSequence and all its elements
- are of type InteractionProcess<T>
+ Create a new InteractionProcess, e.g. for XYModel, via
+ @code
+ class XYModel : public InteractionProcess<XYModel> {};
+ @endcode
+
+ and provide the two necessary interface methods
+ @code
+ template <typename TSecondaryView>
+ void XYModel::doInteraction(TSecondaryView&);
+
+ template <typename TParticle>
+ GrammageType XYModel::getInteractionLength(TParticle const&)
+ @endcode
+
+ Where, of course, SecondaryView and Particle are the valid
+ classes to access particles 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.
+
+ SecondaryView allows to retrieve the properties of the projectile
+ particles, AND to store new particles (secondaries) which then
+ subsequently can be processes by SecondariesProcess. This is how
+ the output of interactions can be studied right away.
*/
template <typename TDerived>
class InteractionProcess : public BaseProcess<TDerived> {
+
public:
using BaseProcess<TDerived>::ref;
- /// here starts the interface-definition part
- // -> enforce TDerived to implement DoInteraction...
template <typename TParticle>
- void doInteraction(TParticle&);
+ InverseGrammageType getInverseInteractionLength(TParticle const& particle) {
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle const&);
+ // interface checking on TProcess1
+ static_assert(
+ has_method_getInteractionLength_v<TDerived, GrammageType, TParticle const&>,
+ "TDerived has no method with correct signature \"GrammageType "
+ "getInteractionLength(TParticle const&)\" required for "
+ "InteractionProcess<TDerived>. ");
- template <typename TParticle>
- InverseGrammageType getInverseInteractionLength(TParticle const& particle) {
return 1. / ref().getInteractionLength(particle);
}
};
+ /**
+ * ProcessTraits specialization to flag InteractionProcess objects
+ **/
+ template <typename TProcess>
+ struct is_interaction_process<
+ TProcess, std::enable_if_t<
+ std::is_base_of_v<InteractionProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>>
+ : std::true_type {};
+
+ /** @} */
+
} // namespace corsika
diff --git a/corsika/framework/process/NullModel.hpp b/corsika/framework/process/NullModel.hpp
index 0e4746f31a924f4bf597f93eada980696d493efb..cdaa3753e0369574eb067e218096613f4c3a3c5d 100644
--- a/corsika/framework/process/NullModel.hpp
+++ b/corsika/framework/process/NullModel.hpp
@@ -8,19 +8,48 @@
#pragma once
-#include <corsika/framework/process/BaseProcess.hpp>
+//#include <corsika/framework/process/BaseProcess.hpp>
+#include <corsika/framework/process/ProcessTraits.hpp>
namespace corsika {
/**
- * Process that does nothing
+ @ingroup Processes
+ @{
+
+ Process that does nothing. It is not even derived from
+ BaseProcess. But it can be added to a ProcessSequence.
*/
- class NullModel : public BaseProcess<NullModel> {
+ class NullModel {
public:
NullModel() = default;
~NullModel() = default;
+
+ static bool const is_process_sequence = false;
+ static bool const is_switch_process_sequence = false;
+
+ //! Default number of processes is zero, obviously
+ static unsigned int constexpr getNumberOfProcesses() { return 0; }
+ };
+
+ /**
+ is_process traits specialization to indicate compatibility with BaseProcess
+ */
+ template <typename TNull>
+ struct is_process<
+ TNull, std::enable_if_t<std::is_base_of_v<NullModel, typename std::decay_t<TNull>>>>
+ : std::true_type {};
+
+ /**
+ count_processes traits specialization to increase process count by one.
+ */
+ template <int N>
+ struct count_processes<NullModel, N, void> {
+ static unsigned int constexpr count = N;
};
+ //! @}
+
} // namespace corsika
diff --git a/corsika/framework/process/ProcessReturn.hpp b/corsika/framework/process/ProcessReturn.hpp
index eef7058d09c171df53686fb4850a77d9e92c91c5..35c4096e53bdb8caba3727ff96c58bb159a361db 100644
--- a/corsika/framework/process/ProcessReturn.hpp
+++ b/corsika/framework/process/ProcessReturn.hpp
@@ -15,6 +15,9 @@
namespace corsika {
/**
+ @ingroup
+ @{
+
since in a process sequence many status updates can accumulate
for a single particle, this enum should define only bit-flags
that can be accumulated easily with "|="
@@ -59,4 +62,6 @@ namespace corsika {
return static_cast<int>(a & ProcessReturn::Interacted);
}
+ //! @}
+
} // namespace corsika
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 30c270bcbf30acaf0c8a539aada621eb0618886e..05191c248d5e1b17ceb5cfe39fc5dfe2fcabe4f6 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -9,7 +9,7 @@
#pragma once
/**
- * \file ProcessSequence.hpp
+ * @file ProcessSequence.hpp
*/
#include <corsika/framework/process/BaseProcess.hpp>
@@ -28,68 +28,137 @@
namespace corsika {
+ /**
+ count_processes traits specialization to increase process count by
+ getNumberOfProcesses(). This is used to statically count processes in the sequence
+ **/
template <typename TProcess, int N>
- struct count_continuous<TProcess, N,
- typename std::enable_if_t<is_process_sequence_v<TProcess>>> {
- static unsigned int constexpr count =
- N + std::decay_t<TProcess>::getNumberOfProcesses();
- };
-
- template <typename TProcess, int N>
- struct count_processes<TProcess, N,
- typename std::enable_if_t<is_process_v<TProcess> &&
- is_process_sequence_v<TProcess>>> {
+ struct count_processes<
+ TProcess, N,
+ typename std::enable_if_t<is_process_v<std::decay_t<TProcess>> &&
+ std::decay_t<TProcess>::is_process_sequence>> {
static unsigned int constexpr count =
N + std::decay_t<TProcess>::getNumberOfProcesses();
};
/**
- *
- * Definition of a static process list/sequence
- *
- * A compile time static list of processes. The compiler will
- * generate a new type based on template logic containing all the
- * elements provided by the user.
- *
- * TProcess1 and TProcess2 must both be derived from BaseProcess,
- * and are both references if possible (lvalue), otherwise (rvalue)
- * they are just classes. This allows us to handle both, rvalue as
- * well as lvalue Processes in the ProcessSequence.
- *
- * (The sequence, and the processes use the CRTP, curiously recurring template
- * pattern).
- *
- * Template parameters:
- * - TProcess1 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
- * - TProcess2 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
- * - ProcessIndexOffset, IndexOfProcess1, IndexOfProcess2 are to count and index each
- *ContinuousProcess in the entire process-chain
- **/
+ @defgroup Processes Physics Processes and Modules
+
+ Physics processes in CORSIKA 8 are clustered in ProcessSequence and
+ SwitchProcessSequence containers. The former is a mere (ordered) collection, while
+ the latter has the option to switch between two alternative ProcessSequences.
+
+ Depending on the type of data to act on and on the allowed actions of processes there
+ are several interface options:
+ - InteractionProcess
+ - DecayProcess
+ - ContinuousProcess
+ - StackProcess
+ - SecondariesProcess
+ - BoundaryCrossingProcess
+
+ And all processes (including ProcessSequence and SwitchProcessSequence) are derived
+ from BaseProcess.
+
+ Processes of any type (e.g. p1, p2, p3,...) can be assembled into a ProcessSequence
+ using the `make_sequence` factory function.
+
+ @code{.cpp}
+ auto sequence1 = make_sequence(p1, p2, p3);
+ auto sequence2 = make_sequence(p4, p5, p6, p7);
+ auto sequence3 = make_sequence(sequence1, sequemce2, p8, p9);
+ @endcode
+
+ Note, if the order of processes
+ matters, the order of occurence
+ in the ProcessSequence determines
+ the executiion order.
+
+ SecondariesProcess alyways act on
+ new secondaries produced (i.e. in
+ InteractionProcess and
+ DecayProcess) in the scope of
+ their ProcessSequence. For
+ example if i1 and i2 are
+ InteractionProcesses and s1 is a
+ SecondariesProcess, then
+
+ @code{.cpp}
+ auto sequence = make_sequence(i1, make_sequence(i2, s1))
+ @endcode
+
+ will result in s1 acting only on
+ the particles produced by i2 and
+ not by i1. This can be very
+ useful, e.g. to fine tune thinning.
+
+ A special type of ProcessSequence
+ is SwitchProcessSequence, which
+ has two branches and a functor
+ that can select between these two
+ branches.
+
+ @code{.cpp}
+ auto sequence = make_switch(sequence1, sequence2, selector);
+ @endcode
+
+ where the only requirement to
+ `selector` is that it
+ provides a `SwitchResult operator()(Particle const& particle) const` method. Thus,
+ based on the dynamic properties
+ of `particle` the functor
+ can make its decision. This is
+ clearly important for switching
+ between low-energy and
+ high-energy models, but not
+ limited to this. The selection
+ can even be done with a lambda
+ function.
+
+
+ @class ProcessSequence
+ @ingroup Processes
+
+ Definition of a static process list/sequence
+
+ A compile time static list of processes. The compiler will
+ generate a new type based on template logic containing all the
+ elements provided by the user.
+
+ TProcess1 and TProcess2 must both be derived from BaseProcess,
+ and are both references if possible (lvalue), otherwise (rvalue)
+ they are just classes. This allows us to handle both, rvalue as
+ well as lvalue Processes in the ProcessSequence.
+
+ (For your potential interest,
+ the static version of the
+ ProcessSequence and all Process
+ types are based on the CRTP C++
+ design pattern)
+
+ Template parameters:
+ @tparam TProcess1 is of type BaseProcess, either a dedicatd process, or a
+ ProcessSequence
+ @tparam TProcess2 is of type BaseProcess, either a dedicatd process, or a
+ ProcessSequence
+ @tparam IndexFirstProcess to count and index each Process in the entire
+ process-chain. The offset is the starting value for this ProcessSequence
+ @tparam IndexOfProcess1 index of TProcess1 (counting of Process)
+ @tparam IndexOfProcess2 index of TProcess2 (counting of Process)
+ **/
template <typename TProcess1, typename TProcess2 = NullModel,
int ProcessIndexOffset = 0,
- int IndexOfProcess1 = count_processes<
- TProcess1, count_processes<TProcess2, ProcessIndexOffset>::count>::count,
- int IndexOfProcess2 = count_processes<TProcess2, ProcessIndexOffset>::count>
+ int IndexOfProcess1 = corsika::count_processes<
+ TProcess1,
+ corsika::count_processes<TProcess2, ProcessIndexOffset>::count>::count,
+ int IndexOfProcess2 =
+ corsika::count_processes<TProcess2, ProcessIndexOffset>::count>
class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> {
using process1_type = typename std::decay_t<TProcess1>;
using process2_type = typename std::decay_t<TProcess2>;
- static bool constexpr t1ProcSeq = is_process_sequence_v<process1_type>;
- static bool constexpr t2ProcSeq = is_process_sequence_v<process2_type>;
-
- static bool constexpr t1SwitchProcSeq = is_switch_process_sequence_v<process1_type>;
- static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<process2_type>;
-
- // make sure only BaseProcess types TProcess1/2 are passed
- static_assert(is_process_v<process1_type>,
- "can only use process derived from BaseProcess in "
- "ProcessSequence, for Process 1");
- static_assert(is_process_v<process2_type>,
- "can only use process derived from BaseProcess in "
- "ProcessSequence, for Process 2");
-
public:
// resource management
ProcessSequence() = delete; // only initialized objects
@@ -98,19 +167,19 @@ namespace corsika {
ProcessSequence& operator=(ProcessSequence const&) = default;
~ProcessSequence() = default;
+ static bool const is_process_sequence = true;
+
/**
- * Only valid user constructor will create fully initialized object
- *
- * ProcessSequence supports and encourages move semantics. You can
- * use object, l-value references or r-value references to
- * construct sequences.
- *
- * \param in_A process/list A
- * \param in_A process/list B
+ Only valid user constructor will create fully initialized object
+
+ ProcessSequence supports and encourages move semantics. You can
+ use object, l-value references or r-value references to
+ construct sequences.
+
+ @param in_A BaseProcess or switch/process list
+ @param in_B BaseProcess or switch/process list
**/
- ProcessSequence(TProcess1 in_A, TProcess2 in_B)
- : A_(in_A)
- , B_(in_B) {}
+ ProcessSequence(TProcess1 in_A, TProcess2 in_B);
template <typename TParticle>
ProcessReturn doBoundaryCrossing(TParticle& particle,
@@ -147,7 +216,7 @@ namespace corsika {
* The maximum allowed step length is the minimum of the allowed track lenght over all
* ContinuousProcess-es in the ProcessSequence.
*
- * \return: ContinuousProcessStepLength which contains the step length itself in
+ * @return ContinuousProcessStepLength which contains the step length itself in
* LengthType, and a unique identifier of the related ContinuousProcess.
**/
@@ -187,6 +256,11 @@ namespace corsika {
**/
static unsigned int constexpr getNumberOfProcesses() { return numberOfProcesses_; }
+#ifdef CORSIKA_UNIT_TESTING
+ TProcess1 getProcess1() const { return A_; }
+ TProcess2 getProcess2() const { return B_; }
+#endif
+
private:
TProcess1 A_; /// process/list A, this is a reference, if possible
TProcess2 B_; /// process/list B, this is a reference, if possible
@@ -195,33 +269,36 @@ namespace corsika {
};
/**
- * Factory function to create ProcessSequence
- *
- * to construct ProcessSequences in a flexible and dynamic way the
- * `sequence` factory functions are provided
- *
- * Any objects of type
- * - BaseProcess,
- * - ContinuousProcess, and
- * - Interaction/DecayProcess,
- * - StackProcess,
- * - SecondariesProcess
- * can be assembled into a ProcessSequence, all
- * combinatorics are allowed.
-
- * The sequence function checks that all its arguments are all of
- * types derived from BaseProcess. Also the ProcessSequence itself
- * is derived from type BaseProcess
- *
- * \param vA needs to derive from BaseProcess or ProcessSequence
- * \param vB paramter-pack, needs to derive BaseProcess or ProcessSequence
- *
- **/
-
+ @fn make_sequence
+ @ingroup Processes
+
+ Factory function to create a ProcessSequence
+
+ to construct ProcessSequences in a flexible and dynamic way the
+ `sequence` factory functions are provided
+
+ Any objects of type
+ - BaseProcess
+ - ContinuousProcess and
+ - InteractionProcess/DecayProcess
+ - StackProcess
+ - SecondariesProcess
+ - BoundaryCrossingProcess
+
+ can be assembled into a ProcessSequence, all
+ combinatorics are allowed.
+
+ The sequence function checks that all its arguments are all of
+ types derived from BaseProcess. Also the ProcessSequence itself
+ is derived from type BaseProcess
+
+ @tparam TProcesses parameter pack with objects of type BaseProcess
+ @tparam TProcess1 another BaseProcess
+ @param vA needs to derive from BaseProcess
+ @param vB paramter-pack, needs to derive BaseProcess
+ */
template <typename... TProcesses, typename TProcess1>
- typename std::enable_if_t<
- is_process_v<typename std::decay_t<TProcess1>>,
- ProcessSequence<TProcess1, decltype(make_sequence(std::declval<TProcesses>()...))>>
+ ProcessSequence<TProcess1, decltype(make_sequence(std::declval<TProcesses>()...))>
make_sequence(TProcess1&& vA, TProcesses&&... vBs) {
return ProcessSequence<TProcess1,
decltype(make_sequence(std::declval<TProcesses>()...))>(
@@ -229,47 +306,40 @@ namespace corsika {
}
/**
- * Factory function to create ProcessSequence
- *
- * specialization for two input objects (no paramter pack in vB).
- *
- * \param vA needs to derive from BaseProcess or ProcessSequence
- * \param vB needs to derive BaseProcess or ProcessSequence
- **/
+ @fn make_sequence
+ @ingroup Processes
+
+ Factory function to create ProcessSequence
+
+ specialization for two input objects (no paramter pack in vB).
+
+ @tparam TProcess1 another BaseProcess
+ @tparam TProcess2 another BaseProcess
+ @param vA needs to derive from BaseProcess
+ @param vB needs to derive BaseProcess
+ */
template <typename TProcess1, typename TProcess2>
- typename std::enable_if_t<is_process_v<typename std::decay_t<TProcess1>> &&
- is_process_v<typename std::decay_t<TProcess2>>,
- ProcessSequence<TProcess1, TProcess2>>
- make_sequence(TProcess1&& vA, TProcess2&& vB) {
+ ProcessSequence<TProcess1, TProcess2> make_sequence(TProcess1&& vA, TProcess2&& vB) {
return ProcessSequence<TProcess1, TProcess2>(vA, vB);
}
/**
- * Factory function to create ProcessSequence from a single BaseProcess
- *
- * also allow a single Process in ProcessSequence, accompany by
- * `NullModel`
- *
- * \param vA needs to derive from BaseProcess or ProcessSequence
- **/
+ @fn make_sequence
+ @ingroup Processes
+
+ Factory function to create ProcessSequence from a single BaseProcess
+
+ also allow a single Process in ProcessSequence, accompany by
+ `NullModel`
+
+ @tparam TProcess1 another BaseProcess
+ @param vA needs to derive from BaseProcess
+ */
template <typename TProcess>
- typename std::enable_if_t<is_process_v<typename std::decay_t<TProcess>>,
- ProcessSequence<TProcess, NullModel>>
- make_sequence(TProcess&& vA) {
+ ProcessSequence<TProcess, NullModel> make_sequence(TProcess&& vA) {
return ProcessSequence<TProcess, NullModel>(vA, NullModel());
}
- /**
- * traits marker to identify objectas ProcessSequence
- **/
- template <typename TProcess1, typename TProcess2>
- struct is_process_sequence<ProcessSequence<TProcess1, TProcess2>> : std::true_type {
- // only switch on for BaseProcesses
- template <typename std::enable_if_t<
- is_process_v<TProcess1> && is_process_v<TProcess2>, int>>
- is_process_sequence() {}
- };
-
} // namespace corsika
#include <corsika/detail/framework/process/ProcessSequence.inl>
diff --git a/corsika/framework/process/ProcessTraits.hpp b/corsika/framework/process/ProcessTraits.hpp
index b08318093e7be4147e5408661b16c8feb1c3921e..e675fb1717e0c848b91b78ea6b3e0730d2063cf4 100644
--- a/corsika/framework/process/ProcessTraits.hpp
+++ b/corsika/framework/process/ProcessTraits.hpp
@@ -9,7 +9,7 @@
#pragma once
/**
- * \file ProcessTraits.hpp
+ * @file ProcessTraits.hpp
*/
#include <type_traits>
@@ -19,7 +19,7 @@ namespace corsika {
/**
* A traits marker to identify BaseProcess, thus any type of process
*/
- template <typename TProcess, typename Enable = void>
+ template <typename TProcess, typename TEnable = void>
struct is_process : std::false_type {};
template <typename TProcess>
@@ -35,23 +35,49 @@ namespace corsika {
bool constexpr is_continuous_process_v = is_continuous_process<TProcess>::value;
/**
- * A traits marker to track which BaseProcess is also a ProcessSequence
- **/
- template <typename TClass>
- struct is_process_sequence : std::false_type {};
+ * A traits marker to identify DecayProcess
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_decay_process : std::false_type {};
- template <typename TClass>
- bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value;
+ template <typename TProcess>
+ bool constexpr is_decay_process_v = is_decay_process<TProcess>::value;
/**
- * A traits marker to identiy a BaseProcess that is also SwitchProcessesSequence
- **/
+ * A traits marker to identify StackProcess
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_stack_process : std::false_type {};
- template <typename TClass>
- struct is_switch_process_sequence : std::false_type {};
+ template <typename TProcess>
+ bool constexpr is_stack_process_v = is_stack_process<TProcess>::value;
- template <typename TClass>
- bool constexpr is_switch_process_sequence_v = is_switch_process_sequence<TClass>::value;
+ /**
+ * A traits marker to identify SecondariesProcess
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_secondaries_process : std::false_type {};
+
+ template <typename TProcess>
+ bool constexpr is_secondaries_process_v = is_secondaries_process<TProcess>::value;
+
+ /**
+ * A traits marker to identify BoundaryProcess
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_boundary_process : std::false_type {};
+
+ template <typename TProcess>
+ bool constexpr is_boundary_process_v = is_boundary_process<TProcess>::value;
+
+ /**
+ * A traits marker to identify InteractionProcess
+ */
+ template <typename TProcess, typename Enable = void>
+ struct is_interaction_process : std::false_type {};
+
+ template <typename TProcess>
+ bool constexpr is_interaction_process_v = is_interaction_process<TProcess>::value;
/**
* A traits marker to identify ProcessSequence that contain a StackProcess
@@ -62,14 +88,6 @@ namespace corsika {
template <typename TClass>
bool constexpr contains_stack_process_v = contains_stack_process<TClass>::value;
- /**
- * traits class to count ContinuousProcess-es, general version
- **/
- template <typename TProcess, int N = 0, typename Enable = void>
- struct count_continuous {
- static unsigned int constexpr count = N;
- };
-
/**
* traits class to count any type of Process, general version
**/
@@ -78,4 +96,33 @@ namespace corsika {
static unsigned int constexpr count = N;
};
+ namespace detail {
+
+ /**
+ Helper traits class (partial) for static compile time checking.
+
+ Note, this is a poor replacement for C++20 concepts... they are
+ eagerly awaited!
+
+ It defines the default body of a generic test function returning
+ std::false_type.
+
+ In addition it defines the pattern for class-method matching with a
+ return type TReturn and function arguments TArgs... . Right now
+ both method signatures, "const" and "not const", are matched.
+ */
+ template <typename TReturn, typename... TArgs>
+ struct has_method_signature {
+
+ // the non-const version
+ template <class T>
+ static std::true_type testSignature(TReturn (T::*)(TArgs...));
+
+ // the const version
+ template <class T>
+ static std::true_type testSignature(TReturn (T::*)(TArgs...) const);
+ };
+
+ } // namespace detail
+
} // namespace corsika
diff --git a/corsika/framework/process/SecondariesProcess.hpp b/corsika/framework/process/SecondariesProcess.hpp
index daa45cc396f817e623ef107c9fdffd24099b8cbe..975b1e7851174b3be25e43a99ef06c5155a2912c 100644
--- a/corsika/framework/process/SecondariesProcess.hpp
+++ b/corsika/framework/process/SecondariesProcess.hpp
@@ -11,24 +11,46 @@
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/detail/framework/process/SecondariesProcess.hpp> // for extra traits, method/interface checking
+
namespace corsika {
/**
- Process that modifies a list of secondaries of other processes
+ @ingroup Processes
+ @{
+
+ Processes acting on the secondaries produced by other processes.
+
+ Create a new SecondariesProcess, e.g. for XYModel, via
+ @code{.cpp}
+ class XYModel : public SecondariesProcess<XYModel> {};
+ @endcode
- The structural base type of a process object in a
- ProcessSequence. Both, the ProcessSequence and all its elements
- are of type SecondariesProcess<T>
+ and provide the necessary interface method:
+ @code{.cpp}
+ template <typename TStackView>
+ void doSecondaries(TStackView& StackView);
+ @endcode
+ where StackView is an object that can store secondaries on a
+ stack and also iterate over these secondaries.
*/
template <typename TDerived>
class SecondariesProcess : public BaseProcess<TDerived> {
public:
- /// here starts the interface-definition part
- // -> enforce TDerived to implement DoSecondaries...
- template <typename TSecondaries>
- void doSecondaries(TSecondaries&);
};
+ /**
+ * ProcessTraits specialization to flag SecondariesProcess objects
+ **/
+ template <typename TProcess>
+ struct is_secondaries_process<
+ TProcess, std::enable_if_t<
+ std::is_base_of_v<SecondariesProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>>
+ : std::true_type {};
+
+ //! @}
+
} // namespace corsika
diff --git a/corsika/framework/process/StackProcess.hpp b/corsika/framework/process/StackProcess.hpp
index 7bf23c82c4b5e5a380e64ee8a3eb8632c0d2be2b..773265f49e102b63ca447191a25040c555179c37 100644
--- a/corsika/framework/process/StackProcess.hpp
+++ b/corsika/framework/process/StackProcess.hpp
@@ -11,15 +11,37 @@
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/detail/framework/process/StackProcess.hpp> // for extra traits, method/interface checking
+
namespace corsika {
/**
+ @ingroup Processes
+ @{
+
Process to act on the entire particle stack
- The structural base type of a process object in a
- ProcessSequence. Both, the ProcessSequence and all its elements
- are of type StackProcess<T>
+ Create a new StackProcess, e.g. for XYModel, via
+ @code{.cpp}
+ class XYModel : public StackProcess<XYModel> {};
+ @endcode
+
+ and provide the necessary interface method
+ @code{.cpp}
+ template <typename TStack>
+ void XYModel::doStack(TStack&);
+ @endcode
+
+ Where, of course, Stack is the valid
+ class to access particles on the Stack. This methods does
+ not need to be templated, they could use the types
+ e.g. corsika::setup::Stack directly -- but by the cost of
+ loosing all flexibility otherwise provided.
+ A StackProcess has only one constructor `StackProcess::StackProcess(unsigned int
+ const nStep)` where nStep is the number of steps of the cascade stepping after which
+ the stack process should be run. Good values are on the order of 1000, which will not
+ compromise run time in the end, but provide all the benefits of the StackProcess.
*/
template <typename TDerived>
@@ -31,19 +53,20 @@ namespace corsika {
StackProcess(const unsigned int nStep)
: nStep_(nStep) {}
- /// here starts the interface-definition part
- // -> enforce TDerived to implement DoStack...
- template <typename TStack>
- void doStack(TStack&);
+ static bool const is_stack_process = true;
+ //! return the current Cascade step counter
int getStep() const { return iStep_; }
+
+ //! check if current step is where StackProcess should be executed, this also
+ //! increases the internal step counter implicitly
bool checkStep() { return !((++iStep_) % nStep_); }
private:
/**
@name The number of "steps" during the cascade processing after
which this StackProcess is going to be executed. The logic is
- "fIStep modulo fNStep"
+ "iStep_ modulo nStep_"
@{
*/
unsigned int nStep_ = 0;
@@ -51,4 +74,16 @@ namespace corsika {
//! @}
};
+ /**
+ * ProcessTraits specialization to flag StackProcess objects
+ **/
+ template <typename TProcess>
+ struct is_stack_process<
+ TProcess,
+ std::enable_if_t<std::is_base_of_v<StackProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>>>
+ : std::true_type {};
+
+ //! @}
+
} // namespace corsika
diff --git a/corsika/framework/process/SwitchProcessSequence.hpp b/corsika/framework/process/SwitchProcessSequence.hpp
index 22da0facc5c5df27f38e03fadf02db014138cd60..09edeea43e551f7f37b51e8b105c78cc92ba6f6c 100644
--- a/corsika/framework/process/SwitchProcessSequence.hpp
+++ b/corsika/framework/process/SwitchProcessSequence.hpp
@@ -9,7 +9,7 @@
#pragma once
/**
- * \file SwitchProcessSequence.hpp
+ @file SwitchProcessSequence.hpp
**/
#include <corsika/framework/process/BaseProcess.hpp>
@@ -31,26 +31,24 @@
namespace corsika {
/**
- * enum for the process switch selection: identify if First or
- * Second process branch should be used.
- **/
- enum class SwitchResult { First, Second };
+ @ingroup Processes
+ @{
- /**
Class to switch between two process branches
A compile-time static list of processes that uses an internal
- TSelect class to switch between different versions of processes
+ TCondition class to switch between different versions of processes
(or process sequence).
- TProcess1 and TProcess2 must be derived from BaseProcess and are
+ TSequence and USequence must be derived from BaseProcess and are
both references if possible (lvalue), otherwise (rvalue) they are
just classes. This allows us to handle both, rvalue as well as
lvalue Processes in the SwitchProcessSequence.
+ Please use the `corsika::make_select(condition, sequence, alt_sequence)`
+ factory function for best results.
- TSelect has to implement a `operator()(const Particle&)` and has to
- return either SwitchResult::First or SwitchResult::Second. Note:
- TSelect may absolutely also use random numbers to sample between
+ TCondition has to implement a `bool operator()(Particle const&)`. Note:
+ TCondition may absolutely also use random numbers to sample between
its results. This can be used to achieve arbitrarily smooth
transition or mixtures of processes.
@@ -59,28 +57,30 @@ namespace corsika {
particle stack and not on indiviidual particles.
Template parameters:
- - TProcess1 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
- - TProcess2 is of type BaseProcess, either a dedicatd process, or a ProcessSequence
- - IndexFirstProcess, IndexOfProcess1, IndexOfProcess2 are to count and index each
- ContinuousProcess in the entire process-chain
-
- See also class \sa ProcessSequence
+ @tparam TCondition selector functor/function
+ @tparam TSequence is of type BaseProcess, either a dedicatd process, or a
+ ProcessSequence
+ @tparam USequence is of type BaseProcess, either a dedicatd process, or a
+ ProcessSequence
+ @tparam IndexFirstProcess to count and index each Process in the entire
+ process-chain
+ @tparam IndexOfProcess1 index of TSequence (counting of Process)
+ @tparam IndexOfProcess2 index of USequence (counting of Process)
+
+ See also class ProcessSequence.
**/
- template <typename TProcess1, typename TProcess2, typename TSelect,
+ template <typename TCondition, typename TSequence, typename USequence,
int IndexFirstProcess = 0,
- int IndexOfProcess1 = count_processes<TProcess1, IndexFirstProcess>::count,
- int IndexOfProcess2 = count_processes<TProcess2, IndexOfProcess1>::count>
+ int IndexOfProcess1 = count_processes<TSequence, IndexFirstProcess>::count,
+ int IndexOfProcess2 = count_processes<USequence, IndexOfProcess1>::count>
class SwitchProcessSequence
- : public BaseProcess<SwitchProcessSequence<TProcess1, TProcess2, TSelect>> {
+ : public BaseProcess<SwitchProcessSequence<TCondition, TSequence, USequence>> {
- using process1_type = typename std::decay_t<TProcess1>;
- using process2_type = typename std::decay_t<TProcess2>;
+ using process1_type = typename std::decay_t<TSequence>;
+ using process2_type = typename std::decay_t<USequence>;
- static bool constexpr t1ProcSeq = is_process_sequence_v<process1_type>;
- static bool constexpr t2ProcSeq = is_process_sequence_v<process2_type>;
-
- // make sure only BaseProcess types TProcess1/2 are passed
+ // make sure only BaseProcess types TSequence/2 are passed
static_assert(is_process_v<process1_type>,
"can only use process derived from BaseProcess in "
"SwitchProcessSequence, for Process 1");
@@ -88,13 +88,13 @@ namespace corsika {
"can only use process derived from BaseProcess in "
"SwitchProcessSequence, for Process 2");
- // make sure none of TProcess1/2 is a StackProcess
+ // make sure none of TSequence/2 is a StackProcess
static_assert(!std::is_base_of_v<StackProcess<process1_type>, process1_type>,
"cannot use StackProcess in SwitchProcessSequence, for Process 1");
static_assert(!std::is_base_of_v<StackProcess<process2_type>, process2_type>,
"cannot use StackProcess in SwitchProcessSequence, for Process 2");
- // if TProcess1/2 are already ProcessSequences, make sure they do not contain
+ // if TSequence/2 are already ProcessSequences, make sure they do not contain
// any StackProcess
static_assert(!contains_stack_process_v<process1_type>,
"cannot use StackProcess in SwitchProcessSequence, remove from "
@@ -111,25 +111,29 @@ namespace corsika {
SwitchProcessSequence& operator=(SwitchProcessSequence const&) = default;
~SwitchProcessSequence() = default;
+ static bool const is_process_sequence = true;
+ static bool const is_switch_process_sequence = true;
+
/**
- * Only valid user constructor will create fully initialized object
- *
- * SwitchProcessSequence supports and encourages move semantics. You can
- * use object, l-value references or r-value references to
- * construct sequences.
- *
- * \param in_A process branch A
- * \param in_A process branch B
- * \param sel functor to swtich between branch A and B
+ Only valid user constructor will create fully initialized object
+
+ SwitchProcessSequence supports and encourages move semantics. You can
+ use object, l-value references or r-value references to
+ construct sequences.
+
+ @param sel functor to switch between branch A and B
+ @param in_A process branch A
+ @param in_A process branch B
**/
- SwitchProcessSequence(TProcess1 in_A, TProcess2 in_B, TSelect sel)
+ SwitchProcessSequence(TCondition sel, TSequence in_A, USequence in_B)
: select_(sel)
, A_(in_A)
, B_(in_B) {}
- template <typename TParticle, typename TVTNType>
- ProcessReturn doBoundaryCrossing(TParticle& particle, TVTNType const& from,
- TVTNType const& to);
+ template <typename TParticle>
+ ProcessReturn doBoundaryCrossing(TParticle& particle,
+ typename TParticle::node_type const& from,
+ typename TParticle::node_type const& to);
template <typename TParticle, typename TTrack>
ProcessReturn doContinuous(TParticle& particle, TTrack& vT,
@@ -174,35 +178,42 @@ namespace corsika {
**/
static unsigned int constexpr getNumberOfProcesses() { return numberOfProcesses_; }
+#ifdef CORSIKA_UNIT_TESTING
+ TCondition getCondition() const { return select_; }
+ TSequence getSequence() const { return A_; }
+ USequence getAltSequence() const { return B_; }
+#endif
+
private:
- TSelect select_; /// selector functor to switch between branch a and b, this is a
- /// reference, if possible
+ TCondition select_; /// selector functor to switch between branch a and b, this is a
+ /// reference, if possible
- TProcess1 A_; /// process branch a, this is a reference, if possible
- TProcess2 B_; /// process branch b, this is a reference, if possible
+ TSequence A_; /// process branch a, this is a reference, if possible
+ USequence B_; /// process branch b, this is a reference, if possible
static unsigned int constexpr numberOfProcesses_ = IndexOfProcess2; // static counter
};
/**
- *
- * the functin `make_select(proc1,proc1,selector)` assembles many
- * BaseProcesses, and ProcessSequences into a SwitchProcessSequence,
- * all combinatorics must be allowed, this is why we define a macro
- * to define all combinations here:
- *
- *
- * Both, Processes1 and Processes2, must derive from BaseProcesses
+ the functin `make_select(select, proc1, proc1)` assembles many
+ BaseProcesses, and ProcessSequences into a SwitchProcessSequence,
+ all combinatorics are allowed.
+
+ @param selector must provide `bool operator()(Particle const&) const`
+ @param vA needs to derive from BaseProcess or ProcessSequence
+ @param vB needs to derive from BaseProcess or ProcessSequence
**/
- template <typename TProcess1, typename TProcess2, typename TSelect>
- typename std::enable_if_t<is_process_v<typename std::decay_t<TProcess1>> &&
- is_process_v<typename std::decay_t<TProcess2>>,
- SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
- make_select(TProcess1&& vA, TProcess2&& vB, TSelect selector) {
- return SwitchProcessSequence<TProcess1, TProcess2, TSelect>(vA, vB, selector);
+ template <typename TCondition, typename TSequence, typename USequence,
+ typename = std::enable_if_t<is_process_v<typename std::decay_t<TSequence>> &&
+ is_process_v<typename std::decay_t<USequence>>>>
+ SwitchProcessSequence<TCondition, TSequence, USequence> make_select(
+ TCondition&& selector, TSequence&& vA, USequence&& vB) {
+ return SwitchProcessSequence<TCondition, TSequence, USequence>(selector, vA, vB);
}
+ //! @}
+
} // namespace corsika
#include <corsika/detail/framework/process/SwitchProcessSequence.inl>
diff --git a/corsika/framework/units/quantity.hpp b/corsika/framework/units/quantity.hpp
index 3e897cfc80151beebdd2fb700af35fd27ae94e47..7b5bd35877b755fe1d751965975c148d42a8f20e 100644
--- a/corsika/framework/units/quantity.hpp
+++ b/corsika/framework/units/quantity.hpp
@@ -468,6 +468,9 @@ namespace phys {
template <typename D, typename X>
friend detail::Root<D, 2, X> sqrt(quantity<D, X> const& x);
+ template <typename D, typename X>
+ friend detail::Root<D, 3, X> cbrt(quantity<D, X> const& x);
+
// comparison
template <typename D, typename X, typename Y>
@@ -676,7 +679,17 @@ namespace phys {
static_assert(detail::root<D, 2, X>::all_even_multiples,
"root result dimensions must be integral");
- return detail::Root<D, 2, X>(std::pow(x.m_value, X(1.0) / 2));
+ return detail::Root<D, 2, X>(std::sqrt(x.m_value));
+ }
+
+ /// cubic root.
+
+ template <typename D, typename X>
+ detail::Root<D, 3, X> cbrt(quantity<D, X> const& x) {
+ static_assert(detail::root<D, 3, X>::all_even_multiples,
+ "root result dimensions must be integral");
+
+ return detail::Root<D, 3, X>(std::cbrt(x.m_value));
}
// Comparison operators
diff --git a/corsika/framework/units/quantity_io_celsius.hpp b/corsika/framework/units/quantity_io_celsius.hpp
index cc145ecdfe073814c089e8a746c357a8ee189813..354ea74d64c2db4284b7fc6d29df59180b2bc7c5 100644
--- a/corsika/framework/units/quantity_io_celsius.hpp
+++ b/corsika/framework/units/quantity_io_celsius.hpp
@@ -22,14 +22,14 @@
namespace phys { namespace units {
/**
- * celsius, [°C].
+ * celsius, [C].
*/
template<>
struct unit_info< thermodynamic_temperature_d >
{
static bool single() { return true; }
static std::string name() { return "celsius"; }
- static std::string symbol() { return "°C"; }
+ static std::string symbol() { return "C"; }
};
namespace literals {
diff --git a/corsika/framework/units/quantity_io_engineering.hpp b/corsika/framework/units/quantity_io_engineering.hpp
index 16b4be8f59ddc8f9676604063054029fb894913e..b6665370da9c1a8462c3da4af8fd15148dbc565f 100644
--- a/corsika/framework/units/quantity_io_engineering.hpp
+++ b/corsika/framework/units/quantity_io_engineering.hpp
@@ -31,7 +31,7 @@
#include <sstream>
/*
- * Note: micro, µ, may not work everywhere, so you can define a glyph yourself:
+ * Note: micro, may not work everywhere, so you can define a glyph yourself:
*/
#ifndef ENG_FORMAT_MICRO_GLYPH
# define ENG_FORMAT_MICRO_GLYPH "u"
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index befd3025035105a184f33c1c2a872c750364665e..eac8dde53f3eadc461e26c559e8b8ed36463b647 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -19,14 +19,36 @@
namespace corsika {
/**
+ @defgroup Utilities
+
+ Collection of classes and methods to perform recurring tasks.
+ **/
+
+ /**
+ @class COMBoost
+ @ingroup Utilities
+
This utility class handles Lorentz boost between different
referenence frames, using FourVector.
+
+ The class is initialized with projectile and optionally target
+ energy/momentum data. During initialization, a rotation matrix is
+ calculated to represent the projectile movement (and thus the
+ boost) along the z-axis. Also the inverse of this rotation is
+ calculated. The Lorentz boost matrix and its inverse are
+ determined as 2x2 matrices considering the energy and
+ pz-momentum.
+
+ Different constructors are offered with different specialization
+ for the cases of collisions (projectile-target) or just decays
+ (projectile only).
*/
class COMBoost {
public:
- //! construct a COMBoost given four-vector of projectile and mass of target
+ //! construct a COMBoost given four-vector of projectile and mass of target (target at
+ //! rest)
COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
HEPEnergyType const massTarget);
@@ -41,9 +63,11 @@ namespace corsika {
template <typename FourVector>
FourVector fromCoM(FourVector const& p) const;
+ //! returns the rotated coordinate system
CoordinateSystemPtr getRotatedCS() const;
protected:
+ //! internal method
void setBoost(double coshEta, double sinhEta);
private:
diff --git a/corsika/framework/utility/CorsikaData.hpp b/corsika/framework/utility/CorsikaData.hpp
index d117db8da59e3399e77ff7304335c7e9eaa21711..2ab2dec7dabcf7a81bcbe976d689bd92f72df1b3 100644
--- a/corsika/framework/utility/CorsikaData.hpp
+++ b/corsika/framework/utility/CorsikaData.hpp
@@ -8,13 +8,18 @@
#pragma once
-#include <string>
+#include <boost/filesystem/path.hpp>
namespace corsika {
/**
+ * @file CorsikaData.hpp
+ * @ingroup Utilities
+ * @{
* returns the full path of the file \p filename within the CORSIKA_DATA directory
*/
- std::string corsika_data(std::string const& filename);
+ boost::filesystem::path corsika_data(boost::filesystem::path const& filename);
+
+ //! @}
} // namespace corsika
#include <corsika/detail/framework/utility/CorsikaData.inl>
diff --git a/corsika/framework/utility/QuarticSolver.hpp b/corsika/framework/utility/QuarticSolver.hpp
index 0bcdfcedf1a9a9055306acacfd756012726cf8b0..ad9484939d10691e38c7e74d9e2632b3d2e60463 100644
--- a/corsika/framework/utility/QuarticSolver.hpp
+++ b/corsika/framework/utility/QuarticSolver.hpp
@@ -11,6 +11,9 @@
#include <complex>
/**
+ * @file QuarticSolver.hpp
+ * @ingroup Utilities
+ * @{
* \todo convert to class
*/
@@ -56,6 +59,8 @@ namespace corsika::quartic_solver {
// afterwards.
DComplex* solve_quartic(double a, double b, double c, double d);
+ //! @}
+
} // namespace corsika::quartic_solver
#include <corsika/detail/framework/utility/QuarticSolver.inl>
diff --git a/corsika/framework/utility/SaveBoostHistogram.hpp b/corsika/framework/utility/SaveBoostHistogram.hpp
index 14e0379966c447cd0c2f8cfa7ccffcf3c29a0478..16e39180c01668ac0f83694f77fb84ed010dd49d 100644
--- a/corsika/framework/utility/SaveBoostHistogram.hpp
+++ b/corsika/framework/utility/SaveBoostHistogram.hpp
@@ -13,6 +13,9 @@
namespace corsika {
/**
+ * @file SaveBoostHistogram.hpp
+ * @ingroup Utilities
+ *
* This functions saves a boost::histogram into a numpy file. Only rather basic axis
* types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
* counts (i.e. a double or int) are supported, nothing fancy like profiles.
diff --git a/corsika/media/MediumProperties.hpp b/corsika/media/MediumProperties.hpp
index 76434f1f710e3799b5fc74e25ad83334015ac834..8d20fd248bf1b2c1703410124abb3b3a91ad4dbe 100644
--- a/corsika/media/MediumProperties.hpp
+++ b/corsika/media/MediumProperties.hpp
@@ -11,12 +11,56 @@
namespace corsika {
/**
- * Medium types are useful most importantly for effective models
- * like energy losses. a particular medium (mixture of components)
- * may have specif properties not reflected by its mixture of
- * components.
- **/
+ @defgroup MediaProperties Media Properties
+
+ Medium types are useful most importantly for effective models
+ like energy losses. a particular medium (mixture of components)
+ may have specif properties not reflected by its mixture of
+ components.
+
+ The data provided here is automatically parsed from the file
+ properties8.dat from NIST.
+
+ The data of each known medium can be access via the global functions in namespace
+ corsika, or via a static class object with the following interface (here at the
+ example of the class HydrogenGas):
+
+ @code{.cpp}
+ static constexpr Medium medium() { return Medium::HydrogenGas; }
+
+ static std::string const getName() { return data_.getName(); }
+ static std::string const getPrettyName() { return data_.getPrettyName(); }
+ static double getWeight() { return data_.getWeight (); }
+ static int weight_significant_figure() { return data_.weight_significant_figure (); }
+ static int weight_error_last_digit() { return data_.weight_error_last_digit(); }
+ static double Z_over_A() { return data_.Z_over_A(); }
+ static double getSternheimerDensity() { return data_.getSternheimerDensity(); }
+ static double getCorrectedDensity() { return data_.getCorrectedDensity(); }
+ static State getState() { return data_.getState(); }
+ static MediumType getType() { return data_.getType(); }
+ static std::string const getSymbol() { return data_.getSymbol(); }
+
+ static double getIeff() { return data_.getIeff(); }
+ static double getCbar() { return data_.getCbar(); }
+ static double getX0() { return data_.getX0(); }
+ static double getX1() { return data_.getX1(); }
+ static double getAA() { return data_.getAA(); }
+ static double getSK() { return data_.getSK(); }
+ static double getDlt0() { return data_.getDlt0(); }
+
+ inline static const MediumData data_ { "hydrogen_gas", "hydrogen gas (H%2#)", 1.008,
+ 3, 7, 0.99212,
+ 8.3748e-05, 8.3755e-05, State::DiatomicGas,
+ MediumType::Element, "H", 19.2, 9.5835, 1.8639, 3.2718, 0.14092, 5.7273, 0.0 };
+ @endcode
+
+ The numeric data known to CORSIKA 8 (and obtained from NIST) can be browsed below.
+ @ingroup MediaProperties
+ @{
+ */
+
+ //! General type of medium
enum class MediumType {
Unknown,
Element,
@@ -28,15 +72,16 @@ namespace corsika {
BiologicalDosimetry
};
+ //! Physical state of medium
enum class State { Unknown, Solid, Liquid, Gas, DiatomicGas };
enum class Medium : int16_t;
+
using MediumIntType = std::underlying_type<Medium>::type;
- /** \todo documentation needs update ...
- * \struct MediumData
+ /**
*
- * Simple object to group together a number of properties
+ * Simple object to group together the properties of a medium.
*
**/
struct MediumData {
@@ -59,34 +104,59 @@ namespace corsika {
double sk_;
double dlt0_;
- std::string getName() const { return name_; }
- std::string getPrettyName() const { return pretty_name_; }
- double getWeight() const { return weight_; }
- const int& weight_significant_figure() const { return weight_significant_figure_; }
- const int& weight_error_last_digit() const { return weight_error_last_digit_; }
- const double& Z_over_A() const { return Z_over_A_; }
- double getSternheimerDensity() const { return sternheimer_density_; }
- double getCorrectedDensity() const { return corrected_density_; }
- State getState() const { return state_; }
- MediumType getType() const { return type_; }
- std::string getSymbol() const { return symbol_; }
- double getIeff() const { return Ieff_; }
- double getCbar() const { return Cbar_; }
- double getX0() const { return x0_; }
- double getX1() const { return x1_; }
- double getAA() const { return aa_; }
- double getSK() const { return sk_; }
- double getDlt0() const { return dlt0_; }
+ //! @name MediumDataInterface Interface methods
+ //! Interface functions for MediumData
+ //! @{
+ std::string getName() const { return name_; } /// returns name
+ std::string getPrettyName() const { return pretty_name_; } /// returns pretty name
+ double getWeight() const { return weight_; } /// return weight
+ const int& weight_significant_figure() const {
+ return weight_significant_figure_;
+ } /// return significnat figures of weight
+ const int& weight_error_last_digit() const {
+ return weight_error_last_digit_;
+ } /// return error of weight
+ const double& Z_over_A() const { return Z_over_A_; } /// Z_over_A_
+ double getSternheimerDensity() const {
+ return sternheimer_density_;
+ } /// Sternheimer density
+ double getCorrectedDensity() const {
+ return corrected_density_;
+ } /// corrected density
+ State getState() const { return state_; } /// state
+ MediumType getType() const { return type_; } /// type
+ std::string getSymbol() const { return symbol_; } /// symbol
+ double getIeff() const { return Ieff_; } /// Ieff
+ double getCbar() const { return Cbar_; } /// Cbar
+ double getX0() const { return x0_; } /// X0
+ double getX1() const { return x1_; } /// X1
+ double getAA() const { return aa_; } /// AA
+ double getSK() const { return sk_; } /// Sk
+ double getDlt0() const { return dlt0_; } /// Delta0
+ //! @}
};
+ //! @}
+
} // namespace corsika
#include <corsika/media/GeneratedMediaProperties.inc>
namespace corsika {
+ /**
+ @file MediaProperties.hpp
+
+ @ingroup MediaProperties
+ @{
+
+ Returns MediumData object for medium identifed by enum Medium.
+ */
+
constexpr MediumData const& mediumData(Medium const m) {
return corsika::detail::medium_data[static_cast<MediumIntType>(m)];
}
+ //! @}
+
} // namespace corsika
diff --git a/corsika/media/NuclearComposition.hpp b/corsika/media/NuclearComposition.hpp
index 5341554cd3f132d7c6415112a6eca72058de025b..e4ad9513a24684eed2ca1c3b9862441003f2c0fa 100644
--- a/corsika/media/NuclearComposition.hpp
+++ b/corsika/media/NuclearComposition.hpp
@@ -53,9 +53,9 @@ namespace corsika {
**/
size_t getSize() const;
- /// Returns a const reference to the fraction
+ //! Returns a const reference to the fraction
std::vector<float> const& getFractions() const;
- /// Returns a const reference to the fraction
+ //! Returns a const reference to the fraction
std::vector<Code> const& getComponents() const;
double const getAverageMassNumber() const;
@@ -67,6 +67,9 @@ namespace corsika {
// must be updated, too!
size_t getHash() const;
+ //! based on hash value
+ bool operator==(NuclearComposition const& v) const;
+
private:
void updateHash();
diff --git a/corsika/modules/StackInspector.hpp b/corsika/modules/StackInspector.hpp
index 9616079540a21792148f54186bbe5890a1400b43..174e8e05c335af7e4c5f89ed4e621767a13abea6 100644
--- a/corsika/modules/StackInspector.hpp
+++ b/corsika/modules/StackInspector.hpp
@@ -26,7 +26,7 @@ namespace corsika {
StackInspector(const int vNStep, const bool vReportStack, const HEPEnergyType vE0);
~StackInspector();
- void doStack(const TStack&);
+ void doStack(TStack const&);
/**
* To set a new E0, for example when a new shower event is started
diff --git a/corsika/modules/conex/CONEXhybrid.hpp b/corsika/modules/conex/CONEXhybrid.hpp
index b8927f6356bea06dc89845b011517f594a48b928..299f85cbda7ae0298002ae8213ae5c8065950735 100644
--- a/corsika/modules/conex/CONEXhybrid.hpp
+++ b/corsika/modules/conex/CONEXhybrid.hpp
@@ -39,6 +39,9 @@ namespace corsika {
CoordinateSystemPtr const& getObserverCS() const { return conexObservationCS_; }
+ HEPEnergyType getEnergyEM() const;
+ void reset();
+
private:
// data members
//! CONEX e.m. particle codes
@@ -52,6 +55,7 @@ namespace corsika {
CoordinateSystemPtr const conexObservationCS_; //!< CONEX observation frame
DirectionVector const x_sf_,
y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction
+ HEPEnergyType energy_em_;
};
} // namespace corsika
diff --git a/corsika/modules/qgsjetII/qgsjet-II-04.hpp b/corsika/modules/qgsjetII/qgsjet-II-04.hpp
index f95d011ae1c1a7784a40aef7012eb2d05e32a6cf..9d3b884fb71a5f0e15801f994021a325bdbaba49 100644
--- a/corsika/modules/qgsjetII/qgsjet-II-04.hpp
+++ b/corsika/modules/qgsjetII/qgsjet-II-04.hpp
@@ -60,40 +60,32 @@ void qgaini_(
const char* datdir); // Note: there is a length limiation 132 from fortran-qgsjet here
/**
- @function qgini_
-
additional initialization procedure per event
- @parameter e0n - interaction energy (per hadron/nucleon),
- @parameter icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), +-4 - K+-, +-5 -
+ @param e0n - interaction energy (per hadron/nucleon),
+ @param icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), +-4 - K+-, +-5 -
K_l/s),
- @parameter iap - projectile mass number (1 - for a hadron),
- @parameter iat - target mass number
+ @param iap - projectile mass number (1 - for a hadron),
+ @param iat - target mass number
*/
void qgini_(const double& e0n, const int& icp0, const int& iap, const int& iat);
/**
- @function qgconf_
-
generate one event configuration
*/
void qgconf_();
/**
- @function qgsect_
-
hadron-nucleus (hadron-nucleus) particle production cross section
- @parameter e0n lab. energy per projectile nucleon (hadron)
- @parameter icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
- @parameter iap projectile mass number (1=<iap<=iapmax),
- @parameter iat target mass number (1=<iat<=iapmax)
+ @param e0n lab. energy per projectile nucleon (hadron)
+ @param icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
+ @param iap projectile mass number (1=<iap<=iapmax),
+ @param iat target mass number (1=<iat<=iapmax)
*/
double qgsect_(const double& e0n, const int& icz, const int& iap0, const int& iat0);
/**
- @function qgran
-
link to random number generation
*/
double qgran_(int&);
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index 5eb28160818abf42f7a4a743165515f2254d89dc..513ffbaf509d2c05a7020f038d0ab320592c4712 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -30,7 +30,6 @@ namespace corsika {
/**
* \file TrackingLeapFrogCurved.hpp
- * \function LeapFrogStep
*
* Performs one leap-frog step consistent of two halve-steps with steplength/2
* The step is caluculated analytically precisely to reach to the next volume
@@ -73,6 +72,9 @@ namespace corsika {
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
+ static std::string getName() { return "LeapFrog-curved"; }
+ static std::string getVersion() { return "1.0.0"; }
+
protected:
/**
* Use internally stored class tracking_line::Tracking to
diff --git a/corsika/modules/tracking/TrackingLeapFrogStraight.hpp b/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
index f54e355ba3ca568e75caff18be305bf665f73557..2ea6ebd34ee3d9b0eaa48563cae7c8a41aed1df9 100644
--- a/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
@@ -61,6 +61,9 @@ namespace corsika {
template <typename Particle>
auto getTrack(Particle& particle);
+ static std::string getName() { return "LeapFrogStraight"; }
+ static std::string getVersion() { return "1.0.0"; }
+
protected:
double firstFraction_;
};
diff --git a/corsika/modules/tracking/TrackingStraight.hpp b/corsika/modules/tracking/TrackingStraight.hpp
index 4b498bfa522d5856b2251435a437c7ec42681687..a5800cd7ed64fa4fee8afe89a32fe927366ad1de 100644
--- a/corsika/modules/tracking/TrackingStraight.hpp
+++ b/corsika/modules/tracking/TrackingStraight.hpp
@@ -51,6 +51,9 @@ namespace corsika::tracking_line {
//! find intersection of Plane with Track
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
+
+ static std::string getName() { return "Tracking-Straight"; }
+ static std::string getVersion() { return "1.0.0"; }
};
} // namespace corsika::tracking_line
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 52439d7234a287718c06a471859e8e35b31fe8fd..3f6d3e9936f1fd13fcb831359f8faee6768f2948 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -12,21 +12,28 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
#include <corsika/setup/SetupStack.hpp>
+#include <boost/filesystem/path.hpp>
+#include <boost/multi_array.hpp>
+
#include <array>
#include <utility>
+#include <string>
namespace corsika::urqmd {
class UrQMD : public InteractionProcess<UrQMD> {
public:
- UrQMD();
+ UrQMD(boost::filesystem::path const& path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"));
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
+ CrossSectionType getTabulatedCrossSection(Code, Code, HEPEnergyType) const;
+
template <typename TParticle>
CrossSectionType getCrossSection(TParticle const&, Code) const;
@@ -35,13 +42,16 @@ namespace corsika::urqmd {
bool canInteract(Code) const;
+ void blob(int) {}
+
private:
static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int);
+ void readXSFile(boost::filesystem::path const&);
// data members
default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("urqmd");
-
std::uniform_int_distribution<int> booleanDist_{0, 1};
+ boost::multi_array<CrossSectionType, 3> xs_interp_support_table_;
};
/**
diff --git a/corsika/output/BaseOutput.hpp b/corsika/output/BaseOutput.hpp
index 023b9b884ccb0594f4d95f9beefb9eebe6c4a5d5..db998f336c8b445c3a8466fb9f6059bb8c5378ec 100644
--- a/corsika/output/BaseOutput.hpp
+++ b/corsika/output/BaseOutput.hpp
@@ -7,7 +7,7 @@
*/
#pragma once
-#include <boost/filesystem>
+#include <boost/filesystem.hpp>
#include <yaml-cpp/yaml.h>
diff --git a/corsika/output/OutputManager.hpp b/corsika/output/OutputManager.hpp
index aa6c52a8dfaeef500e9e42e3ac6843094e1724b4..2582220520454d028b660ca06f5ac93ed8a18b73 100644
--- a/corsika/output/OutputManager.hpp
+++ b/corsika/output/OutputManager.hpp
@@ -32,7 +32,7 @@ namespace corsika {
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.
+ 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.
diff --git a/corsika/setup/SetupTrajectory.hpp b/corsika/setup/SetupTrajectory.hpp
index c40351d26f1c7bdef9d9b8ac18d3653caaf1927d..95a91115b0aa422c860448fa8057d70c35f95456 100644
--- a/corsika/setup/SetupTrajectory.hpp
+++ b/corsika/setup/SetupTrajectory.hpp
@@ -38,15 +38,15 @@ namespace corsika::setup {
The default tracking algorithm.
*/
- // typedef corsika::tracking_leapfrog_curved::Tracking Tracking;
+ typedef corsika::tracking_leapfrog_curved::Tracking Tracking;
// typedef corsika::tracking_leapfrog_straight::Tracking Tracking;
- typedef corsika::tracking_line::Tracking Tracking;
+ // typedef corsika::tracking_line::Tracking Tracking;
/**
The default trajectory.
*/
/// definition of Trajectory base class, to be used in tracking and cascades
- typedef StraightTrajectory Trajectory;
- // typedef corsika::LeapFrogTrajectory Trajectory;
+ // typedef StraightTrajectory Trajectory;
+ typedef corsika::LeapFrogTrajectory Trajectory;
} // namespace corsika::setup
diff --git a/corsika/stack/WeightStackExtension.hpp b/corsika/stack/WeightStackExtension.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b2e6081f2638c0fe116e5ad87d18af0502e71605
--- /dev/null
+++ b/corsika/stack/WeightStackExtension.hpp
@@ -0,0 +1,102 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <tuple>
+#include <utility>
+#include <vector>
+
+namespace corsika::weights {
+
+ /**
+ * Describe "particle weights" on a Stack.
+ *
+ * Corresponding defintion of a stack-readout object, the iteractor
+ * dereference operator will deliver access to these function
+ * defintion of a stack-readout object, the iteractor dereference
+ * operator will deliver access to these function
+ */
+
+ /**
+ * \tparam TParentStack The stack to be exteneded here with weight data
+ */
+ template <typename TParentStack>
+ struct WeightDataInterface : public TParentStack {
+
+ typedef TParentStack super_type;
+
+ public:
+ // default version for particle-creation from input data
+ void setParticleData(std::tuple<double> const v);
+ void setParticleData(WeightDataInterface const& parent, std::tuple<double> const);
+ void setParticleData();
+ void setParticleData(WeightDataInterface const& parent);
+
+ std::string asString() const;
+
+ void setWeight(double const v);
+
+ double getWeight() const;
+ };
+
+ // definition of stack-data object to store geometry information
+
+ /**
+ * @class WeightData
+ *
+ * definition of stack-data object to store geometry information
+ */
+ class WeightData {
+
+ public:
+ typedef std::vector<double> weight_vector_type;
+
+ WeightData() = default;
+ WeightData(WeightData const&) = default;
+ WeightData(WeightData&&) = default;
+ WeightData& operator=(WeightData const&) = default;
+ WeightData& operator=(WeightData&&) = default;
+
+ // these functions are needed for the Stack interface
+ void clear();
+
+ unsigned int getSize() const;
+
+ unsigned int getCapacity() const;
+
+ void copy(int const i1, int const i2);
+
+ void swap(int const i1, int const i2);
+
+ // custom data access function
+ void setWeight(int const i, double const v);
+
+ double getWeight(int const i) const;
+
+ // these functions are also needed by the Stack interface
+ void incrementSize();
+
+ void decrementSize();
+
+ // custom private data section
+ private:
+ weight_vector_type weight_vector_;
+ };
+
+ template <typename TParentStack>
+ struct MakeWeightDataInterface {
+ typedef WeightDataInterface<TParentStack> type;
+ };
+
+} // namespace corsika::weights
+
+#include <corsika/detail/stack/WeightStackExtension.inl>
diff --git a/documentation/CMakeLists.txt b/documentation/CMakeLists.txt
index e62d3c71db4fb27081f910d7623a53eadf85c7bf..5e727190a8bb493e5b5a4f679226970303b084a1 100644
--- a/documentation/CMakeLists.txt
+++ b/documentation/CMakeLists.txt
@@ -1,5 +1,6 @@
find_package (Doxygen OPTIONAL_COMPONENTS dot mscgen dia)
-
+find_package (Sphinx)
+
if (DOXYGEN_FOUND)
if (NOT DOXYGEN_DOT_EXECUTABLE)
message (FATAL_ERROR "Found doxygen but not 'dot' command, please install graphviz or set DOXYGEN_DOT_EXECUTABLE")
@@ -10,23 +11,43 @@ if (DOXYGEN_FOUND)
set (DOXYGEN_GENERATE_HTML YES)
set (DOXYGEN_GENERATE_MAN YES)
configure_file (${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
- message ("Start doxygen with \"make doxygen\"")
+ message ("Start doxygen with \"make docs\"")
# note the option ALL which allows to build the docs together with the application
- add_custom_target (doxygen # ALL
+ add_custom_target (docs # ALL
COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
COMMENT "Generating API documentation with Doxygen"
VERBATIM
)
- add_custom_command (TARGET doxygen POST_BUILD
+ add_custom_command (TARGET docs POST_BUILD
COMMAND cd ${CMAKE_CURRENT_BINARY_DIR}/latex; pdflatex refman.tex
)
install (DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html DESTINATION share/doc OPTIONAL)
install (FILES ${CMAKE_CURRENT_BINARY_DIR}/latex/refman.pdf DESTINATION share/doc OPTIONAL)
-
+
+ if (SPHINX_FOUND)
+
+ set (SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR})
+ set (SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR}/docs/sphinx)
+
+ add_custom_command (TARGET docs POST_BUILD
+ COMMAND ${SPHINX_EXECUTABLE} -b html
+ # Tell Breathe where to find the Doxygen output
+ -Dbreathe_projects.CORSIKA8=${CMAKE_CURRENT_BINARY_DIR}/xml/
+ ${SPHINX_SOURCE} ${SPHINX_BUILD}
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+ COMMENT "Generating documentation with Sphinx"
+ )
+
+ else (SPHINX_FOUND)
+
+ message ("Sphinx not found;")
+
+ endif (SPHINX_FOUND)
+
else (DOXYGEN_FOUND)
message ("Doxygen need to be installed to generate the doxygen documentation")
diff --git a/documentation/Doxyfile.in b/documentation/Doxyfile.in
index 1378999e9fd140e792d576fdceefce4556bff5e4..ded50d6a5ee5068ce65aa1a15338654f98f3ec72 100644
--- a/documentation/Doxyfile.in
+++ b/documentation/Doxyfile.in
@@ -7,12 +7,12 @@ GENERATE_LATEX = YES
GENERATE_XML = YES
OUTPUT_DIRECTORY = @CMAKE_CURRENT_BINARY_DIR@/
-INPUT = @PROJECT_SOURCE_DIR@/corsika @PROJECT_SOURCE_DIR@/src
+INPUT = @PROJECT_SOURCE_DIR@/corsika @PROJECT_SOURCE_DIR@/src @CMAKE_BINARY_DIR@/corsika
EXCLUDE_PATTERNS = *.inl
FULL_PATH_NAMES = YES
STRIP_FROM_PATH = @PROJECT_SOURCE_DIR@
-FILE_PATTERNS = *.cpp *.hpp *.dox *.md
+FILE_PATTERNS = *.cpp *.hpp *.dox *.md *.inc
EXTENSION_MAPPING = inc=C++
RECURSIVE = YES
@@ -34,13 +34,16 @@ GENERATE_LEGEND = YES
USE_MATHJAX = YES
MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols
-#JAVADOC_BANNER = YES
+#JAVADOC_BANNER = YES
JAVADOC_AUTOBRIEF = YES
-REPEAT_BRIEF = YES
+REPEAT_BRIEF = YES
ABBREVIATE_BRIEF = YES
-MULTILINE_CPP_IS_BRIEF = YES
-MARKDOWN_SUPPORT = YES
-BUILTIN_STL_SUPPORT = YES
-CLANG_ASSISTED_PARSING = YES
+MULTILINE_CPP_IS_BRIEF = YES
+MARKDOWN_SUPPORT = YES
+BUILTIN_STL_SUPPORT = YES
+
+#CLANG_ASSISTED_PARSING = YES
+#CLANG_DATABASE_PATH = @CMAKE_CURRENT_BINARY_DIR@/
+#CLANG_OPTIONS = "-std=c++17"
SEARCHENGINE = YES
diff --git a/documentation/api.rst b/documentation/api.rst
new file mode 100644
index 0000000000000000000000000000000000000000..f8436257b4f5bfb33b642b0a7ac4bd452774634c
--- /dev/null
+++ b/documentation/api.rst
@@ -0,0 +1,11 @@
+Full API
+========
+
+Consider using Doxygen directly for a full API reference....
+
+..
+ .. doxygenindex::
+ :project: CORSIKA8
+ :path: c8
+ :outline:
+ :no-link:
diff --git a/documentation/conf.py b/documentation/conf.py
new file mode 100644
index 0000000000000000000000000000000000000000..ab28020201a0af1b306226ac3afad29d2a606b85
--- /dev/null
+++ b/documentation/conf.py
@@ -0,0 +1,155 @@
+import sys
+import subprocess, os
+
+def configureDoxyfile(input_dir, output_dir):
+ with open('Doxyfile.in', 'r') as file :
+ filedata = file.read()
+
+ filedata = filedata.replace('@PROJECT_SOURCE_DIR@', input_dir)
+ filedata = filedata.replace('@CMAKE_CURRENT_BINARY_DIR@', output_dir)
+ filedata = filedata.replace('@CMAKE_BINARY_DIR@', output_dir)
+
+ with open('Doxyfile', 'w') as file:
+ file.write(filedata)
+
+read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True'
+
+breathe_projects = {}
+
+if read_the_docs_build:
+
+ input_dir = '../'
+ output_dir = 'build'
+ configureDoxyfile(input_dir, output_dir)
+
+ subprocess.call('mkdir -p build/corsika/framework/core; cd build/corsika/framework/core && ../../../../../src/framework/core/pdxml_reader.py ../../../../../src/framework/core/ParticleData.xml ../../../../../src/framework/core/NuclearData.xml ../../../../../src/framework/core/ParticleClassNames.xml', shell=True)
+
+ subprocess.call('mkdir -p build/corsika/media; cd build/corsika/media && ../../../../src/media/readProperties.py ../../../../src/media/properties8.dat', shell=True)
+
+ subprocess.call('doxygen', shell=True)
+ breathe_projects['CORSIKA8'] = output_dir + '/xml'
+
+
+# -- Project information -----------------------------------------------------
+
+project = u'CORSIKA8'
+copyright = u'2021, CORSIKA 8 Collaboration'
+author = u'CORSIKA 8 Collaboration'
+
+# The short X.Y version
+version = u'0.0.0'
+# The full version, including alpha/beta/rc tags
+release = u'prototype-0.1.0'
+
+extensions = [
+ 'sphinx.ext.todo',
+ 'breathe',
+ 'sphinx.ext.mathjax',
+ 'sphinx.ext.autodoc',
+ 'sphinx.ext.extlinks',
+ 'sphinx.ext.autosummary',
+ 'sphinx.ext.napoleon',
+ 'recommonmark'
+]
+
+breathe_default_project = "CORSIKA8"
+
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The master toctree document.
+master_doc = 'index'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set "language" from the command line for these cases.
+language = None
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path .
+exclude_patterns = []
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+
+# -- Options for HTML output -------------------------------------------------
+
+html_theme = 'sphinx_rtd_theme'
+# html_theme_options = {}
+
+# html_static_path = []
+# html_sidebars = {}
+
+htmlhelp_basename = 'C8doc'
+
+# -- Options for LaTeX output ------------------------------------------------
+
+latex_elements = {
+ # The paper size ('letterpaper' or 'a4paper').
+ #
+ # 'papersize': 'letterpaper',
+
+ # The font size ('10pt', '11pt' or '12pt').
+ #
+ # 'pointsize': '10pt',
+
+ # Additional stuff for the LaTeX preamble.
+ #
+ # 'preamble': '',
+
+ # Latex figure (float) alignment
+ #
+ # 'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+# author, documentclass [howto, manual, or own class]).
+latex_documents = [
+ (master_doc, 'C8-GUIDE.tex', u'CORSIKA 8 Documentation',
+ u'CORSIKA 8 Collaboration', 'manual'),
+]
+
+
+# -- Options for manual page output ------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+ (master_doc, 'CORSIKA 8', u'CORSIKA 8 Documentation',
+ [author], 1)
+]
+
+
+# -- Options for Texinfo output ----------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+# dir menu entry, description, category)
+texinfo_documents = [
+ (master_doc, 'CORSIKA 8', u'CORSIKA 8 Documentation',
+ author, 'CORSIKA 8', 'The air shower simulation framework',
+ 'Miscellaneous'),
+]
+
+
+# -- Extension configuration -------------------------------------------------
+
+napoleon_google_docstring = True
+napoleon_include_init_with_doc = True
+napoleon_use_keyword = True
+napoleon_use_param = True
+napoleon_use_rtype = False
+napoleon_use_admonition_for_examples = True
+# autodoc_docstring_signature=False
diff --git a/documentation/environment.rst b/documentation/environment.rst
new file mode 100644
index 0000000000000000000000000000000000000000..9761e3918a8b19d697697820516be3d575d89aca
--- /dev/null
+++ b/documentation/environment.rst
@@ -0,0 +1,4 @@
+Geometry and Environment
+========================
+
+Not yet documented in sphinx. Check doxygen, examples, tests.
diff --git a/documentation/index.rst b/documentation/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..79440c503bd0b500d8235d5c1edd77c2bf459822
--- /dev/null
+++ b/documentation/index.rst
@@ -0,0 +1,19 @@
+CORSIKA 8
+*********
+
+Welcome to the CORSIKA 8 air shower simulation framework.
+
+.. toctree::
+ :maxdepth: 2
+
+ readme_link
+ modules
+ particles
+ media
+ units
+ environment
+ stack
+ utilities
+ api
+
+
diff --git a/documentation/media.rst b/documentation/media.rst
new file mode 100644
index 0000000000000000000000000000000000000000..6a24954c02ee6158c54015bd1be12592c4625d1f
--- /dev/null
+++ b/documentation/media.rst
@@ -0,0 +1,9 @@
+Media Properties
+================
+
+.. toctree::
+ media_classes
+
+.. doxygengroup:: MediaProperties
+ :project: CORSIKA8
+ :members:
diff --git a/documentation/media_classes.rst b/documentation/media_classes.rst
new file mode 100644
index 0000000000000000000000000000000000000000..0cd1607cc38d3d2ea8e5357e6bb6244e713e262f
--- /dev/null
+++ b/documentation/media_classes.rst
@@ -0,0 +1,7 @@
+Media Classes
+=============
+
+.. doxygengroup:: MediaPropertiesClasses
+ :project: CORSIKA8
+ :members:
+
diff --git a/documentation/modules.rst b/documentation/modules.rst
new file mode 100644
index 0000000000000000000000000000000000000000..bcb292faadf034b78507d07eedbebde8071525c5
--- /dev/null
+++ b/documentation/modules.rst
@@ -0,0 +1,9 @@
+Physics modules and processes
+=============================
+
+.. doxygengroup:: Processes
+ :project: CORSIKA8
+ :members:
+
+
+
diff --git a/documentation/particle_classes.rst b/documentation/particle_classes.rst
new file mode 100644
index 0000000000000000000000000000000000000000..63e184dc5017929e78d16bedeaa69476c2af76f3
--- /dev/null
+++ b/documentation/particle_classes.rst
@@ -0,0 +1,7 @@
+Particle Classes
+================
+
+.. doxygengroup:: ParticleClasses
+ :project: CORSIKA8
+ :members:
+
diff --git a/documentation/particles.rst b/documentation/particles.rst
new file mode 100644
index 0000000000000000000000000000000000000000..33f22f996ceebbd735665fc2d8d088e7e0b7dc7b
--- /dev/null
+++ b/documentation/particles.rst
@@ -0,0 +1,12 @@
+Particle Properties
+===================
+
+.. toctree::
+ particle_classes
+
+.. doxygengroup:: Particles
+ :project: CORSIKA8
+ :members:
+
+
+
diff --git a/documentation/readme_link.rst b/documentation/readme_link.rst
new file mode 100644
index 0000000000000000000000000000000000000000..72a33558153fb57def85612b021ec596ef2a51b9
--- /dev/null
+++ b/documentation/readme_link.rst
@@ -0,0 +1 @@
+.. include:: ../README.rst
diff --git a/documentation/rtfd-requirements.txt b/documentation/rtfd-requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3087817421ed52d662b781f58b78ca9e2ae4a38a
--- /dev/null
+++ b/documentation/rtfd-requirements.txt
@@ -0,0 +1,4 @@
+sphinx_rtd_theme == 0.5.1
+sphinx == 3.5.3
+breathe == 4.27.0
+recommonmark == 0.7.1
diff --git a/documentation/stack.rst b/documentation/stack.rst
new file mode 100644
index 0000000000000000000000000000000000000000..3aa4a5bf1835e937cf2ffa00e06ca2006e3c0cad
--- /dev/null
+++ b/documentation/stack.rst
@@ -0,0 +1,4 @@
+Particle storage in memory
+==========================
+
+Not yet documented in sphinx. Check doxygen, examples, tests.
diff --git a/documentation/units.rst b/documentation/units.rst
new file mode 100644
index 0000000000000000000000000000000000000000..af92b833fb27512fc03c67eeeb1c70f15d18d4b3
--- /dev/null
+++ b/documentation/units.rst
@@ -0,0 +1,4 @@
+Physics Units
+=============
+
+Not yet documented in sphinx. Check doxygen, examples, tests.
diff --git a/documentation/utilities.rst b/documentation/utilities.rst
new file mode 100644
index 0000000000000000000000000000000000000000..37c289fde6a68f36548ff2fcdab0ce0dbda2114b
--- /dev/null
+++ b/documentation/utilities.rst
@@ -0,0 +1,9 @@
+Utilities
+==========
+
+.. doxygengroup:: Utilities
+ :project: CORSIKA8
+ :members:
+
+
+
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 0444935039776970b24a3c1e6180001b8c996042..ed4324b3e066d69d0ffd75c34d91da4908b14da7 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -229,15 +229,12 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
- SwitchResult operator()(const setup::Stack::particle_type& p) {
- if (p.getEnergy() < cutE_)
- return SwitchResult::First;
- else
- return SwitchResult::Second;
+ bool operator()(const setup::Stack::particle_type& p) {
+ return (p.getEnergy() < cutE_);
}
};
- auto hadronSequence = make_select(
- urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted), EnergySwitch(55_GeV));
+ auto hadronSequence = make_select(EnergySwitch(55_GeV), urqmdCounted,
+ make_sequence(sibyllNucCounted, sibyllCounted));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
auto sequence = make_sequence(hadronSequence, reset_particle_mass, decaySequence, eLoss,
cut, conex_model, longprof, observationLevel);
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 3b759abfadb2547282e1115010f2b137254fe2be..3ada2b875dbf5b8b3c7627fcf1883e2f91027a1b 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -95,8 +95,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
// argv : 1.number of nucleons, 2.number of protons,
// 3.total energy in GeV, 4.number of showers,
-// 5.output directory
-// 6.seed (0 by default to generate random values for all)
+// 5.seed (0 by default to generate random values for all)
int main(int argc, char** argv) {
@@ -106,7 +105,7 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("vertical_EAS");
if (argc < 5) {
- std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed] \n"
+ std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] \n"
" if A=0, Z is interpreted as PDG code \n"
" if no seed is given, a random seed is chosen \n"
<< std::endl;
@@ -114,15 +113,10 @@ int main(int argc, char** argv) {
}
feenableexcept(FE_INVALID);
- string output_dir = "";
int seed = 0;
int number_showers = std::stoi(std::string(argv[4]));
- if (argc > 5) {
- output_dir = argv[5];
- if (output_dir.back() != '/' && output_dir != "") output_dir += '/';
- }
- if (argc > 6) { seed = std::stoi(std::string(argv[6])); }
+ if (argc > 5) { seed = std::stoi(std::string(argv[5])); }
// initialize random number sequence(s)
registerRandomStreams(seed);
@@ -165,17 +159,131 @@ int main(int argc, char** argv) {
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;
+ 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 = 20.;
+ double phi = 180.;
+ auto const thetaRad = theta / 180. * M_PI;
+ auto const phiRad = phi / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, phiRad, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << ",phi=" << phi << 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 = -observationHeight * cos(thetaRad) +
+ sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+ static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore + DirectionVector{rootCS,
+ {-sin(thetaRad) * cos(phiRad),
+ -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
+ t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+
+ // 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};
+
+ // 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, 60_GeV, 60_GeV, 60_GeV, true};
+ corsika::proposal::Interaction emCascade(env);
+ corsika::proposal::ContinuousProcess emContinuous(env);
+ InteractionCounter emCascadeCounted(emCascade);
+
+ OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+
+ LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+
+ corsika::urqmd::UrQMD urqmd;
+ InteractionCounter urqmdCounted{urqmd};
+ StackInspector<setup::Stack> stackInspect(50000, false, E0);
+
+ // assemble all processes into an ordered process list
+ struct EnergySwitch {
+ HEPEnergyType cutE_;
+ EnergySwitch(HEPEnergyType cutE)
+ : cutE_(cutE) {}
+ bool operator()(const Particle& p) { return (p.getEnergy() < cutE_); }
+ };
+ auto hadronSequence = make_select(EnergySwitch(55_GeV), urqmdCounted,
+ make_sequence(sibyllNucCounted, sibyllCounted));
+ auto decaySequence = make_sequence(decayPythia, decaySibyll);
+
for (int i_shower = 1; i_shower < number_showers + 1; i_shower++) {
// directory for outputs
- string labHist_dir =
- output_dir + "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
- string cMSHist_dir =
- output_dir + "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
- string longprof_dir =
- output_dir + "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
- string tracks_dir = output_dir + "tracks_" + to_string(i_shower) + ".dat";
- string particles_dir = output_dir + "particles_" + to_string(i_shower) + ".dat";
+ 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;
@@ -183,48 +291,6 @@ 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 = -observationHeight * cos(thetaRad) +
- sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
- static_pow<2>(injectionHeight));
- 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));
@@ -351,9 +417,9 @@ int main(int argc, char** argv) {
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
- save_hist(hists.labHist(), labHist_dir, true);
- save_hist(hists.CMSHist(), cMSHist_dir, true);
- longprof.save(longprof_dir);
+ 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/CMakeLists.txt b/modules/qgsjetII/CMakeLists.txt
index bef449e65c6e91cbc2066c73dcb389a1e2e1f87f..c566ea2eb61e91c7e4935e73f2bf18898b208dc7 100644
--- a/modules/qgsjetII/CMakeLists.txt
+++ b/modules/qgsjetII/CMakeLists.txt
@@ -10,20 +10,6 @@ set (
)
enable_language (Fortran)
-add_library (QGSJetII SHARED ${MODEL_SOURCES})
-set_target_properties (
- QGSJetII
- PROPERTIES
- POSITION_INDEPENDENT_CODE 1
- )
-
-target_include_directories (
- QGSJetII
- PUBLIC
- $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
- $<INSTALL_INTERFACE:include/corsika_modules/qgsjetII>
- )
-target_link_libraries (QGSJetII CorsikaData)
add_library (QGSJetII_static STATIC ${MODEL_SOURCES})
set_target_properties (
@@ -46,12 +32,6 @@ install (
DESTINATION include/corsika_modules/qgsjetII
)
-install (
- TARGETS QGSJetII
- EXPORT CORSIKA8PublicTargets
- DESTINATION lib/corsika
- )
-
install (
TARGETS QGSJetII_static
EXPORT CORSIKA8PublicTargets
diff --git a/modules/qgsjetII/qgsjet-II-04.hpp b/modules/qgsjetII/qgsjet-II-04.hpp
index 7b41f01ed1413874a27a36b1f8a0b41f4ba2f5f2..cf58e8dd72a008a385bb7e8a0cae3bf2a9d691b7 100644
--- a/modules/qgsjetII/qgsjet-II-04.hpp
+++ b/modules/qgsjetII/qgsjet-II-04.hpp
@@ -99,8 +99,8 @@ void qgconf_();
@parameter e0n lab. energy per projectile nucleon (hadron)
@parameter icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
- @parameter iap projectile mass number (1=<iap<=iapmax),
- @parameter iat target mass number (1=<iat<=iapmax)
+ @parameter iap0 projectile mass number (1=<iap0<=iapmax),
+ @parameter iat0 target mass number (1=<iat0<=iapmax)
*/
double qgsect_(const double& e0n, const int& icz, const int& iap0, const int& iat0);
diff --git a/modules/sibyll/CMakeLists.txt b/modules/sibyll/CMakeLists.txt
index c974237befbc23ad84283a71759a94867aa5a7e6..7f03881a4462e66016290e218fe52d46021a95e3 100644
--- a/modules/sibyll/CMakeLists.txt
+++ b/modules/sibyll/CMakeLists.txt
@@ -14,28 +14,14 @@ set (
)
enable_language (Fortran)
-add_library (Sibyll SHARED ${MODEL_SOURCES})
add_library (Sibyll_static STATIC ${MODEL_SOURCES})
-set_target_properties (
- Sibyll
- PROPERTIES
- POSITION_INDEPENDENT_CODE 1
- )
-
set_target_properties (
Sibyll_static
PROPERTIES
POSITION_INDEPENDENT_CODE 1
)
-target_include_directories (
- Sibyll
- PUBLIC
- $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
- $<INSTALL_INTERFACE:include/corsika_modules/sibyll>
- )
-
target_include_directories (
Sibyll_static
PUBLIC
@@ -55,12 +41,6 @@ install (
DESTINATION include/corsika_modules/sibyll
)
-install (
- TARGETS Sibyll
- EXPORT CORSIKA8PublicTargets
- LIBRARY DESTINATION lib/corsika
- )
-
install (
TARGETS Sibyll_static
EXPORT CORSIKA8PublicTargets
diff --git a/modules/urqmd/CMakeLists.txt b/modules/urqmd/CMakeLists.txt
index b6048484806a24e2760f0684074def6726de1e7d..77bb20787191dffe145b5ae85648d32f98a397ca 100644
--- a/modules/urqmd/CMakeLists.txt
+++ b/modules/urqmd/CMakeLists.txt
@@ -41,17 +41,6 @@ set (
enable_language (Fortran)
-add_library (UrQMD SHARED ${MODEL_SOURCES})
-target_include_directories (UrQMD PUBLIC
- $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
- $<INSTALL_INTERFACE:include/corsika_modules/urqmd>
- )
-set_target_properties (
- UrQMD
- PROPERTIES
- POSITION_INDEPENDENT_CODE 1
- )
-
add_library (UrQMD_static STATIC ${MODEL_SOURCES})
target_include_directories (UrQMD_static PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
@@ -69,13 +58,6 @@ install (
DESTINATION include/corsika_modules/urqmd
)
-install (
- TARGETS UrQMD
- EXPORT CORSIKA8PublicTargets
- LIBRARY DESTINATION lib/corsika
- INCLUDES DESTINATION include/corsika_modules/urqmd
- )
-
install (
TARGETS UrQMD_static
EXPORT CORSIKA8PublicTargets
diff --git a/modules/urqmd/urqmdInterface.F b/modules/urqmd/urqmdInterface.F
index ec2bd092a8386393368a6789df73b3efa3253160..446221dfc2b2f311b11a4c0080e42ba02fab7296 100644
--- a/modules/urqmd/urqmdInterface.F
+++ b/modules/urqmd/urqmdInterface.F
@@ -432,16 +432,3 @@ c~ xsegymin=0.25d0
end
-c
-c M. Reininghaus, 2020-04-08
-c
- integer function ReadSigmaLn(ia, ib, ic)
- implicit none
-
- include 'comres.f'
-
- integer :: ia, ib, ic
-
- ReadSigmaLn = SigmaLn(ia, ib, ic)
-
- end function ReadSigmaLn
diff --git a/src/framework/core/ParticleData.xml b/src/framework/core/ParticleData.xml
index db89d04628abe4613e0c4a673575bddb550aed9e..1727681045ba9b15c92a0206607ea758725cc60e 100644
--- a/src/framework/core/ParticleData.xml
+++ b/src/framework/core/ParticleData.xml
@@ -138,7 +138,6 @@
m0="0.00000">
</particle>
-<!--
<particle id="23" name="Z0" spinType="3" chargeType="0" colType="0"
m0="91.18760" mWidth="2.49520" mMin="10.00000" mMax="0.00000">
<channel onMode="1" bRatio="0.1539950" products="1 -1"/>
@@ -166,7 +165,7 @@
<channel onMode="1" bRatio="0.1081660" products="-13 14"/>
<channel onMode="1" bRatio="0.1080870" products="-15 16"/>
</particle>
- -->
+
<particle id="25" name="h0" spinType="1" chargeType="0" colType="0"
m0="125.00000" mWidth="0.00374" mMin="50.00000" mMax="0.00000">
<channel onMode="1" bRatio="0.0000009" products="1 -1"/>
diff --git a/src/framework/core/pdxml_reader.py b/src/framework/core/pdxml_reader.py
index 2520cdfe4adb3f75bf316043425ca47e3e3742d7..a099dc427ff835ee6bad75bdd5169513ea100d05 100755
--- a/src/framework/core/pdxml_reader.py
+++ b/src/framework/core/pdxml_reader.py
@@ -275,7 +275,8 @@ def gen_conversion_PDG_ngc(particle_db):
# return string with enum of all internal particle codes
#
def gen_internal_enum(particle_db):
- string = ("enum class Code : CodeIntType {\n"
+ string = ("//! @cond EXCLUDE_DOXY\n"
+ "enum class Code : CodeIntType {\n"
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
@@ -284,7 +285,7 @@ def gen_internal_enum(particle_db):
string += " {key:s} = {code:d},\n".format(key = k, code = last_ngc_id)
string += (" LastParticle = {:d},\n" # identifier for eventual loops...
- "}};").format(last_ngc_id + 1)
+ "}}; //! @endcond").format(last_ngc_id + 1)
if last_ngc_id > 0x7fff: # does not fit into int16_t
raise Exception("Integer overflow in internal particle code definition prevented!")
@@ -297,13 +298,14 @@ def gen_internal_enum(particle_db):
# return string with enum of all PDG particle codes
#
def gen_pdg_enum(particle_db):
- string = "enum class PDGCode : PDGCodeType {\n"
+ string = ("//! @cond EXCLUDE_DOXY\n"
+ "enum class PDGCode : PDGCodeType {\n")
for cId in particle_db:
pdgCode = particle_db[cId]['pdg']
string += " {key:s} = {code:d},\n".format(key = cId, code = pdgCode)
- string += " };\n"
+ string += " }; //! @endcond \n"
return string
@@ -407,7 +409,9 @@ def gen_properties(particle_db):
#
def gen_classes(particle_db):
- string = "// list of C++ classes to access particle properties"
+ string = ("// list of C++ classes to access particle properties\n"
+ "/** @defgroup ParticleClasses \n"
+ " @{ */\n")
for cname in particle_db:
if cname == "Nucleus":
@@ -433,6 +437,7 @@ def gen_classes(particle_db):
string += " * - nuclear Z=" + str(particle_db[cname]['Z']) + "\n"
string += "*/\n\n"
string += "class " + cname + " {\n"
+ string += " /** @cond EXCLUDE_DOXY */ \n"
string += " public:\n"
string += " " + cname + "() = delete;\n"
string += " static constexpr Code code{Code::" + cname + "};\n"
@@ -447,8 +452,11 @@ def gen_classes(particle_db):
string += " static constexpr int nucleus_Z{corsika::get_nucleus_Z(code)};\n"
string += " private:\n"
string += " static constexpr CodeIntType TypeIndex = static_cast<CodeIntType>(code);\n"
+ string += " /** @endcond */ \n"
string += "};\n"
+ string += " //! @}\n";
+
return string
@@ -459,7 +467,11 @@ def inc_start():
string = ('// generated by pdxml_reader.py\n'
'// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. \n'
'\n'
- 'namespace corsika {\n')
+ 'namespace corsika {\n'
+ '/** @ingroup Particles \n'
+ ' @{ \n'
+ ' */ \n'
+ )
return string
@@ -467,7 +479,8 @@ def inc_start():
#
#
def detail_start():
- string = ('namespace particle::detail {\n\n')
+ string = ('/** @cond EXCLUDE_DOXY */\n'
+ ' namespace particle::detail {\n\n')
return string
@@ -475,14 +488,14 @@ def detail_start():
#
#
def detail_end():
- string = "\n}//end namespace particle::detail\n"
+ string = "\n}//end namespace particle::detail\n /** @endcond */"
return string
###############################################################
#
#
def inc_end():
- string = "} // end namespace corsika"
+ string = "/** @} */\n} // end namespace corsika"
return string
diff --git a/src/media/readProperties.py b/src/media/readProperties.py
index 44eb3f839b068d7d26ca457f5d327cabec6b9591..eb89a0c34cc06231d0ead176b0caaca6d37923f6 100755
--- a/src/media/readProperties.py
+++ b/src/media/readProperties.py
@@ -219,6 +219,7 @@ def read_data(filename):
def gen_code(media_db):
string = """
+ /** @cond EXCLUDE_DOXY */
// enum for all media
enum class Medium : MediumIntType {
Unkown,
@@ -231,7 +232,7 @@ def gen_code(media_db):
imedium += 1
string += " {} = {} ,\n".format("First", 0)
string += " {} = {} ,\n".format("Last", imedium-1)
- string += " };\n\n"
+ string += " };/** @endcond */\n\n"
return string
@@ -242,11 +243,11 @@ def gen_code(media_db):
def gen_classes(media_db):
string = """
- // list of C++ classes to access media properties"
-
-// typedef std::map<std::string, double> Properties;
-// typedef std::array<Properties, static_cast<MediumIntType>(Medium::Last)+1> Constituents; this is wrong> num_elements
+ /** @defgroup MediaPropertiesClasses
+ list of C++ classes to access media properties
+ @{
+ */
"""
for entry in media_db:
@@ -278,13 +279,17 @@ def gen_classes(media_db):
class_string = """
/**
- * \class {cname}
+ * @class {cname}
*
* Media properties from properties8.dat file from NIST:
- * - Sternheimer index {stern_index}, label {stern_label}
+ * - Sternheimer index: {stern_index}, label: {stern_label}, name: {name}, nice_name: {nice_name}, symbol: {symbol}
+ * - weight: {weight}, weight_significant_figure: {weight_significant_figure}, weight_error_last_digit: {weight_error_last_digit}
+ * - Z_over_A: {Z_over_A}, sternheimhers_density: {sternheimer_density}, corrected_density: {corrected_density},
+ * - State::{state}, MediumType::{type}, Ieff={Ieff}, Cbar={Cbar}, X0={x0}, x1={x1}, aa={aa}, sk={sk}, dlt0={dlt0}
**/
class {cname} {{
+ /** @cond EXCLUDE_DOXY */
public:
static constexpr Medium medium() {{ return Medium::{cname}; }}
@@ -315,6 +320,7 @@ def gen_classes(media_db):
{weight_significant_figure}, {weight_error_last_digit}, {Z_over_A},
{sternheimer_density}, {corrected_density}, State::{state},
MediumType::{type}, "{symbol}", {Ieff}, {Cbar}, {x0}, {x1}, {aa}, {sk}, {dlt0} }};
+ /** @endcond */
}};
""".format(cname=cname,
@@ -342,32 +348,10 @@ def gen_classes(media_db):
properties=properties);
- # static std::string const name() {{ return "{name}"; }}
- # static std::string const pretty_name() {{ return "{nice_name}"; }}
- # static constexpr double weight() {{ return {weight} ; }}
- # static constexpr int weight_significant_figure() {{ return {weight_significant_figure} ; }}
- # static constexpr int weight_error_last_digit() {{ return {weight_error_last_digit}; }}
- # static constexpr double Z_over_A() {{ return {Z_over_A}; }}
- # static constexpr double sternheimer_density() {{ return {sternheimer_density}; }}
- # static constexpr double corrected_density() {{ return {corrected_density}; }}
- # static constexpr State state() {{ return State::{state}; }}
- # static constexpr MediumType type() {{ return MediumType::{type}; }}
- # static std::string const symbol() {{ return "{symbol}"; }}
-
- # static constexpr double Ieff() {{ return {Ieff}; }}
- # static constexpr double Cbar() {{ return {Cbar}; }}
- # static constexpr double x0() {{ return {x0}; }}
- # static constexpr double x1() {{ return {x1}; }}
- # static constexpr double aa() {{ return {aa}; }}
- # static constexpr double sk() {{ return {sk}; }}
- # static constexpr double dlt0() {{ return {dlt0}; }}
-
-
string += class_string
-
-# private:\n"
-# static constexpr CodeIntType TypeIndex = static_cast<CodeIntType const>(Type);\n"
+ string += " /** @} */"
+
return string
###############################################################
@@ -419,7 +403,9 @@ def inc_start():
#pragma once
namespace corsika {
-
+ /** @ingroup MediaProperties
+ @{
+ */
"""
return string
@@ -428,7 +414,8 @@ namespace corsika {
#
#
def detail_start():
- string = ('namespace detail {\n\n')
+ string = (' /** @} */ \n'
+ 'namespace detail {\n\n')
return string
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 16ed61826867efb14e0636bfe468b68f70a85411..4d703880103150310330657e6fd11c1431acef3c 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -79,6 +79,8 @@ public:
// trajectory: just go ahead forever
particle.getNode()); // next volume node
}
+ static std::string getName() { return "DummyTracking"; }
+ static std::string getVersion() { return "1.0.0"; }
};
class ProcessSplit : public InteractionProcess<ProcessSplit> {
@@ -92,7 +94,7 @@ public:
}
template <typename TView>
- ProcessReturn doInteraction(TView& view) {
+ void doInteraction(TView& view) {
calls_++;
auto vP = view.getProjectile();
const HEPEnergyType E = vP.getEnergy();
@@ -100,7 +102,6 @@ public:
vP.getPosition(), vP.getTime()));
vP.addSecondary(std::make_tuple(vP.getPID(), E / 2, vP.getMomentum(),
vP.getPosition(), vP.getTime()));
- return ProcessReturn::Interacted;
}
int getCalls() const { return calls_; }
diff --git a/tests/framework/testNullModel.cpp b/tests/framework/testNullModel.cpp
index 0589b1782a85d3385f4bd666a598d655f8f93872..5b1f351e13638a1ac2358bd5f0b9d2f295563350 100644
--- a/tests/framework/testNullModel.cpp
+++ b/tests/framework/testNullModel.cpp
@@ -23,5 +23,10 @@ TEST_CASE("NullModel", "[processes]") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
- SECTION("interface") { [[maybe_unused]] NullModel model; }
+ SECTION("interface") {
+ [[maybe_unused]] NullModel nm;
+
+ CHECK(is_process_v<decltype(nm)>);
+ CHECK(count_processes<decltype(nm)>::count == 0);
+ }
}
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index efdaa2fc45fe55945b331d6dd741cfaf294dee09..02105eb9e343f8cf51fc034a3d708b6f88ec2b76 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -5,6 +5,7 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
+#define CORSIKA_UNIT_TESTING
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
@@ -21,11 +22,44 @@
#include <boost/type_index.hpp>
+/*
+ Unit test for testing all Process types and their arrangement in
+ containers ProcessSequence and SwitchProcessSequence
+ */
+
using namespace corsika;
using namespace std;
static int const nData = 10;
+// DummyNode is only needed for BoundaryCrossingProcess
+struct DummyNode {
+ DummyNode(int v)
+ : data_(v) {}
+ int data_ = 0;
+};
+
+// The stack is non-existent for this example
+struct DummyStack {};
+
+// our data object (particle) is a simple arrary of doubles
+struct DummyData {
+ double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ typedef DummyNode node_type; // for BoundaryCrossingProcess
+};
+
+// there is no real trajectory/track
+struct DummyTrajectory {};
+
+// since there is no stack, there is also no view. This is a simplistic dummy object
+// sufficient here.
+struct DummyView {
+ DummyView(DummyData& p)
+ : p_(p) {}
+ DummyData& p_;
+ DummyData& parent() { return p_; }
+};
+
int globalCount = 0; // simple counter
int checkDecay = 0; // use this as a bit field
@@ -49,7 +83,7 @@ public:
void setStep(LengthType const v) { step_ = v; }
template <typename D, typename T>
- ProcessReturn doContinuous(D& d, T&, bool const flag) const {
+ ProcessReturn doContinuous(D& d, T&, bool flag) const {
flag_ = flag;
CORSIKA_LOG_TRACE("ContinuousProcess1::DoContinuous");
checkCont |= 1;
@@ -281,10 +315,7 @@ public:
TimeType getLifetime(Particle&) const {
return 2_s;
}
- template <typename TView>
- void doDecay(TView&) const {
- checkDecay |= 2;
- }
+ void doDecay(DummyView&) const { checkDecay |= 2; }
};
class Stack1 : public StackProcess<Stack1> {
@@ -292,9 +323,8 @@ public:
Stack1(int const n)
: StackProcess(n) {}
template <typename TStack>
- ProcessReturn doStack(TStack&) {
+ void doStack(TStack const&) {
count_++;
- return ProcessReturn::Ok;
}
int getCount() const { return count_; }
@@ -302,21 +332,21 @@ private:
int count_ = 0;
};
-// The stack is non-existent for this example
-struct DummyStack {};
-// our data object (particle) is a simple arrary of doubles
-struct DummyData {
- double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-};
-// there is no real trajectory/track
-struct DummyTrajectory {};
-// since there is no stack, there is also no view. This is a simplistic dummy object
-// sufficient here.
-struct DummyView {
- DummyView(DummyData& p)
- : p_(p) {}
- DummyData& p_;
- DummyData& parent() { return p_; }
+class Boundary1 : public BoundaryCrossingProcess<Boundary1> {
+public:
+ Boundary1(double const v = 1.0)
+ : v_(v) {}
+
+ template <typename Particle>
+ ProcessReturn doBoundaryCrossing(Particle& p, typename Particle::node_type const& from,
+ typename Particle::node_type const& to) {
+
+ for (int i = 0; i < nData; ++i) { p.data_[i] += v_ * (from.data_ - to.data_); }
+ return ProcessReturn::Ok;
+ }
+
+private:
+ double v_ = 0.0;
};
TEST_CASE("ProcessSequence General", "ProcessSequence") {
@@ -330,7 +360,6 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
const Process4 m4(3);
CHECK(is_process_v<Process1>);
- CHECK_FALSE(is_process_v<DummyData>);
CHECK(is_process_v<decltype(m4)>);
CHECK(is_process_v<decltype(Decay1(1))>);
CHECK(is_process_v<decltype(ContinuousProcess3{3, 3_m})>);
@@ -350,19 +379,39 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
auto sequence1 = make_sequence(m1, m2, m3, m4);
CHECK(is_process_v<decltype(sequence1)>);
CHECK(is_process_v<decltype(m2)>);
- CHECK(is_process_sequence_v<decltype(sequence1)>);
- CHECK_FALSE(is_process_sequence_v<decltype(m2)>);
- CHECK_FALSE(is_switch_process_sequence_v<decltype(sequence1)>);
- CHECK_FALSE(is_switch_process_sequence_v<decltype(m2)>);
+ CHECK(decltype(sequence1)::is_process_sequence);
+ CHECK_FALSE(decltype(m2)::is_process_sequence);
+ CHECK_FALSE(decltype(sequence1)::is_switch_process_sequence);
+ CHECK_FALSE(decltype(m2)::is_switch_process_sequence);
- CHECK_FALSE(is_process_sequence_v<decltype(Decay1(7))>);
- CHECK_FALSE(is_switch_process_sequence_v<decltype(Decay1(7))>);
+ CHECK_FALSE(decltype(Decay1(7))::is_process_sequence);
+ CHECK_FALSE(decltype(Decay1(7))::is_switch_process_sequence);
auto sequence2 = make_sequence(m1, m2, m3);
- CHECK(is_process_sequence_v<decltype(sequence2)> == true);
+ CHECK(decltype(sequence2)::is_process_sequence == true);
auto sequence3 = make_sequence(m4);
- CHECK(is_process_sequence_v<decltype(sequence3)> == true);
+ CHECK(decltype(sequence3)::is_process_sequence == true);
+
+ CHECK(std::is_reference_v<decltype(sequence3.getProcess1())>); // Process4&
+ CHECK(!std::is_reference_v<decltype(sequence3.getProcess2())>); // NullModel
+
+ CHECK(std::is_reference_v<decltype(sequence2.getProcess1())>); // Process1&
+ CHECK(!std::is_reference_v<decltype(sequence2.getProcess2())>); // ProcessSequence
+ CHECK(std::is_reference_v<decltype(
+ sequence2.getProcess2().getProcess1())>); // Process2&
+ CHECK(std::is_reference_v<decltype(
+ sequence2.getProcess2().getProcess2())>); // Process3&
+
+ // and now with rvalue initialization
+
+ auto sequence2_rv = make_sequence(Process1(0), m2, Process3(0));
+ CHECK(!std::is_reference_v<decltype(sequence2_rv.getProcess1())>); // Process1
+ CHECK(!std::is_reference_v<decltype(sequence2_rv.getProcess2())>); // ProcessSequence
+ CHECK(std::is_reference_v<decltype(
+ sequence2_rv.getProcess2().getProcess1())>); // Process2&
+ CHECK(!std::is_reference_v<decltype(
+ sequence2_rv.getProcess2().getProcess2())>); // Process3
}
SECTION("interaction length") {
@@ -479,6 +528,9 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
auto sequence1 = make_sequence(s1, s2);
+ std::cout << boost::typeindex::type_id<decltype(sequence1)>().pretty_name()
+ << std::endl;
+
DummyStack stack;
int const nLoop = 20;
@@ -492,11 +544,37 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
auto sequence2 = make_sequence(cp2, m2);
auto sequence3 = make_sequence(cp2, m2, s1);
- CHECK(is_process_sequence_v<decltype(sequence2)> == true);
- CHECK(is_process_sequence_v<decltype(sequence3)> == true);
+ CHECK(decltype(sequence2)::is_process_sequence == true);
+ CHECK(decltype(sequence3)::is_process_sequence == true);
CHECK(contains_stack_process_v<decltype(sequence2)> == false);
CHECK(contains_stack_process_v<decltype(sequence3)> == true);
}
+
+ SECTION("BoundaryCrossingProcess") {
+
+ globalCount = 0;
+ Boundary1 b1;
+
+ auto sequence1 = make_sequence(b1);
+
+ DummyData particle;
+ DummyNode node_from(5);
+ DummyNode node_to(4);
+
+ int const nLoop = 20;
+ for (int i = 0; i < nLoop; ++i) {
+ sequence1.doBoundaryCrossing(particle, node_from, node_to);
+ }
+
+ for (int i = 0; i < nData; i++) {
+ CORSIKA_LOG_DEBUG("data_[{}]={}", i, particle.data_[i]);
+ CHECK(particle.data_[i] == Approx(nLoop).margin(1e-9));
+ }
+
+ CHECK(decltype(sequence1)::is_process_sequence == true);
+ CHECK(contains_stack_process_v<decltype(sequence1)> == false);
+ CHECK(count_processes<decltype(sequence1)>::count == 1);
+ }
}
TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
@@ -510,10 +588,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
*/
struct SwitchSelect {
- SwitchResult operator()(DummyData const& p) const {
- if (p.data_[0] > 0) return SwitchResult::First;
- return SwitchResult::Second;
- }
+ bool operator()(DummyData const& p) const { return (p.data_[0] > 0); }
};
SwitchSelect select1;
@@ -521,39 +596,48 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
auto cp2 = ContinuousProcess2(0, 2_m);
auto cp3 = ContinuousProcess3(0, 3_m);
- auto sequence1 = make_sequence(Process1(0), cp2, Decay1(0));
- auto sequence2 = make_sequence(cp3, Process2(0), Decay2(0));
+ auto sequence1 = make_sequence(Process1(0), cp2, Decay1(0), Boundary1(1.0));
+ auto sequence2 = make_sequence(cp3, Process2(0), Boundary1(-1.0), Decay2(0));
auto sequence3 = make_sequence(cp1, Process3(0),
- SwitchProcessSequence(sequence1, sequence2, select1));
+ SwitchProcessSequence(select1, sequence1, sequence2));
+
+ auto sequence4 =
+ make_sequence(cp1, Boundary1(2.0), Process3(0),
+ SwitchProcessSequence(select1, sequence1, Boundary1(-1.0)));
SECTION("Check construction") {
- auto sequence_alt =
- make_sequence(cp1, Process3(0),
- make_select(make_sequence(Process1(0), cp2, Decay1(0)),
- make_sequence(cp3, Process2(0), Decay2(0)), select1));
+ auto switch_seq = make_select(select1, sequence1, sequence2);
+ CHECK(decltype(switch_seq)::is_process_sequence);
+ CHECK(decltype(switch_seq)::is_switch_process_sequence);
+ CHECK(decltype(SwitchProcessSequence(select1, sequence1,
+ sequence2))::is_switch_process_sequence);
- auto switch_seq = SwitchProcessSequence(sequence1, sequence2, select1);
- CHECK(is_process_sequence_v<decltype(switch_seq)>);
- CHECK(is_switch_process_sequence_v<decltype(switch_seq)>);
- // CHECK(is_switch_process_sequence_v<decltype(&switch_seq)>);
- CHECK(is_switch_process_sequence_v<decltype(
- SwitchProcessSequence(sequence1, sequence2, select1))>);
+ CHECK(decltype(sequence3)::is_process_sequence);
+ CHECK_FALSE(decltype(sequence3)::is_switch_process_sequence);
- CHECK(is_process_sequence_v<decltype(sequence3)>);
- CHECK_FALSE(is_switch_process_sequence_v<decltype(sequence3)>);
+ auto sps1 = SwitchProcessSequence(select1, sequence1, sequence2);
+ CHECK(decltype(sps1)::is_process_sequence);
+ CHECK(decltype(sps1)::is_switch_process_sequence);
- CHECK(is_process_sequence_v<decltype(
- SwitchProcessSequence(sequence1, sequence2, select1))>);
- CHECK(is_switch_process_sequence_v<decltype(
- SwitchProcessSequence(sequence1, sequence2, select1))>);
+ std::cout << boost::typeindex::type_id<decltype(sequence3)>().pretty_name()
+ << std::endl;
+
+ CHECK(decltype(sequence3)::is_process_sequence);
+ auto sps2 = SwitchProcessSequence(select1, sequence1, sequence2);
+ CHECK(decltype(sps2)::is_process_sequence);
+
+ CHECK(std::is_reference_v<decltype(switch_seq.getCondition())>); //
+ CHECK(std::is_reference_v<decltype(switch_seq.getSequence())>); //
+ CHECK(std::is_reference_v<decltype(switch_seq.getAltSequence())>); //
- // check that same process sequence can be build in different ways
- CHECK(typeid(sequence3) == typeid(sequence_alt));
- CHECK(is_process_sequence_v<decltype(sequence3)>);
- CHECK(is_process_sequence_v<decltype(
- SwitchProcessSequence(sequence1, sequence2, select1))>);
+ // check with rvalue init
+ auto switch_seq_rv =
+ make_select(SwitchSelect(), make_sequence(Process1(0)), Process3(0));
+ CHECK(!std::is_reference_v<decltype(switch_seq_rv.getCondition())>);
+ CHECK(!std::is_reference_v<decltype(switch_seq_rv.getSequence())>);
+ CHECK(!std::is_reference_v<decltype(switch_seq_rv.getAltSequence())>);
}
SECTION("Check interfaces") {
@@ -628,6 +712,33 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
CHECK(checkDecay == 0);
CHECK(checkCont == 0);
CHECK(checkSec == 0);
+
+ // check the SwitchProcessSequence where no process is selected in
+ // selected branch (fallthrough)
+
+ checkDecay = 0;
+ checkInteract = 0;
+ checkSec = 0;
+ checkCont = 0;
+ particle.data_[0] = -100; // data positive
+ sequence4.selectInteraction(view, lambda_select);
+ sequence4.doSecondaries(view);
+ sequence4.selectDecay(view, time_select);
+ sequence4.doSecondaries(view);
+ CHECK(checkInteract == 0);
+ CHECK(checkDecay == 0);
+ CHECK(checkCont == 0);
+ CHECK(checkSec == 0);
+
+ // check that large "select" value will correctly ignore the call
+ lambda_select = 1e5 * square(1_cm) / 1_g;
+ time_select = 1e5 / second;
+ checkDecay = 0;
+ checkInteract = 0;
+ sequence3.selectInteraction(view, lambda_select);
+ sequence3.selectDecay(view, time_select);
+ CHECK(checkInteract == 0);
+ CHECK(checkDecay == 0);
}
SECTION("Check ContinuousProcesses in SwitchProcessSequence") {
@@ -718,6 +829,26 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
ContinuousProcessIndex(step4).getIndex(),
boost::typeindex::type_id<decltype(sequence3)>().pretty_name());
}
+
+ SECTION("Check BoundaryCrossingProcess in SwitchProcessSequence") {
+
+ DummyData particle;
+ DummyNode node_from(1);
+ DummyNode node_to(2);
+
+ particle.data_[0] =
+ 100; // data positive, selects particular branch on SwitchProcessSequence
+
+ sequence4.doBoundaryCrossing(particle, node_from, node_to);
+
+ CHECK(particle.data_[0] == 97); // 100 - 2*1 - 1*1
+
+ particle.data_[0] =
+ -100; // data positive, selects particular branch on SwitchProcessSequence
+
+ sequence4.doBoundaryCrossing(particle, node_from, node_to);
+ CHECK(particle.data_[0] == -101); // -100 - 2*1 + 1*1
+ }
}
TEST_CASE("ProcessSequence Indexing", "ProcessSequence") {
@@ -727,40 +858,38 @@ TEST_CASE("ProcessSequence Indexing", "ProcessSequence") {
SECTION("Indexing") {
- int const n0 = count_continuous<Decay2>::count;
- int const n1 = count_continuous<ContinuousProcess3>::count;
- int const n2 = count_continuous<ContinuousProcess2,
- count_continuous<ContinuousProcess3>::count>::count;
+ int const n0 = count_processes<Decay2>::count;
+ int const n1 = count_processes<ContinuousProcess3>::count;
+ int const n2 = count_processes<ContinuousProcess2,
+ count_processes<ContinuousProcess3>::count>::count;
int const n1_b =
- count_continuous<Process2, count_continuous<ContinuousProcess3>::count>::count;
+ count_processes<Process2, count_processes<ContinuousProcess3>::count>::count;
int const n1_c =
- count_continuous<ContinuousProcess3, count_continuous<Process2>::count>::count;
+ count_processes<ContinuousProcess3, count_processes<Process2>::count>::count;
int const n12 =
- count_continuous<ContinuousProcess2,
- count_continuous<ContinuousProcess3, 10>::count>::count;
+ count_processes<ContinuousProcess2,
+ count_processes<ContinuousProcess3, 10>::count>::count;
int const n11_b =
- count_continuous<Process1,
- count_continuous<ContinuousProcess3, 10>::count>::count;
- int const n11_c = count_continuous<ContinuousProcess3,
- count_continuous<Process1, 10>::count>::count;
+ count_processes<Process1, count_processes<ContinuousProcess3, 10>::count>::count;
+ int const n11_c =
+ count_processes<ContinuousProcess3, count_processes<Process1, 10>::count>::count;
- CHECK(n0 == 0);
+ CHECK(n0 == 1);
CHECK(n1 == 1);
- CHECK(n1_b == 1);
- CHECK(n1_c == 1);
+ CHECK(n1_b == 2);
+ CHECK(n1_c == 2);
CHECK(n2 == 2);
- CHECK(n11_b == 11);
- CHECK(n11_c == 11);
+ CHECK(n11_b == 12);
+ CHECK(n11_c == 12);
CHECK(n12 == 12);
- std::cout << count_continuous<ContinuousProcess3>::count << std::endl;
- std::cout << count_continuous<Process3>::count << std::endl;
+ std::cout << count_processes<ContinuousProcess3>::count << std::endl;
+ std::cout << count_processes<Process3>::count << std::endl;
struct SwitchSelect {
- SwitchResult operator()(DummyData const& p) const {
+ bool operator()(DummyData const& p) const {
std::cout << "SwitchSelect data=" << p.data_[0] << std::endl;
- if (p.data_[0] > 0) return SwitchResult::First;
- return SwitchResult::Second;
+ return (p.data_[0] > 0);
}
};
@@ -769,21 +898,21 @@ TEST_CASE("ProcessSequence Indexing", "ProcessSequence") {
ContinuousProcess1(0, 1_m));
SwitchSelect select1;
- auto switch_seq = SwitchProcessSequence(sequence1, sequence2, select1);
+ auto switch_seq = SwitchProcessSequence(select1, sequence1, sequence2);
auto sequence3 = make_sequence(ContinuousProcess1(0, 1_m), Process3(0), switch_seq);
auto sequence4 = make_sequence(ContinuousProcess1(0, 1_m), Process3(0),
- SwitchProcessSequence(sequence1, sequence2, select1));
+ SwitchProcessSequence(select1, sequence1, sequence2));
- int const switch_seq_n = count_continuous<decltype(switch_seq)>::count;
- int const sequence3_n = count_continuous<decltype(sequence3)>::count;
+ int const switch_seq_n = count_processes<decltype(switch_seq)>::count;
+ int const sequence3_n = count_processes<decltype(sequence3)>::count;
CHECK(decltype(sequence1)::getNumberOfProcesses() == 3);
- CHECK(count_continuous<decltype(sequence1)>::count == 3);
- CHECK(count_continuous<decltype(sequence2)>::count == 4);
+ CHECK(count_processes<decltype(sequence1)>::count == 3);
+ CHECK(count_processes<decltype(sequence2)>::count == 4);
CHECK(switch_seq_n == 7);
CHECK(sequence3_n == 9);
- CHECK(count_continuous<decltype(sequence4)>::count == 9);
+ CHECK(count_processes<decltype(sequence4)>::count == 9);
std::cout << "switch_seq "
<< boost::typeindex::type_id<decltype(switch_seq)>().pretty_name()
diff --git a/tests/framework/testUnits.cpp b/tests/framework/testUnits.cpp
index 67db321d8734fc488a833348dca02c120a831b52..8d15544aca0d5c2dc1923af5853e6de7cafeb4f4 100644
--- a/tests/framework/testUnits.cpp
+++ b/tests/framework/testUnits.cpp
@@ -92,6 +92,9 @@ TEST_CASE("PhysicalUnits", "[Units]") {
const auto E3 = E2 + 100_GeV + pow(10, lgE) * 1_GeV;
CHECK(E3 == 180_GeV);
+
+ CHECK(sqrt(5_GeV * 5_GeV) / 5_GeV == Approx(1));
+ CHECK(cbrt(static_pow<3>(5_GeV)) / 5_GeV == Approx(1));
}
SECTION("Output") {
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index ca16d9f589d0e84932e004b4df440335be90fd9f..3eab2cd6b9e5cdce3efc87ae48a9281145fbcb72 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -46,6 +46,9 @@ TEST_CASE("HomogeneousMedium") {
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
std::vector<float>{1.f});
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
+
+ CHECK_THROWS(NuclearComposition({Code::Proton}, {1.1}));
+ CHECK_THROWS(NuclearComposition({Code::Proton}, {0.99}));
}
TEST_CASE("FlatExponential") {
@@ -216,6 +219,8 @@ TEST_CASE("InhomogeneousMedium") {
inhMedium.getIntegratedGrammage(trajectory, l));
CHECK(rho.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
inhMedium.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
+ CHECK(inhMedium.getNuclearComposition() == composition);
+ CHECK(inhMedium.getMassDensity({gCS, {0_m, 0_m, 0_m}}) == 1_kg / static_pow<3>(1_m));
}
}
@@ -234,6 +239,8 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
builder.addLinearLayer(2_km, 20_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 30_km);
+ CHECK_THROWS(builder.addLinearLayer(0.5_km, 5_km));
+
CHECK(builder.getSize() == 3);
auto const builtEnv = builder.assemble();
diff --git a/tests/media/testShowerAxis.cpp b/tests/media/testShowerAxis.cpp
index 8ba454c342eb3ba1f972c4c6129ab0876fb1bcf9..19d7635aedebd4ff901396af86cc4651ba94b3e2 100644
--- a/tests/media/testShowerAxis.cpp
+++ b/tests/media/testShowerAxis.cpp
@@ -62,8 +62,8 @@ TEST_CASE("Homogeneous Density") {
Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos), *env,
- false, // -> do not throw exceptions
- 20}; // -> number of bins
+ true, // -> do not throw exceptions
+ 20}; // -> number of bins
CHECK(showerAxis.getSteplength() == 500_m);
@@ -80,4 +80,7 @@ TEST_CASE("Homogeneous Density") {
const Vector<dimensionless_d> dir{cs, {0, 0, -1}};
CHECK(showerAxis.getDirection().getComponents(cs) == dir.getComponents(cs));
CHECK(showerAxis.getStart().getCoordinates() == injectionPos.getCoordinates());
+
+ CHECK_THROWS(showerAxis.getX(-1_m));
+ CHECK_THROWS(showerAxis.getX((injectionPos - showerCore).getNorm() + 1_m));
}
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 8b039f87ca2979ebdef9831b287c73ef8a2c0478..4a604fd8bbf32f440db2eb0e073aed684b412fe4 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -18,7 +18,7 @@
#include <string>
#include <cstdlib>
-#include <experimental/filesystem>
+#include <boost/filesystem.hpp>
/*
NOTE, WARNING, ATTENTION
@@ -56,13 +56,12 @@ TEST_CASE("CORSIKA_DATA", "[processes]") {
const char* data = std::getenv("CORSIKA_DATA");
// these CHECKS are needed:
CHECK(data != 0);
- CHECK(std::experimental::filesystem::is_directory(
- std::experimental::filesystem::path(std::string(data) + "/QGSJetII")));
+ CHECK(boost::filesystem::is_directory(boost::filesystem::path(data) / "QGSJetII"));
CORSIKA_LOG_INFO(
"data: {}"
" isDir: {}"
"/QGSJetII",
- data, std::experimental::filesystem::is_directory(std::string(data)));
+ data, boost::filesystem::is_directory(data));
}
}
@@ -80,6 +79,8 @@ TEST_CASE("QgsjetII", "[processes]") {
SECTION("QgsjetII -> Corsika") {
CHECK(Code::PiPlus == corsika::qgsjetII::convertFromQgsjetII(
corsika::qgsjetII::QgsjetIICode::PiPlus));
+ CHECK_THROWS(
+ corsika::qgsjetII::convertFromQgsjetII(corsika::qgsjetII::QgsjetIICode::Unknown));
}
SECTION("Corsika -> QgsjetII") {
@@ -170,4 +171,94 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
Approx(0).margin(1e-2));
}
+
+ SECTION("InteractionInterface Nuclei") {
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Nucleus, 60, 30, 20100_GeV,
+ (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ auto particle = stackPtr->first();
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
+
+ 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
+ int countFragments = 0;
+ for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); }
+ CHECK(countFragments == Approx(2).margin(1)); // this is not physics validation
+ [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
+
+ CHECK(length / (1_g / square(1_cm)) ==
+ Approx(12).margin(2)); // this is not physics validation
+ }
+
+ SECTION("Heavy nuclei") {
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Nucleus, 1000, 1000, 1100_GeV,
+ (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ auto particle = stackPtr->first();
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
+
+ corsika::qgsjetII::Interaction model;
+
+ CHECK_THROWS(
+ model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 10., 1000.));
+ CHECK_THROWS(
+ model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 1000., 10.));
+ CHECK_THROWS(model.doInteraction(view));
+ CHECK_THROWS(model.getInteractionLength(particle));
+ }
+
+ SECTION("Allowed Particles") {
+ { // electron
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ auto particle = stackPtr->first();
+ corsika::qgsjetII::Interaction model;
+ GrammageType const length = model.getInteractionLength(particle);
+ CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
+ }
+ { // pi0 is internally converted into pi+/pi-
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(18).margin(2)); // this is not physics validation
+ }
+ { // rho0 is internally converted into pi-/pi+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(7).margin(2)); // 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());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(25).margin(3)); // 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,
+ (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ }
+ }
}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 28ed39a159c9ceb8bb0960d8adac735fc67a95a2..0f9f4a86243ec0f7f5b0c16e989b9459a5e0a293 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -77,24 +77,41 @@ TEST_CASE("UrQMD") {
auto const& cs = *csPtr;
{ [[maybe_unused]] auto const& env_dummy = env; }
- Code validProjectileCodes[] = {Code::PiPlus, Code::PiMinus, Code::Proton,
- Code::Neutron, Code::KPlus, Code::KMinus,
- Code::K0, Code::K0Bar, Code::K0Long};
+ Code validProjectileCodes[] = {Code::PiPlus, Code::PiMinus, Code::Proton,
+ Code::AntiProton, Code::AntiNeutron, Code::Neutron,
+ Code::KPlus, Code::KMinus, Code::K0,
+ Code::K0Bar, Code::K0Long};
for (auto code : validProjectileCodes) {
- auto [stack, view] = setup::testing::setup_stack(
- code, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ auto [stack, view] = setup::testing::setup_stack(code, 0, 0, 100_GeV, nodePtr, cs);
CHECK(stack->getEntries() == 1);
CHECK(view->getEntries() == 0);
// simple check whether the cross-section is non-vanishing
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Proton) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Nitrogen) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Oxygen) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Argon) / 1_mb > 0);
+ // only nuclei with available tabluated data so far
+ CHECK(urqmd.getInteractionLength(stack->getNextParticle()) > 1_g / square(1_cm));
}
}
+ SECTION("targets options") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
+ auto const& cs = *csPtr;
+ { [[maybe_unused]] auto const& env_dummy = env; }
+ auto [stack, view] =
+ setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs);
+ CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) ==
+ Approx(105).margin(5));
+ }
+
+ SECTION("invalid targets options") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Omega);
+ auto const& cs = *csPtr;
+ { [[maybe_unused]] auto const& env_dummy = env; }
+ auto [stack, view] =
+ setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs);
+ CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
+ }
+
SECTION("nucleus projectile") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
@@ -146,6 +163,45 @@ TEST_CASE("UrQMD") {
Approx(0).margin(1e-2));
}
+ SECTION("\"special\" projectile and target") {
+ {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ [[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::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target
+ }
+
+ {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+ [[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::PiPlus, 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();
+
+ urqmd.doInteraction(*secViewPtr);
+
+ CHECK(sumCharge(*secViewPtr) ==
+ get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+ Approx(0).margin(1e-2));
+ }
+ }
+
SECTION("K0Long projectile") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
@@ -171,4 +227,22 @@ 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/stack/CMakeLists.txt b/tests/stack/CMakeLists.txt
index 8fe1070f199e6cf0ca529d267b6189dc6a990cc6..17305b6a4dbbf330c8a28a6bb11e6a771a52858a 100644
--- a/tests/stack/CMakeLists.txt
+++ b/tests/stack/CMakeLists.txt
@@ -3,6 +3,7 @@ set (test_stack_sources
testHistoryStack.cpp
testHistoryView.cpp
testGeometryNodeStackExtension.cpp
+ testWeightStackExtension.cpp
testDummyStack.cpp
testVectorStack.cpp
testNuclearStackExtension.cpp
diff --git a/tests/stack/testWeightStackExtension.cpp b/tests/stack/testWeightStackExtension.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2754c367913e336325491945f6570e4bd014899b
--- /dev/null
+++ b/tests/stack/testWeightStackExtension.cpp
@@ -0,0 +1,85 @@
+/*
+ * (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.
+ */
+
+#include <corsika/framework/stack/CombinedStack.hpp>
+#include <corsika/stack/DummyStack.hpp>
+#include <corsika/stack/WeightStackExtension.hpp>
+
+using namespace corsika;
+
+#include <catch2/catch.hpp>
+
+#include <iostream>
+using namespace std;
+
+// the Weight stack:
+template <typename TStackIter>
+using DummyWeightDataInterface =
+ typename weights::MakeWeightDataInterface<TStackIter>::type;
+
+// combine dummy stack with geometry information for tracking
+template <typename TStackIter>
+using StackWithGeometryInterface =
+ CombinedParticleInterface<dummy_stack::DummyStack::pi_type, DummyWeightDataInterface,
+ TStackIter>;
+
+using TestStack =
+ CombinedStack<typename dummy_stack::DummyStack::stack_implementation_type,
+ weights::WeightData, StackWithGeometryInterface>;
+
+TEST_CASE("WeightStackExtension", "[stack]") {
+
+ logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+ dummy_stack::NoData noData;
+
+ SECTION("write weights") {
+
+ double const weight = 5.1;
+
+ TestStack s;
+ s.addParticle(std::make_tuple(noData), std::tuple<double>{weight});
+
+ CHECK(s.getEntries() == 1);
+ }
+
+ SECTION("write/read weights") {
+ double const weight = 15;
+
+ TestStack s;
+ auto p = s.addParticle(std::make_tuple(noData));
+ p.setWeight(weight);
+ CHECK(s.getEntries() == 1);
+
+ const auto pout = s.getNextParticle();
+ CHECK(pout.getWeight() == weight);
+ }
+
+ SECTION("stack fill and cleanup") {
+
+ double const weight = 16;
+
+ TestStack s;
+ // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
+ for (int i = 0; i < 99; ++i) {
+ auto p = s.addParticle(std::tuple<dummy_stack::NoData>{noData});
+ p.setWeight(weight);
+ }
+
+ CHECK(s.getEntries() == 99);
+ double v = 0;
+ for (int i = 0; i < 99; ++i) {
+ auto p = s.getNextParticle();
+ v += p.getWeight();
+ p.erase();
+ }
+ CHECK(v == Approx(99 * weight));
+ CHECK(s.getEntries() == 0);
+ }
+}