IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Showing
with 596 additions and 650 deletions
This diff is collapsed.
......@@ -3,9 +3,8 @@
#
# See file AUTHORS for a list of contributors.
#
# This software is distributed under the terms of the GNU General Public
# Licence version 3 (GPL Version 3). See file LICENSE for a full version of
# the license.
# This software is distributed under the terms of the 3-clause BSD license.
# See file LICENSE for a full version of the license.
#
set (CORSIKA8_VERSION @c8_version@)
......@@ -31,10 +30,16 @@ set (COMPILE_OPTIONS @COMPILE_OPTIONS@)
set (CMAKE_VERBOSE_MAKEFILE @CMAKE_VERBOSE_MAKEFILE@)
#+++++++++++++++++++++++++++++
# Setup external dependencies via conan
# external dependencies
# same list as top-level CML.txt, except for Catch2 (not needed here)
#
include (${CMAKE_CURRENT_LIST_DIR}/conanbuildinfo.cmake)
conan_basic_setup (TARGETS)
find_package(Boost COMPONENTS filesystem REQUIRED)
find_package(CLI11 REQUIRED)
find_package(Eigen3 REQUIRED)
find_package(spdlog REQUIRED)
find_package(yaml-cpp REQUIRED)
find_package(Arrow REQUIRED)
find_package(PROPOSAL REQUIRED)
#+++++++++++++++++++++++++++++
# Import Pythia8
......@@ -51,6 +56,35 @@ set_target_properties (
set (Pythia8_FOUND @Pythia8_FOUND@)
message (STATUS "Pythia8 at: @Pythia8_PREFIX@")
#+++++++++++++++++++++++++++++
# Import TAUOLA
# since we always import TAUOLA (ExternalProject_Add) we have to
# import here, too.
#
add_library (C8::ext::tauola::CxxInterface STATIC IMPORTED GLOBAL)
add_library (C8::ext::tauola::Fortran STATIC IMPORTED GLOBAL)
add_library(C8::ext::tauola INTERFACE IMPORTED)
set_property(TARGET C8::ext::tauola
PROPERTY
INTERFACE_LINK_LIBRARIES
C8::ext::tauola::CxxInterface
C8::ext::tauola::Fortran)
set_target_properties (
C8::ext::tauola::CxxInterface PROPERTIES
IMPORTED_LOCATION @TAUOLA_LIBDIR@/libTauolaCxxInterface.a
IMPORTED_LINK_INTERFACE_LIBRARIES dl
INTERFACE_INCLUDE_DIRECTORIES @TAUOLA_INCDIR@
)
set_target_properties (
C8::ext::tauola::Fortran PROPERTIES
IMPORTED_LOCATION @TAUOLA_LIBDIR@/libTauolaFortran.a
IMPORTED_LINK_INTERFACE_LIBRARIES dl
INTERFACE_INCLUDE_DIRECTORIES @TAUOLA_INCDIR@
)
set (TAUOLA_FOUND @TAUOLA_FOUND@)
message (STATUS "TAUOLA at: @TAUOLA_PREFIX@")
#++++++++++++++++++++++++++++++
# import CORSIKA8
......
......@@ -3,9 +3,8 @@
#
# See file AUTHORS for a list of contributors.
#
# This software is distributed under the terms of the GNU General Public
# Licence version 3 (GPL Version 3). See file LICENSE for a full version of
# the license.
# This software is distributed under the terms of the 3-clause BSD license.
# See file LICENSE for a full version of the license.
#
......@@ -24,7 +23,7 @@ set (CMAKE_CXX_FLAGS_COVERAGE "-g --coverage")
set (CMAKE_EXE_LINKER_FLAGS_COVERAGE "--coverage")
set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage")
# set a flag to inform code that we are in debug mode
set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG")
set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -D_C8_DEBUG_")
#+++++++++++++++++++++++++++++
......
......@@ -3,9 +3,8 @@
#
# See file AUTHORS for a list of contributors.
#
# This software is distributed under the terms of the GNU General Public
# Licence version 3 (GPL Version 3). See file LICENSE for a full version of
# the license.
# This software is distributed under the terms of the 3-clause BSD license.
# See file LICENSE for a full version of the license.
#
#################################################
......@@ -55,9 +54,10 @@ function (CORSIKA_ADD_TEST)
else ()
set(sanitize ${C8_ADD_TEST_SANITIZE})
endif ()
find_package(Catch2 REQUIRED)
add_executable (${name} ${sources})
target_link_libraries (${name} CORSIKA8 CONAN_PKG::catch2 CorsikaTestingCommon)
target_link_libraries (${name} CORSIKA8 Catch2::Catch2WithMain CorsikaTestingCommon)
target_compile_options (${name} PRIVATE -g) # do not skip asserts
target_include_directories (${name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
file (MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/test_outputs/)
......
#!/bin/bash
function show_usage(){
printf "\n\nUsage:"
printf "$0 options [parameters]\n"
printf "\n"
printf "Options:\n"
printf "\n -s or --source-directory:\n Corsika 8 download directory, which contains the 'conanfile.py' recipe. Default is the current directory."
printf "\n -d or --debug:\n Specify 'Debug' as build type for the installed dependences. This should be matched when building CORSIKA 8."
printf "\n -r or --release:\n Specify 'Release' as build type for the installed dependences. This should be matched when building CORSIKA 8."
printf "\n -rd or --release-with-debug:\n Specify 'RelWithDebInfo' as build type for the installed dependences. This should be matched when building CORSIKA 8."
printf "\n\nExample: ./conan-install.sh --source-directory /some_path/corsika --debug"
printf "\n -h or --help:\n Prints this message.\n"
exit 0
}
echo "|---------------------------------------------------|"
echo "|-----------------[ CORSIKA 8 ]---------------------|"
echo "|-----[ CONAN2 DEPENDENCIES INSTALL SCRIPT ]------- |"
echo "|---------------------------------------------------|"
echo "|-------------------- BEGIN ------------------------| "
echo " "
SOURCE=${BASH_SOURCE[0]}
while [ -L "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink
DIR=$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )
SOURCE=$(readlink "$SOURCE")
[[ $SOURCE != /* ]] && SOURCE=$DIR/$SOURCE # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located
done
SCRIPT_DIR=$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )
BUILD_TYPE="RelWithDebInfo"
CORSIKA_DIR=${CURRENT_DIR}
CORSIKA_DIR_INFORMED=${CURRENT_DIR}
CONAN2_OUTPUT_FOLDER_NAME="conan_cmake"
printf "[ conan-install | info > This script is located at the directory: ${SCRIPT_DIR}\n"
if [[ "$1" == "--help" ]] || [[ "$1" == "-h" ]];then
show_usage
fi
while [ ! -z "$1" ]; do
case "$1" in
--source-directory|-s)
shift
CORSIKA_DIR=$(readlink -e ${1})
CORSIKA_DIR_INFORMED=${1}
;;
--debug|-d)
#shift
BUILD_TYPE="Debug"
;;
--release|-r)
#shift
BUILD_TYPE="Release"
;;
--release-with-debug|-rd)
#shift
BUILD_TYPE="RelWithDebInfo"
;;
--help|-h)
#shift
show_usage
;;
*)
show_usage
;;
esac
if [ $# -gt 0 ]; then
shift
fi
done
if [[ ${#CORSIKA_DIR} -eq 0 ]]; then
printf "[ conan-install | warning > Output folder for cmake scripts generated by conan2 not found.\n"
printf "[ conan-install | warning > Directory '${CORSIKA_DIR_INFORMED}' does not exist.\n"
exit 1
fi
if [ -d "${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}" ]; then
printf "[ conan-install | info > Output folder for cmake scripts generated by conan2 is: ${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}\n"
else
printf "[ conan-install | warning > Output folder for cmake scripts generated by conan2 not found.\n"
printf "[ conan-install | warning > Creating directory ${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}.\n"
mkdir -p "${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}"
fi
if [ -f "${CORSIKA_DIR}/conanfile.py" ]; then
printf "[ conan-install | info > conan2 recipe: ${CORSIKA_DIR}/conanfile.py\n"
else
printf "[ conan-install | warning > conan2 recipe not found.\n"
printf "[ conan-install | warning > File '${CORSIKA_DIR}/conanfile.py' does not exist.\n"
exit 1
fi
# Conan2 variabes
CONAN2_HOME=$(readlink -e `conan config home`)
CONAN2_PROFILE_NAME="corsika8"
printf "[ conan-install | info > conan2 home: ${CONAN2_HOME}\n"
# Conan2 commands
CONAN2_DEFAULT_PROFILE_COMMAND="conan profile detect --force"
CONAN2_PROFILE_COMMAND="conan profile detect --name ${CONAN2_PROFILE_NAME} --force"
CONAN2_INSTALL_COMMAND="conan install ${CORSIKA_DIR} --output-folder=${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME} --build=missing --settings=build_type=${BUILD_TYPE} --profile=${CONAN2_PROFILE_NAME}"
CONAN2_SHOW_PROFLE_COMMAND="conan profile show -pr ${CONAN2_PROFILE_NAME}"
printf "[ conan-install | info > Creating default profile...\n\n"
eval $CONAN2_DEFAULT_PROFILE_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_DEFAULT_PROFILE_COMMAND}.\n"
exit 126
fi
printf "[ conan-install | info > Creating '${CONAN2_PROFILE_NAME}' profile...\n\n"
eval $CONAN2_PROFILE_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_PROFILE_COMMAND}.\n"
exit 126
fi
printf "\n\n[ conan-install | info > Editing '${CONAN2_PROFILE_NAME}' profile\n"
#========== cppstd setting ============
STD_NUMBER=`grep -n -m 1 "compiler.cppstd=" ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME} | cut -d: -f1`
SED_STD_COMMAND="sed -i '${STD_NUMBER}s/.*/compiler.cppstd=gnu17/' ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME}"
eval $SED_STD_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${SED_STD_COMMAND}\n"
exit 126
fi
#========== libcxx setting ============
CXX_NUMBER=`grep -n -m 1 "compiler.libcxx=" ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME} | cut -d: -f1`
SED_CXX_COMMAND="sed -i '${CXX_NUMBER}s/.*/compiler.libcxx=libstdc++11/' ${CONAN2_HOME}/profiles/${CONAN2_PROFILE_NAME}"
eval $SED_CXX_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${SED_CXX_COMMAND}\n"
exit 126
fi
eval $CONAN2_SHOW_PROFLE_COMMAND
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_SHOW_PROFLE_COMMAND}\n"
exit 126
fi
printf "[ conan-install | info > ${CONAN2_INSTALL_COMMAND}\n"
eval ${CONAN2_INSTALL_COMMAND}
if [ ! $? -eq 0 ]; then
printf "[ conan-install | error > Exit code 126 (Command invoked cannot execute):\n ${CONAN2_INSTALL_COMMAND}\n"
exit 126
fi
CORSIKA_CMAKE_SCRIPT="#!/bin/bash
function show_usage(){
printf \"\n\nUsage:\"
printf \"\$0 options [parameters]\n\"
printf \"\n\"
printf \"Options:\n\"
printf \"\n -c or --cmake-flags:\n Additional flags and settings to cmake base command. Default is empty string.\"
printf \"\n\nExample: ./corsika-cmake.sh --cmake-flags '-DUSE_Pythia8_C8=C8' \"
printf \"\n\nNote: the source directory (the one containing CMakeLists.txt), CMAKE_BUILD_TYPE, CMAKE_POLICY_DEFAULT_CMP0091 and CMAKE_TOOLCHAIN_FILE are already set. Do not repeat them.\"
printf \"\n -h or --help:\n Prints this message.\n\"
exit 0
}
echo \"|---------------------------------------------------|\"
echo \"|-----------------[ CORSIKA 8 ]---------------------|\"
echo \"|----------[ CMAKE CONFIGURATION SCRIPT ]---------- |\"
echo \"|---------------------------------------------------|\"
echo \"|-------------------- BEGIN ------------------------|\"
CMAKE_BASE_SETTINGS=\"cmake -S ${CORSIKA_DIR} -D CONAN_CMAKE_DIR=${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME} -D CMAKE_TOOLCHAIN_FILE=${CORSIKA_DIR}/${CONAN2_OUTPUT_FOLDER_NAME}/conan_toolchain.cmake -D CMAKE_POLICY_DEFAULT_CMP0091=NEW -DCMAKE_BUILD_TYPE=${BUILD_TYPE}\"
CMAKE_ADDITIONAL_SETTINGS=\"\"
if [[ \"\$1\" == \"--help\" ]] || [[ \"\$1\" == \"-h\" ]];then
show_usage
fi
while [ ! -z \"\$1\" ]; do
case \"\$1\" in
--cmake-flags|-c)
shift
CMAKE_ADDITIONAL_SETTINGS=\${1}
;;
--help|-h)
#shift
show_usage
;;
*)
show_usage
;;
esac
if [ \$# -gt 0 ]; then
shift
fi
done
printf \"[corsika-cmake | info > Issuing CMake command :\n\"
printf '[corsika-cmake | info > \e[1;36m%s\e[0m\n\n\n' \"\${CMAKE_BASE_SETTINGS} \${CMAKE_ADDITIONAL_SETTINGS}\"
eval \"\${CMAKE_BASE_SETTINGS} \${CMAKE_ADDITIONAL_SETTINGS}\"
if [ ! \$? -eq 0 ]; then
printf \"[ corsika | error > Project configuration (CMake) failed. \n\"
exit 1
fi
echo \"|---------------------------------------------------|\"
echo \"|-----------------[ CORSIKA 8 ]---------------------|\"
echo \"|----------[ CMAKE CONFIGURATION SCRIPT ]---------- |\"
echo \"|---------------------------------------------------|\"
echo \"|-------------------- END --------------------------|\"
"
printf "%s" "${CORSIKA_CMAKE_SCRIPT}" > ${SCRIPT_DIR}/corsika-cmake.sh
chmod +x ${SCRIPT_DIR}/corsika-cmake.sh
printf "\n\n[ conan-install | info > Copy and paste the commands below in the corsika-build directory:\n\n"
printf '\e[1;36m%s\e[0m\n' "> ${SCRIPT_DIR}/corsika-cmake.sh "
printf '\e[1;36m%s\e[0m\n' "> make -j8"
echo " "
echo "|---------------------------------------------------|"
echo "|-----------------[ CORSIKA 8 ]---------------------|"
echo "|-----[ CONAN2 DEPENDENCIES INSTALL SCRIPT ]------- |"
echo "|---------------------------------------------------|"
echo "|-------------------- END --------------------------| "
echo " "
exit 0
from conan import ConanFile
from conan.tools.cmake import cmake_layout,CMakeToolchain, CMakeDeps
class Pkg(ConanFile):
generators = "CMakeDeps", #"CMakeToolchain",
settings = "os", "arch", "compiler", "build_type"
default_options = {
'readline*:shared': 'True',
'arrow*:shared': 'False',
'arrow*:parquet': 'True',
'arrow*:fPIC': 'False',
'arrow*:with_re2': 'True',
'arrow*:with_protobuf': 'False',
'arrow*:with_openssl': 'False',
'arrow*:with_gflags': 'False',
'arrow*:with_glog': 'False',
'arrow*:with_grpc': 'False',
'arrow*:with_utf8proc': 'False',
'arrow*:with_zstd': 'False',
'arrow*:with_bz2': 'False',
'arrow*:with_lz4': 'True',
'arrow*:with_thrift': 'True',
'arrow*:with_boost': 'True',
'boost*:without_container': 'True',
'boost*:without_context': 'True',
'boost*:without_contract': 'True',
'boost*:without_coroutine': 'True',
'boost*:without_date_time': 'True',
'boost*:without_fiber': 'True',
'boost*:without_filesystem': 'False',
'boost*:without_graph': 'True',
'boost*:without_graph_parallel': 'True',
'boost*:without_iostreams': 'False',
'boost*:without_json': 'True',
'boost*:without_locale': 'True',
'boost*:without_log': 'True',
'boost*:without_math': 'False',
'boost*:without_mpi': 'True',
'boost*:without_nowide': 'True',
'boost*:without_program_options': 'True',
'boost*:without_python': 'True',
'boost*:without_serialization': 'False',
'boost*:without_stacktrace': 'True',
'boost*:without_system': 'False',
'boost*:without_test': 'True',
'boost*:without_thread': 'True',
'boost*:without_timer': 'True',
'boost*:without_type_erasure': 'True',
'boost*:without_wave': 'True'
}
def configure(self):
self.options['arrow'].with_boost = True
self.options['arrow'].parquet = True
self.options['arrow'].with_thrift = True
def requirements(self):
self.requires("spdlog/1.14.1", force=True)
self.requires("catch2/3.6.0")
self.requires("bzip2/1.0.8")
self.requires("boost/1.85.0", force=True)
self.requires("eigen/3.4.0")
self.requires("zlib/1.3.1")
self.requires("yaml-cpp/0.8.0")
self.requires("cli11/1.9.1")
self.requires("arrow/16.1.0")
self.requires("proposal/7.6.2")
def build_requirements(self):
self.tool_requires("readline/8.0")
self.tool_requires("bison/[>1.0]")
def generate(self):
tc = CMakeToolchain(self)
tc.absolute_paths = True
tc.generate()
[requires]
spdlog/1.9.2
catch2/2.13.8
eigen/3.3.8
boost/1.76.0
zlib/1.2.11
proposal/7.3.1
yaml-cpp/0.6.3
arrow/2.0.0
cli11/1.9.1
[build_requires]
readline/8.0 # this is only needed to fix a missing dependency in "bison" which is pulled-in by arrow
bison/[>1.0] # needed for arrow, and we HAVE to compile it
[generators]
cmake
[options]
readline:shared=True
arrow:shared=False
arrow:parquet=True
arrow:fPIC=False
arrow:with_re2=False
arrow:with_protobuf=False
arrow:with_openssl=False
arrow:with_gflags=False
arrow:with_glog=False
arrow:with_grpc=False
arrow:with_utf8proc=False
arrow:with_zstd=False
arrow:with_bz2=False
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
// Another possibility:
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Step.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
......@@ -42,10 +42,15 @@ namespace corsika {
, sequence_(pl)
, output_(out)
, stack_(stack)
, forceInteraction_(false) {
, forceInteraction_(false)
, forceDecay_(false) {
CORSIKA_LOG_INFO(c8_ascii_);
CORSIKA_LOG_INFO("This is CORSIKA {}.{}.{}.{}", CORSIKA_RELEASE_NUMBER,
CORSIKA_MAJOR_NUMBER, CORSIKA_MINOR_NUMBER, CORSIKA_PATCH_NUMBER);
CORSIKA_LOG_INFO(
"The C8 author list can be found at: "
"https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/"
"Current-CORSIKA-8-author-list");
CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(),
TTracking::getVersion());
if constexpr (stack_view_type::has_event) {
......@@ -92,6 +97,21 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() {
forceInteraction_ = true;
if (forceDecay_) {
CORSIKA_LOG_ERROR("Cannot set forceInteraction when forceDecay is already set");
throw std::runtime_error(
"Cannot set forceInteraction when forceDecay is already set");
}
}
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceDecay() {
forceDecay_ = true;
if (forceInteraction_) {
CORSIKA_LOG_ERROR("Cannot set forceDecay when forceInteraction is already set");
throw std::runtime_error(
"Cannot set forceDecay when forceInteraction is already set");
}
}
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
......@@ -115,21 +135,37 @@ namespace corsika {
FourMomentum const projectileP4{Elab, particle.getMomentum()};
// determine combined full inelastic cross section of the particles in the material
auto const targetMomentum = MomentumVector{
particle.getMomentum().getCoordinateSystem(), {0_GeV, 0_GeV, 0_GeV}};
auto const xs_function = [&](Code const targetId) -> CrossSectionType {
FourMomentum const targetP4{get_mass(targetId), targetMomentum};
return sequence_.getCrossSection(particle, targetId, targetP4);
};
CrossSectionType const total_cx =
composition.getWeightedSum([=](Code const targetId) -> CrossSectionType {
FourMomentum const targetP4(
get_mass(targetId),
MomentumVector(particle.getMomentum().getCoordinateSystem(),
{0_GeV, 0_GeV, 0_GeV}));
return sequence_.getCrossSection(particle, targetId, targetP4);
});
CrossSectionType const total_cx_pre = composition.getWeightedSum(xs_function);
if (forceInteraction_) {
CORSIKA_LOG_TRACE("forced interaction!");
forceInteraction_ = false; // just one (first) interaction
stack_view_type secondaries(particle);
interaction(secondaries, projectileP4, composition, total_cx);
interaction(secondaries, projectileP4, composition, total_cx_pre);
sequence_.doSecondaries(secondaries);
particle.erase(); // primary particle is done
return;
}
if (forceDecay_) {
CORSIKA_LOG_TRACE("forced decay!");
forceDecay_ = false; // just one decay
stack_view_type secondaries(particle);
decay(secondaries, sequence_.getInverseLifetime(particle));
if (secondaries.getSize() == 1 && secondaries.getProjectile().getPID() ==
secondaries.getNextParticle().getPID()) {
throw std::runtime_error(
fmt::format("Particle {} decays into itself!",
get_name(secondaries.getProjectile().getPID())));
}
sequence_.doSecondaries(secondaries);
particle.erase(); // primary particle is done
return;
......@@ -137,10 +173,10 @@ namespace corsika {
// calculate interaction length in medium
GrammageType const total_lambda =
(composition.getAverageMassNumber() * constants::u) / total_cx;
(composition.getAverageMassNumber() * constants::u) / total_cx_pre;
// sample random exponential step length in grammage
ExponentialDistribution expDist(total_lambda);
ExponentialDistribution expDist{total_lambda};
GrammageType const next_interact = expDist(rng_);
CORSIKA_LOG_DEBUG("total_lambda={} g/cm2, next_interact={} g/cm2",
......@@ -148,31 +184,31 @@ namespace corsika {
double(next_interact / 1_g * 1_cm * 1_cm));
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime = sequence_.getInverseLifetime(particle);
InverseTimeType const total_inv_lifetime_pre = sequence_.getInverseLifetime(particle);
// sample random exponential decay time
ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
ExponentialDistribution expDistDecay(1 / total_inv_lifetime_pre);
TimeType const next_decay = expDistDecay(rng_);
CORSIKA_LOG_DEBUG("total_lifetime={} ns, next_decay={} ns",
(1 / total_inv_lifetime) / 1_ns, next_decay / 1_ns);
(1 / total_inv_lifetime_pre) / 1_ns, next_decay / 1_ns);
// convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * particle.getMomentum().getNorm() /
particle.getEnergy() * constants::c;
// determine geometric tracking
auto [step, nextVol] = tracking_.getTrack(particle);
auto geomMaxLength = step.getLength(1);
auto [track, nextVol] = tracking_.getTrack(particle);
auto geomMaxLength = track.getLength(1);
// convert next_step from grammage to length
LengthType const distance_interact =
currentLogicalNode->getModelProperties().getArclengthFromGrammage(step,
currentLogicalNode->getModelProperties().getArclengthFromGrammage(track,
next_interact);
// determine the maximum geometric step length
ContinuousProcessStepLength const continuousMaxStep =
sequence_.getMaxStepLength(particle, step);
sequence_.getMaxStepLength(particle, track);
LengthType const continuous_max_dist = continuousMaxStep;
// take minimum of geometry, interaction, decay for next step
......@@ -180,16 +216,14 @@ namespace corsika {
LengthType const min_non_continuous = std::min(min_discrete, geomMaxLength);
LengthType const min_distance = std::min(min_non_continuous, continuous_max_dist);
bool const isContinuous = continuous_max_dist < min_non_continuous;
// inform ContinuousProcesses (if applicable) that it is responsible for step-limit
// this would become simpler if we follow the idea of Max to enumerate ALL types of
// processes. Then non-continuous are included and no further logic is needed to
// distinguish between continuous and non-continuous limit.
ContinuousProcessIndex limitingId;
bool const isContinuous = continuous_max_dist < min_non_continuous;
if (isContinuous) {
limitingId =
continuousMaxStep; // the current step IS limited by a known continuous process
}
auto const limitingId = isContinuous ? continuousMaxStep : ContinuousProcessIndex{};
// // the current step IS limited by a known continuous process
CORSIKA_LOG_DEBUG(
"transport particle by : {} m "
......@@ -202,11 +236,12 @@ namespace corsika {
// move particle along the trajectory to new position
// also update momentum/direction/time
step.setLength(min_distance);
track.setLength(min_distance);
Step step{particle, track};
// apply all continuous processes on particle + track
if (sequence_.doContinuous(particle, step, limitingId) ==
ProcessReturn::ParticleAbsorbed) {
if (sequence_.doContinuous(step, limitingId) == ProcessReturn::ParticleAbsorbed) {
CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
particle.getPID(), particle.getEnergy() / 1_GeV);
if (particle.isErased()) {
......@@ -218,9 +253,10 @@ namespace corsika {
}
return; // particle is gone -> return
}
particle.setTime(particle.getTime() + step.getDuration());
particle.setPosition(step.getPosition(1));
particle.setDirection(step.getDirection(1));
particle.setTime(step.getTimePost());
particle.setPosition(step.getPositionPost());
particle.setDirection(step.getDirectionPost());
particle.setKineticEnergy(step.getEkinPost());
if (isContinuous) {
return; // there is nothing further, step is finished
......@@ -267,7 +303,7 @@ namespace corsika {
"Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
if (numericalNodeAfterStep != currentLogicalNode) {
CORSIKA_LOG_ERROR(
CORSIKA_LOG_DEBUG(
"expect to be in node currentLogicalNode={} but are in "
"numericalNodeAfterStep={}. Continue, but without guarantee.",
fmt::ptr(currentLogicalNode), fmt::ptr(numericalNodeAfterStep));
......@@ -284,7 +320,7 @@ namespace corsika {
// secondaries, b) the projectile particle deleted (or
// changed)
stack_view_type secondaries(particle);
stack_view_type secondaries{particle};
/*
Create SecondaryView object on Stack. The data container
......@@ -295,12 +331,19 @@ namespace corsika {
important to use projectile/view (and not particle) for Interaction,
and Decay!
*/
FourMomentum const projectileP4Post{particle.getEnergy(), particle.getMomentum()};
bool eraseParticle =
false; // only erase original particle if it decayed or interacted
if (distance_interact < distance_decay) {
interaction(secondaries, projectileP4, composition, total_cx);
eraseParticle = isInteracted(
interaction(secondaries, projectileP4Post, composition, total_cx_pre));
} else {
[[maybe_unused]] auto projectile = secondaries.getProjectile();
if (decay(secondaries, total_inv_lifetime) == ProcessReturn::Decayed) {
if (decay(secondaries, total_inv_lifetime_pre) == ProcessReturn::Decayed) {
eraseParticle = true;
if (secondaries.getSize() == 1 &&
projectile.getPID() == secondaries.getNextParticle().getPID()) {
throw std::runtime_error(fmt::format("Particle {} decays into itself!",
......@@ -309,8 +352,11 @@ namespace corsika {
}
}
sequence_.doSecondaries(secondaries);
particle.erase();
if (eraseParticle) {
// doSecondaries() makes sense only if there was an actual event
sequence_.doSecondaries(secondaries);
particle.erase();
}
}
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
......@@ -327,7 +373,8 @@ namespace corsika {
auto const returnCode = sequence_.selectDecay(view, sample_process);
if (returnCode != ProcessReturn::Decayed) {
CORSIKA_LOG_DEBUG("Particle did not decay!");
CORSIKA_LOG_ERROR("Particle {} did not decay!",
get_name(view.getProjectile().getPID()));
}
setEventType(view, history::EventType::Decay);
return returnCode;
......
/* -*-c++-*-
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......@@ -14,14 +13,13 @@
namespace corsika {
inline HEPEnergyType constexpr get_kinetic_energy_propagation_threshold(
Code const code) {
inline HEPEnergyType get_kinetic_energy_propagation_threshold(Code const code) {
if (is_nucleus(code)) return particle::detail::threshold_nuclei;
return particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)];
}
inline void constexpr set_kinetic_energy_propagation_threshold(
Code const code, HEPEnergyType const val) {
inline void set_kinetic_energy_propagation_threshold(Code const code,
HEPEnergyType const val) {
if (is_nucleus(code))
particle::detail::threshold_nuclei = val;
else
......@@ -33,12 +31,11 @@ namespace corsika {
return particle::detail::masses[static_cast<CodeIntType>(code)];
}
inline HEPEnergyType constexpr get_energy_production_threshold(Code const p) {
inline HEPEnergyType get_energy_production_threshold(Code const p) {
return particle::detail::production_thresholds[static_cast<CodeIntType>(p)];
}
inline void constexpr set_energy_production_threshold(Code const p,
HEPEnergyType const val) {
inline void set_energy_production_threshold(Code const p, HEPEnergyType const val) {
particle::detail::production_thresholds[static_cast<CodeIntType>(p)] = val;
}
......@@ -134,6 +131,11 @@ namespace corsika {
if (std::abs(k) <= maxPDG) {
return particle::detail::conversionArray[k + maxPDG];
} else {
if (1000000000 <= k && k <= 1009999990) { // nucleus (no L or I)
int const Z = (k - 1000000000) / 10000;
int const A = (k - 1000000000 - 10000 * Z) / 10;
return get_nucleus_code(A, Z);
}
return particle::detail::conversionMap.at(p);
}
}
......@@ -149,7 +151,7 @@ namespace corsika {
return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron);
}
inline std::string_view get_nucleus_name(Code const code) {
inline std::string get_nucleus_name(Code const code) {
size_t const A = get_nucleus_A(code);
size_t const Z = get_nucleus_Z(code);
return fmt::format("Nucleus_A{}_Z{}", A, Z);
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <fmt/format.h>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <boost/filesystem/path.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <Eigen/Dense>
//-----------------------------------
// STD
//-----------------------------------
namespace std {
auto inline format_as(std::_Put_time<char> const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
} // namespace std
//-----------------------------------
// CORSIKA
//-----------------------------------
namespace corsika {
// formatters for particle codes declared on ParticleProperties.hpp
auto inline format_as(Code code) { return get_name(code); }
auto inline format_as(PDGCode code) { return fmt::underlying(code); }
template <typename Type>
auto inline format_as(Type const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
} // namespace corsika
//-----------------------------------
// boost::filesystem
//-----------------------------------
namespace boost::filesystem {
auto inline format_as(path const& fname) { return fname.string().c_str(); }
} // namespace boost::filesystem
//-----------------------------------
// phys::units
//-----------------------------------
namespace phys::units {
template <typename Dimensions>
auto inline format_as(phys::units::quantity<Dimensions> const& arg) {
return io::to_string(arg);
}
} // namespace phys::units
//----------------------------------
// Eigen
//----------------------------------
namespace Eigen {
template <typename Scalar, int M, int N>
auto inline format_as(Matrix<Scalar, M, N> const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
template <typename Matrix1, typename Matrix2>
auto inline format_as(Product<Matrix1, Matrix2> const& arg) {
std::ostringstream os;
os << arg;
return os.str();
}
} // namespace Eigen
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......@@ -81,7 +80,13 @@ namespace corsika {
template <typename TDim>
inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
Vector<TDim> const& vVec) {
auto const a = vVec.normalized().getComponents(cs).getEigenVector();
auto const vVecComp = vVec.getComponents(cs);
if (vVecComp.getX().magnitude() == 0 && vVecComp.getY().magnitude() == 0 &&
vVecComp.getZ().magnitude() == 0) {
return cs;
}
auto const a = vVecComp.normalized().getEigenVector();
auto const a1 = a(0), a2 = a(1), a3 = a(2);
Eigen::Matrix3d A, B;
......
/*
* (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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......
/*
* (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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......@@ -29,8 +28,9 @@ namespace corsika {
}
inline Point LeapFrogTrajectory::getPosition(double const u) const {
Point position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
VelocityVector velocity =
if (u == 0) return initialPosition_;
Point const position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
VelocityVector const velocity =
initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
return position + velocity * timeStep_ * u / 2;
}
......@@ -40,6 +40,8 @@ namespace corsika {
}
inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
if (u == 0) return initialDirection_;
return (initialDirection_ +
initialDirection_.cross(magneticfield_) * timeStep_ * u * k_)
.normalized();
......
/*
* (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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......