diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 53f2e7b4ce54b88a9e45df329eb9a7348829330c..2ef93a896e2a7343746c270b20ea6e32834087c8 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -249,7 +249,6 @@ test-clang-8:
- if: $CI_MERGE_REQUEST_ID
when: manual
allow_failure: true
- allow_failure: true
artifacts:
when: always
expire_in: 3 days
@@ -315,7 +314,6 @@ build_test-clang-8:
- if: $CI_COMMIT_TAG
when: manual
allow_failure: true
- allow_failure: true
artifacts:
when: always
expire_in: 3 days
@@ -373,7 +371,6 @@ build_test_example-clang-8:
- if: $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
- if: $CI_COMMIT_BRANCH
when: manual
- allow_failure: true
cache:
paths:
- ${CI_PROJECT_DIR}/build/
diff --git a/.gitmodules b/.gitmodules
index 5429f2d75f591d52f001615fddde188e8bd012de..9c4be8bab0ad06af04dab930d151388062890e83 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,13 +1,8 @@
[submodule "modules/data"]
- path = modules/data
- url = ../../AirShowerPhysics/corsika-data.git
- branch = master
- shallow = true
-[submodule "modules/proposal"]
- path = modules/proposal
- url = https://github.com/tudo-astroparticlephysics/PROPOSAL.git
- branch = restructure_parametrization
- shallow = true
+ path = modules/data
+ url = ../../AirShowerPhysics/corsika-data.git
+ branch = master
+ shallow = true
[submodule "modules/conex"]
path = modules/conex
url = ../../AirShowerPhysics/cxroot.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d81de38f705fa7f6b5901a5643e74cb2f5f9ab93..7310df1c8f9af36fe354f22a90bdce58efb49d2a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -24,7 +24,7 @@ project (
#+++++++++++++++++++++++++++++
# prevent in-source builds and give warning message
#
-if ("${CMAKE_BINARY_DIR}" STREQUAL "${CMAKE_SOURCE_DIR}")
+if ("${CMAKE_BINARY_DIR}" STREQUAL "${CMAKE_SOURCE_DIR}")
message (FATAL_ERROR "In-source builds are disabled.
Please create a build-dir and use `cmake <source-dir>` inside it.
NOTE: cmake will now create CMakeCache.txt and CMakeFiles/*.
@@ -124,7 +124,7 @@ if (CMAKE_BUILD_TYPE STREQUAL Coverage)
# compile coverage under -O0 to avoid any optimization, function elimation etc.
add_compile_options ("-O0")
-
+
set (GCOV gcov CACHE STRING "gcov executable" FORCE)
set (LCOV_BIN_DIR "${PROJECT_SOURCE_DIR}/externals/lcov/bin")
# collect coverage data
@@ -158,7 +158,7 @@ endif ()
# CTest config
#
enable_testing ()
-if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
+if (${CMAKE_VERSION} VERSION_LESS "3.12.0")
list (APPEND CMAKE_CTEST_ARGUMENTS "--output-on-failure")
else (${CMAKE_VERSION} VERSION_LESS "3.12.0")
set (CTEST_OUTPUT_ON_FAILURE 1) # this is for new versions of cmake
@@ -196,18 +196,19 @@ set_target_properties (
INTERFACE_CORSIKA8_MAJOR_VERSION 0
COMPATIBLE_INTERFACE_STRING CORSIKA8_MAJOR_VERSION
)
-#
+#
target_include_directories (
CORSIKA8
INTERFACE
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>
$<INSTALL_INTERFACE:include>
- )
+ )
# since CORSIKA8 is a header only library we must specify all link dependencies here:
target_link_libraries (
CORSIKA8
INTERFACE
+ CONAN_PKG::proposal
CONAN_PKG::eigen
CONAN_PKG::spdlog
CONAN_PKG::boost
@@ -234,17 +235,10 @@ add_subdirectory (modules/sibyll)
add_subdirectory (modules/qgsjetII)
add_subdirectory (modules/urqmd)
add_subdirectory (modules/conex)
-#
-# we really need a better proposal CMakeList.txt file:
-set (ADD_PYTHON OFF)
-file (READ modules/CMakeLists_PROPOSAL.txt overwrite_CMakeLists_PROPOSAL_txt)
-file (WRITE modules/proposal/CMakeLists.txt ${overwrite_CMakeLists_PROPOSAL_txt})
-add_subdirectory (modules/proposal)
-target_link_libraries (CORSIKA8 INTERFACE PROPOSAL)
#+++++++++++++++++++++++++++++++
# unit testing
-#
+#
add_subdirectory (tests)
#+++++++++++++++++++++++++++++++
@@ -298,8 +292,8 @@ configure_package_config_file (
# config for install tree
# overwrite with install location (if it was build as part of corsika)
-set (Pythia8_INCDIR ${Pythia8_INCDIR_INSTALL})
-set (Pythia8_LIBDIR ${Pythia8_LIBDIR_INSTALL})
+set (Pythia8_INCDIR ${Pythia8_INCDIR_INSTALL})
+set (Pythia8_LIBDIR ${Pythia8_LIBDIR_INSTALL})
configure_package_config_file (
cmake/corsikaConfig.cmake.in
${PROJECT_BINARY_DIR}/cmake/corsikaConfig.cmake
diff --git a/conanfile.txt b/conanfile.txt
index 78cda273ada5bf3c944cb45f6f22df53d0567a8c..645158b384218cbb653f08842c32ea917c9cb1fc 100644
--- a/conanfile.txt
+++ b/conanfile.txt
@@ -5,7 +5,7 @@ catch2/2.13.3
eigen/3.3.8
boost/1.76.0
zlib/1.2.11
-proposal/7.0.2
+proposal/7.0.4
yaml-cpp/0.6.3
arrow/2.0.0
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 430d779c6062549b4e44d5142aca15c8dd85886b..91308390ba4f9b5b8e87bb0d94e79646e9a6f101 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -38,14 +38,19 @@ namespace corsika::proposal {
get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
+ // Use higland multiple scattering and deactivate stochastic deflection by
+ // passing an empty vector
+ static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::Highland;
+ auto s_type = std::vector<PROPOSAL::InteractionType>();
+
// Build displacement integral and scattering object and interpolate them too and
// saved in the calc map by a key build out of a hash of composed of the component and
// particle code.
- auto disp = PROPOSAL::make_displacement(c, true);
- auto scatter =
- PROPOSAL::make_scattering("highland", particle[code], media.at(comp.getHash()));
- calc[std::make_pair(comp.getHash(), code)] =
- std::make_tuple(std::move(disp), std::move(scatter));
+ auto calculator =
+ Calculator{PROPOSAL::make_displacement(c, true),
+ PROPOSAL::make_scattering(ms_type, s_type, particle[code],
+ media.at(comp.getHash()))};
+ calc[std::make_pair(comp.getHash(), code)] = std::move(calculator);
}
template <>
@@ -63,28 +68,25 @@ namespace corsika::proposal {
// Cast corsika vector to proposal vector
auto vP_dir = vP.getDirection();
auto d = vP_dir.getComponents();
- auto direction = PROPOSAL::Vector3D(d.getX().magnitude(), d.getY().magnitude(),
- d.getZ().magnitude());
+ auto direction = PROPOSAL::Cartesian3D(d.getX().magnitude(), d.getY().magnitude(),
+ d.getZ().magnitude());
auto E_f = vP.getEnergy() - loss;
// draw random numbers required for scattering process
std::uniform_real_distribution<double> distr(0., 1.);
- auto rnd = array<double, 4>();
+ auto rnd = std::array<double, 4>();
for (auto& it : rnd) it = distr(RNG_);
// calculate deflection based on particle energy, loss
- auto [mean_dir, final_dir] = get<eSCATTERING>(c->second)->Scatter(
- grammage / 1_g * square(1_cm), vP.getEnergy() / 1_MeV, E_f / 1_MeV, direction,
- rnd);
-
- // TODO: neglect mean direction deflection because Trajectory is a const ref
- (void)mean_dir;
+ [[maybe_unused]] auto deflection = (c->second).scatter->CalculateMultipleScattering(
+ grammage / 1_g * square(1_cm), vP.getEnergy() / 1_MeV, E_f / 1_MeV, rnd);
+ // TODO: multiple scattering is temporary deactivated !!!!!
// update particle direction after continuous loss caused by multiple
// scattering
- auto vec = QuantityVector(final_dir.GetX() * E_f, final_dir.GetY() * E_f,
- final_dir.GetZ() * E_f);
+ auto vec = QuantityVector(direction.GetX() * E_f, direction.GetY() * E_f,
+ direction.GetZ() * E_f);
vP.setMomentum(MomentumVector(vP_dir.getCoordinateSystem(), vec));
}
@@ -103,7 +105,7 @@ namespace corsika::proposal {
// get or build corresponding track integral calculator and solve the
// integral
auto c = getCalculator(vP, calc);
- auto final_energy = get<eDISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
+ auto final_energy = (c->second).disp->UpperLimitTrackIntegral(
vP.getEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
1_MeV;
auto dE = vP.getEnergy() - final_energy;
@@ -126,20 +128,20 @@ namespace corsika::proposal {
// hyper parameter which must be adjusted.
//
auto const energy = vP.getEnergy();
- auto const energy_lim = std::max(
- energy * 0.9, // either 10% relative loss max., or
- get_kinetic_energy_threshold(code) +
- get_mass(code) // energy thresholds globally defined for individual particles
- * 0.99 // need to go 1% below global e-cut to assure removal in
- // ParticleCut. The 1% does not matter since at cut-time the entire
- // energy is removed.
- );
+ auto const energy_lim =
+ std::max(energy * 0.9, // either 10% relative loss max., or
+ get_kinetic_energy_threshold(
+ code) // energy thresholds globally defined for individual particles
+ * 0.9999 // need to go slightly below global e-cut to assure removal
+ // in ParticleCut. This does not matter since at cut-time
+ // the entire energy is removed.
+ );
// solving the track integral for giving energy lim
auto c = getCalculator(vP, calc);
- auto grammage = get<eDISPLACEMENT>(c->second)->SolveTrackIntegral(
- energy / 1_MeV, energy_lim / 1_MeV) *
- 1_g / square(1_cm);
+ auto grammage =
+ (c->second).disp->SolveTrackIntegral(energy / 1_MeV, energy_lim / 1_MeV) * 1_g /
+ square(1_cm);
// return it in distance aequivalent
auto dist = vP.getNode()->getModelProperties().getArclengthFromGrammage(vT, grammage);
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index 26aa7904c6fce84c33046fe6c2505ab5071b407d..019d31270abeb1da6fe2a99121081363762dcac6 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -66,38 +66,41 @@ namespace corsika::proposal {
std::uniform_real_distribution<double> distr(0., 1.);
// sample a interaction-type, loss and component
- auto rates = get<eINTERACTION>(c->second)->Rates(projectile.getEnergy() / 1_MeV);
- auto [type, comp_ptr, v] = get<eINTERACTION>(c->second)->SampleLoss(
+ auto rates =
+ std::get<eINTERACTION>(c->second)->Rates(projectile.getEnergy() / 1_MeV);
+ auto [type, target_hash, v] = std::get<eINTERACTION>(c->second)->SampleLoss(
projectile.getEnergy() / 1_MeV, rates, distr(RNG_));
// Read how much random numbers are required to calculate the secondaries.
// Calculate the secondaries and deploy them on the corsika stack.
- auto rnd =
- vector<double>(get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type));
+ auto rnd = std::vector<double>(
+ std::get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type));
for (auto& it : rnd) it = distr(RNG_);
Point const& place = projectile.getPosition();
CoordinateSystemPtr const& labCS = place.getCoordinateSystem();
- auto point = PROPOSAL::Vector3D(place.getX(labCS) / 1_cm, place.getY(labCS) / 1_cm,
- place.getZ(labCS) / 1_cm);
+ auto point = PROPOSAL::Cartesian3D(
+ place.getX(labCS) / 1_cm, place.getY(labCS) / 1_cm, place.getZ(labCS) / 1_cm);
auto projectile_dir = projectile.getDirection();
auto d = projectile_dir.getComponents(labCS);
- auto direction = PROPOSAL::Vector3D(d.getX().magnitude(), d.getY().magnitude(),
- d.getZ().magnitude());
- auto loss = make_tuple(static_cast<int>(type), point, direction,
- v * projectile.getEnergy() / 1_MeV, 0.);
- auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries(
- projectile.getEnergy() / 1_MeV, loss, *comp_ptr, rnd);
+ auto direction = PROPOSAL::Cartesian3D(d.getX().magnitude(), d.getY().magnitude(),
+ d.getZ().magnitude());
+
+ auto loss = PROPOSAL::StochasticLoss(
+ static_cast<int>(type), v * projectile.getEnergy() / 1_MeV, point, direction,
+ projectile.getTime() / 1_s, 0., projectile.getEnergy() / 1_MeV);
+ auto target = PROPOSAL::Component::component_map->operator[](target_hash);
+ auto sec =
+ std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
for (auto& s : sec) {
- auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
- auto vecProposal = get<PROPOSAL::Loss::DIRECTION>(s);
+ auto E = s.energy * 1_MeV;
+ auto vecProposal = s.direction;
auto vec = QuantityVector(vecProposal.GetX() * E, vecProposal.GetY() * E,
vecProposal.GetZ() * E);
auto p = MomentumVector(labCS, vec);
- auto sec_code =
- convert_from_PDG(static_cast<PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
+ auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
view.addSecondary(
- make_tuple(sec_code, p, projectile.getPosition(), projectile.getTime()));
+ std::make_tuple(sec_code, p, projectile.getPosition(), projectile.getTime()));
}
}
return ProcessReturn::Ok;
@@ -109,7 +112,8 @@ namespace corsika::proposal {
if (canInteract(projectile.getPID())) {
auto c = getCalculator(projectile, calc);
- return get<eINTERACTION>(c->second)->MeanFreePath(projectile.getEnergy() / 1_MeV) *
+ return std::get<eINTERACTION>(c->second)->MeanFreePath(projectile.getEnergy() /
+ 1_MeV) *
1_g / (1_cm * 1_cm);
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
diff --git a/corsika/detail/modules/proposal/ProposalProcessBase.inl b/corsika/detail/modules/proposal/ProposalProcessBase.inl
index 992d24e4d5feefc081ca658e6fd270c534ee5d25..7b07100888bbdb8c7ecd63fdbbc57bdc7668dd4b 100644
--- a/corsika/detail/modules/proposal/ProposalProcessBase.inl
+++ b/corsika/detail/modules/proposal/ProposalProcessBase.inl
@@ -38,7 +38,7 @@ namespace corsika::proposal {
const auto& medium = mediumData(
prop.getMedium(Point(_env.getCoordinateSystem(), 0_cm, 0_cm, 0_cm)));
- auto comp_vec = std::vector<PROPOSAL::Components::Component>();
+ auto comp_vec = std::vector<PROPOSAL::Component>();
const auto& comp = prop.getNuclearComposition();
auto frac_iter = comp.getFractions().cbegin();
for (auto& pcode : comp.getComponents()) {
@@ -54,15 +54,11 @@ namespace corsika::proposal {
}
});
- PROPOSAL::InterpolationDef::order_of_interpolation = 2;
- PROPOSAL::InterpolationDef::nodes_cross_section = 100;
- PROPOSAL::InterpolationDef::nodes_propagate = 1000;
-
//! If corsika data exist store interpolation tables to the corresponding
//! path, otherwise interpolation tables would only stored in main memory if
//! no explicit intrpolation def is specified.
if (auto data_path = std::getenv("CORSIKA_DATA")) {
- PROPOSAL::InterpolationDef::path_to_tables = std::string(data_path) + "/PROPOSAL";
+ PROPOSAL::InterpolationSettings::TABLES_PATH = std::string(data_path) + "/PROPOSAL";
} else {
throw std::runtime_error(
"It is not recommended to run PROPOSAL without its tables in "
diff --git a/corsika/detail/modules/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
index 27fd2a31f43290a26e6c7807f22d29b3debf7e37..163ec37a7864fdd2b2c85022aaafce73e73d6739 100644
--- a/corsika/detail/modules/pythia8/Decay.inl
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -62,8 +62,10 @@ namespace corsika::pythia8 {
// Pythia8::Pythia::particleData.readString("59:m0 = 101.00");
+ // LCOV_EXCL_START, we don't validate pythia8 internals
if (!Pythia8::Pythia::init())
throw std::runtime_error("Pythia::Decay: Initialization failed!");
+ // LCOV_EXCL_STOP
}
inline bool Decay::canHandleDecay(Code const vParticleCode) {
@@ -74,7 +76,7 @@ namespace corsika::pythia8 {
vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
vParticleCode == Code::Electron || vParticleCode == Code::Positron)
return false;
- else if (canDecay(vParticleCode)) // non-zero for particles known to sibyll
+ else if (canDecay(vParticleCode)) // check pythia8 internal
return true;
else
return false;
@@ -207,15 +209,19 @@ namespace corsika::pythia8 {
event.append(pdgCode, 1, 0, 0, // PID, status, col, acol
px, py, pz, en, m);
+ // LCOV_EXCL_START, we don't validate pythia8 internals
if (!Pythia8::Pythia::next())
throw std::runtime_error("Pythia::Decay: decay failed!");
- else
- CORSIKA_LOG_DEBUG("Pythia::Decay: particles after decay: {} ", event.size());
+ // LCOV_EXCL_STOP
+
+ CORSIKA_LOG_DEBUG("Pythia::Decay: particles after decay: {} ", event.size());
+ // LCOV_EXCL_START, we don't validate pythia8 internals
if (print_listing_) {
// list final state
event.list();
}
+ // LCOV_EXCL_STOP
// loop over final state
for (int i = 0; i < event.size(); ++i)
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 4f288918b4215eba6925c2ca70ea6160892628c2..74a910de769ff259aeafacc46e9b2d74d02e32e5 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -45,8 +45,10 @@ namespace corsika::pythia8 {
Pythia8::Pythia::readString("HardQCD:all = on");
Pythia8::Pythia::readString("ProcessLevel:resonanceDecays = off");
+ // we can't test this block, LCOV_EXCL_START
if (!Pythia8::Pythia::init())
throw std::runtime_error("Pythia::Interaction: Initialization failed!");
+ // LCOV_EXCL_STOP
// any decays in pythia? if yes need to define which particles
if (internalDecays_) {
@@ -110,8 +112,10 @@ namespace corsika::pythia8 {
Pythia8::Pythia::readString("Beams:eB = 0.");
// initialize this config
+ // we can't test this block, LCOV_EXCL_START
if (!Pythia8::Pythia::init())
throw std::runtime_error("Pythia::Interaction: Initialization failed!");
+ // LCOV_EXCL_STOP
}
inline bool Interaction::canInteract(Code const pCode) {
@@ -137,9 +141,11 @@ namespace corsika::pythia8 {
return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
- } else
+ } else {
+ // we can't test pythia8 internals, LCOV_EXCL_START
throw std::runtime_error("pythia cross section init failed");
-
+ // we can't test pythia8 internals, LCOV_EXCL_STOP
+ }
} else {
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
std::numeric_limits<double>::infinity() * 1_mb);
@@ -341,17 +347,20 @@ namespace corsika::pythia8 {
configureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
- // create event in pytia
+ // create event in pytia. LCOV_EXCL_START: we don't validate pythia8 internals
if (!Pythia8::Pythia::next())
throw std::runtime_error("Pythia::DoInteraction: failed!");
+ // LCOV_EXCL_STOP
// link to pythia stack
Pythia8::Event& event = Pythia8::Pythia::event;
+ // LCOV_EXCL_START, we don't validate pythia8 internals
if (print_listing_) {
// print final state
event.list();
}
+ // LCOV_EXCL_STOP
MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index 78b237e9484aa3b76dc9395f3bbfe8b24fe1ba40..bbcacf4e5e7e90fca3fa5a21a0aa639848ccdc3c 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -81,8 +81,9 @@ namespace corsika::sibyll {
BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
if (!iBeam)
throw std::runtime_error(
- "Interaction: getCrossSection: interaction of beam hadron not defined in "
- "Sibyll!");
+ fmt::format("Interaction of beam {} not defined in "
+ "Sibyll!",
+ BeamId));
if (!isValidCoMEnergy(CoMenergy)) {
throw std::runtime_error(
"Interaction: getCrossSection: CoM energy outside range for Sibyll!");
diff --git a/corsika/modules/proposal/ContinuousProcess.hpp b/corsika/modules/proposal/ContinuousProcess.hpp
index c8895c8aa752bcd3a3630d4876ace48578d232f3..bf421186c1e68c2522060f601fbc5116bfc2bc7b 100644
--- a/corsika/modules/proposal/ContinuousProcess.hpp
+++ b/corsika/modules/proposal/ContinuousProcess.hpp
@@ -14,6 +14,7 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/random/UniformRealDistribution.hpp>
@@ -32,11 +33,12 @@ namespace corsika::proposal {
: public corsika::ContinuousProcess<proposal::ContinuousProcess>,
ProposalProcessBase {
- enum { eDISPLACEMENT, eSCATTERING };
- using calc_t = std::tuple<std::unique_ptr<PROPOSAL::Displacement>,
- std::unique_ptr<PROPOSAL::Scattering>>;
+ struct Calculator {
+ std::unique_ptr<PROPOSAL::Displacement> disp;
+ std::unique_ptr<PROPOSAL::Scattering> scatter;
+ };
- std::unordered_map<calc_key_t, calc_t, hash>
+ std::unordered_map<calc_key_t, Calculator, hash>
calc; //!< Stores the displacement and scattering calculators.
HEPEnergyType energy_lost_ = 0 * electronvolt;
diff --git a/corsika/modules/proposal/Interaction.hpp b/corsika/modules/proposal/Interaction.hpp
index b492862f0b091aea23251ba449821caa4a8faf54..80a8a7864247e2bc798db5564aa8506869ce1525 100644
--- a/corsika/modules/proposal/Interaction.hpp
+++ b/corsika/modules/proposal/Interaction.hpp
@@ -12,6 +12,7 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/random/UniformRealDistribution.hpp>
@@ -31,8 +32,8 @@ namespace corsika::proposal {
class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
enum { eSECONDARIES, eINTERACTION };
- using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
- unique_ptr<PROPOSAL::Interaction>>;
+ using calculator_t = std::tuple<std::unique_ptr<PROPOSAL::SecondariesCalculator>,
+ std::unique_ptr<PROPOSAL::Interaction>>;
std::unordered_map<calc_key_t, calculator_t, hash>
calc; //!< Stores the secondaries and interaction calculators.
diff --git a/corsika/modules/proposal/ProposalProcessBase.hpp b/corsika/modules/proposal/ProposalProcessBase.hpp
index bb9e3b721c84d57fcaa8956042fbea99a221d91f..3edb90720f2346161f062c45ec79b8f9de85c639 100644
--- a/corsika/modules/proposal/ProposalProcessBase.hpp
+++ b/corsika/modules/proposal/ProposalProcessBase.hpp
@@ -48,19 +48,16 @@ namespace corsika::proposal {
corsika::units::si::HEPEnergyType
emCut) { //!< Stochastic losses smaller than the given cut
//!< will be handeled continuously.
- using namespace corsika::units::si;
- auto p_cut =
- std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, true);
- return PROPOSAL::DefaultCrossSections<T>::template Get<std::false_type>(
- T(), m, p_cut, true);
+ auto p_cut = std::make_shared<const PROPOSAL::EnergyCutSettings>(
+ 0.5 * emCut / 1_MeV, 1, false);
+ return PROPOSAL::GetStdCrossSections(T(), m, p_cut, true);
};
//!
//! PROPOSAL default crosssections are maped to corresponding corsika particle
//! code.
//!
- static std::map<Code, std::function<PROPOSAL::crosssection_list_t<PROPOSAL::ParticleDef,
- PROPOSAL::Medium>(
+ static std::map<Code, std::function<PROPOSAL::crosssection_list_t(
PROPOSAL::Medium&, corsika::units::si::HEPEnergyType)>>
cross = {{Code::Photon, cross_builder<PROPOSAL::GammaDef>},
{Code::Electron, cross_builder<PROPOSAL::EMinusDef>},
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index a8916b55e67a245f7aac1f01cdf8d3a3c3a81689..21354adcfc31d8fec2ec28381a7aeb4ceb39e39e 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -160,7 +160,7 @@ int main(int argc, char** argv) {
mass = get_mass(beamCode);
}
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
- double theta = 20.;
+ double theta = 0.;
double phi = 180.;
auto const thetaRad = theta / 180. * M_PI;
auto const phiRad = phi / 180. * M_PI;
@@ -396,9 +396,10 @@ int main(int argc, char** argv) {
// register the observation plane with the output
output.add("particles", observationLevel);
- auto sequence =
- make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
- emContinuous, cut, trackWriter, observationLevel, longprof);
+ auto sequence = make_sequence( // emCascadeCounted,
+ stackInspect, hadronSequence, reset_particle_mass, decaySequence,
+ // emContinuous,
+ BetheBlochPDG(showerAxis), cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
diff --git a/modules/CMakeLists_PROPOSAL.txt b/modules/CMakeLists_PROPOSAL.txt
deleted file mode 100644
index c487b17cdb880e772e0ac7ccc723f1133fdf71f9..0000000000000000000000000000000000000000
--- a/modules/CMakeLists_PROPOSAL.txt
+++ /dev/null
@@ -1,215 +0,0 @@
-if(CMAKE_PROJECT_NAME STREQUAL corsika)
- message(STATUS "Including PROPOSAL as part of CORSIKA8")
- set (IN_CORSIKA8 ON)
-else(CMAKE_PROJECT_NAME STREQUAL corsika)
- set (IN_CORSIKA8 OFF)
-
- cmake_minimum_required(VERSION 3.8)
- project(PROPOSAL VERSION 6.1.2 LANGUAGES CXX)
-
-endif(CMAKE_PROJECT_NAME STREQUAL corsika)
-
-project(PROPOSAL VERSION 6.1.2 LANGUAGES CXX)
-
-if (NOT IN_CORSIKA8)
-IF(APPLE)
- # In newer version of cmake this will be the default
- SET(CMAKE_MACOSX_RPATH 1)
-ENDIF(APPLE)
-
-# sets standard installtion paths
-include(GNUInstallDirs)
-
-### full RPATH
-### copied from https://cmake.org/Wiki/CMake_RPATH_handling
-### set the RPATH so that for using PROPOSAL in python
-### the DYLD_LIBRARY_PATH must not be set in the bashrc
-### But for using PROPOSAL as c-Library, this path still
-### has to be set
-
-# use, i.e. don't skip the full RPATH for the build tree
-SET(CMAKE_SKIP_BUILD_RPATH FALSE)
-
-# when building, don't use the install RPATH already
-# (but later on when installing)
-OPTION(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
-list(APPEND CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib;${CMAKE_INSTALL_PREFIX}/lib64")
-message(STATUS "CMAKE_INSTALL_RPATH=${CMAKE_INSTALL_RPATH}")
-
-# add the automatically determined parts of the RPATH
-# which point to directories outside the build tree to the install RPATH
-OPTION(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
-### end full RPATH
-
-# Some additional options
-OPTION(ADD_PYTHON "Choose to compile the python wrapper library" ON)
-endif (NOT IN_CORSIKA8)
-
-#################################################################
-#################### PROPOSAL ########################
-#################################################################
-file(GLOB_RECURSE SRC_FILES ${PROJECT_SOURCE_DIR}/src/PROPOSAL/*)
-if (IN_CORSIKA8)
- add_library(PROPOSAL STATIC ${SRC_FILES})
-else (IN_CORSIKA8)
- add_library(PROPOSAL SHARED ${SRC_FILES})
-endif (IN_CORSIKA8)
-
-add_library(PROPOSAL::PROPOSAL ALIAS PROPOSAL)
-#target_compile_features(PROPOSAL PUBLIC cxx_std_11)
-#set_target_properties(PROPOSAL PROPERTIES CXX_EXTENSIONS OFF)
-target_include_directories(
- PROPOSAL PUBLIC
- $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
- $<INSTALL_INTERFACE:include>
-)
-target_compile_options(PROPOSAL PRIVATE -Wall -Wextra -Wnarrowing -Wpedantic -fdiagnostics-show-option -Wno-format-security)
-if (NOT IN_CORSIKA8)
-install(
- TARGETS PROPOSAL
- EXPORT PROPOSALTargets
- LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
- ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
-)
-else (NOT IN_CORSIKA8)
-install(
- TARGETS PROPOSAL
- EXPORT CORSIKA8PublicTargets
- LIBRARY DESTINATION lib/corsika
- ARCHIVE DESTINATION lib/corsika
-)
-endif (NOT IN_CORSIKA8)
-
-# input version from the project call, so the library knows its own version
-configure_file(
- "${PROJECT_SOURCE_DIR}/include/PROPOSAL/version.h.in"
- "${PROJECT_BINARY_DIR}/include/PROPOSAL/version.h"
-)
-if (NOT IN_CORSIKA8)
-install(
- FILES ${PROJECT_BINARY_DIR}/include/PROPOSAL/version.h
- DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/include/PROPOSAL
-)
-target_include_directories(
- PROPOSAL PUBLIC
- $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
- $<INSTALL_INTERFACE:include>
-)
-else (NOT IN_CORSIKA8)
-install(
- FILES ${PROJECT_BINARY_DIR}/include/PROPOSAL/version.h
- DESTINATION include/corsika_modules/PROPOSAL
-)
-target_include_directories(
- PROPOSAL PUBLIC
- $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
- $<INSTALL_INTERFACE:include/corsika_modules>
-)
-endif (NOT IN_CORSIKA8)
-
-# install header files
-file(GLOB_RECURSE INC_FILES ${PROJECT_SOURCE_DIR}/include/*)
-if (NOT IN_CORSIKA8)
-foreach(INC_FILE ${INC_FILES})
- file(RELATIVE_PATH REL_FILE ${PROJECT_SOURCE_DIR}/include ${INC_FILE})
- get_filename_component(DIR ${REL_FILE} DIRECTORY)
- install(FILES include/${REL_FILE} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${DIR})
-endforeach()
-else (NOT IN_CORSIKA8)
-foreach(INC_FILE ${INC_FILES})
- file(RELATIVE_PATH REL_FILE ${PROJECT_SOURCE_DIR}/include/ ${INC_FILE})
- get_filename_component(DIR ${REL_FILE} DIRECTORY)
- install(FILES include/${REL_FILE} DESTINATION include/corsika_modules/${DIR})
-endforeach()
-endif (NOT IN_CORSIKA8)
-
-
-if(NOT IS_SYMLINK ${CMAKE_BINARY_DIR}/resources)
- execute_process(COMMAND ln -sv ${CMAKE_CURRENT_SOURCE_DIR}/resources ${CMAKE_BINARY_DIR}/resources OUTPUT_VARIABLE link_msg OUTPUT_STRIP_TRAILING_WHITESPACE)
- message(STATUS "Symlink to resources created:")
- message(STATUS " ${link_msg}")
-endif()
-
-
-#################################################################
-################# spdlog #######################
-#################################################################
-
-if (NOT IN_CORSIKA8)
-set(CMAKE_POSITION_INDEPENDENT_CODE ON)
-add_subdirectory("vendor/spdlog")
-target_link_libraries(PROPOSAL PRIVATE spdlog::spdlog)
-else (NOT IN_CORSIKA8)
-target_link_libraries(PROPOSAL PRIVATE CONAN_PKG::spdlog)
-endif (NOT IN_CORSIKA8)
-
-
-#################################################################
-################# Tests ########################
-#################################################################
-
-if (NOT IN_CORSIKA8)
-if(CMAKE_PROJECT_NAME STREQUAL PROPOSAL)
- include(CTest)
-endif()
-if(CMAKE_PROJECT_NAME STREQUAL PROPOSAL AND BUILD_TESTING)
- add_subdirectory(tests)
-else()
- MESSAGE(STATUS "No tests will be build.")
-endif()
-endif (NOT IN_CORSIKA8)
-
-#################################################################
-################# Documentation ########################
-#################################################################
-if (NOT IN_CORSIKA8)
-add_subdirectory(doc)
-endif (NOT IN_CORSIKA8)
-
-
-#################################################################
-################# python #########################
-#################################################################
-if (NOT IN_CORSIKA8)
-IF(ADD_PYTHON)
- message(STATUS "Building the python wrapper library.")
- find_package(PythonLibs REQUIRED)
- add_subdirectory("vendor/pybind11")
- file(GLOB_RECURSE PYTHON_SRC_FILES ${PROJECT_SOURCE_DIR}/interface/python/*)
-
- pybind11_add_module(pyproposal SHARED ${PYTHON_SRC_FILES})
- set_target_properties(pyproposal PROPERTIES OUTPUT_NAME proposal)
- target_include_directories(pyproposal PRIVATE ${PYTHON_INCLUDE_DIRS})
- target_include_directories(pyproposal PRIVATE ${PROJECT_SOURCE_DIR}/interface/python)
- target_link_libraries(pyproposal PRIVATE PROPOSAL)
- target_compile_options(pyproposal PRIVATE -fvisibility=hidden)
- install(TARGETS pyproposal EXPORT PROPOSALTargets DESTINATION ${CMAKE_INSTALL_LIBDIR})
-ELSE(ADD_PYTHON)
- MESSAGE(STATUS "No python wrapper library will be build.")
-ENDIF(ADD_PYTHON)
-endif (NOT IN_CORSIKA8)
-
-
-
-#################################################################
-################# INSTALLATION ###################
-#################################################################
-if (NOT IN_CORSIKA8)
-install(
- EXPORT PROPOSALTargets
- FILE PROPOSALConfig.cmake
- NAMESPACE PROPOSAL::
- DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/PROPOSAL
-)
-
-include(CMakePackageConfigHelpers)
-write_basic_package_version_file(
- "PROPOSALConfigVersion.cmake"
- VERSION ${PACKAGE_VERSION}
- COMPATIBILITY SameMajorVersion
-)
-install(
- FILES "${CMAKE_CURRENT_BINARY_DIR}/PROPOSALConfigVersion.cmake"
- DESTINATION lib/cmake/PROPOSAL
-)
-endif (NOT IN_CORSIKA8)
diff --git a/modules/proposal b/modules/proposal
deleted file mode 160000
index 92b9b9ca20826725cb805d135a1a5d5b29a7e06a..0000000000000000000000000000000000000000
--- a/modules/proposal
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 92b9b9ca20826725cb805d135a1a5d5b29a7e06a
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index b2d683a88bbcff787753d434f46b4cc601657ab1..846c734877badaaa8cbb508f30da1716ed24b78d 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -20,7 +20,7 @@
using namespace corsika;
-TEST_CASE("Pythia", "[processes]") {
+TEST_CASE("Pythia8", "modules") {
logging::set_level(logging::level::info);
@@ -64,7 +64,7 @@ TEST_CASE("Pythia", "[processes]") {
std::set<Code> const particleList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
Code::KMinus, Code::K0Long, Code::K0Short};
RNGManager::getInstance().registerRandomStream("pythia");
- corsika::pythia8::Decay model(particleList);
+ corsika::pythia8::Decay decay(particleList);
}
}
@@ -91,10 +91,9 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
return sum;
}
-TEST_CASE("PythiaInterface", "[processes]") {
+TEST_CASE("Pythia8Interface", "modules") {
logging::set_level(logging::level::info);
-
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
@@ -115,14 +114,12 @@ TEST_CASE("PythiaInterface", "[processes]") {
RNGManager::getInstance().registerRandomStream("pythia");
- corsika::pythia8::Decay model(particleList);
+ corsika::pythia8::Decay decay(particleList);
- CORSIKA_LOG_INFO("stack: {} {}", stack.asString(), particle.asString());
- [[maybe_unused]] const TimeType time = model.getLifetime(particle);
+ [[maybe_unused]] const TimeType time = decay.getLifetime(particle);
double const gamma = particle.getEnergy() / get_mass(Code::PiPlus);
CHECK(time == get_lifetime(Code::PiPlus) * gamma);
- model.doDecay(*secViewPtr);
- CORSIKA_LOG_INFO("piplus->{}", stack.asString());
+ decay.doDecay(*secViewPtr);
CHECK(stack.getEntries() == 3); // piplus, muplu, numu
auto const pSum = sumMomentum(view, cs);
CHECK((pSum - plab).getNorm() / 1_GeV == Approx(0).margin(1e-4));
@@ -130,26 +127,33 @@ TEST_CASE("PythiaInterface", "[processes]") {
}
SECTION("pythia decay config") {
- corsika::pythia8::Decay model({Code::PiPlus, Code::PiMinus});
- CHECK(model.isDecayHandled(Code::PiPlus));
- CHECK(model.isDecayHandled(Code::PiMinus));
- CHECK_FALSE(model.isDecayHandled(Code::KPlus));
+ corsika::pythia8::Decay decay({Code::PiPlus, Code::PiMinus});
+ CHECK(decay.isDecayHandled(Code::PiPlus));
+ CHECK(decay.isDecayHandled(Code::PiMinus));
+ CHECK_FALSE(decay.isDecayHandled(Code::KPlus));
const std::vector<Code> particleTestList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
Code::Lambda0Bar, Code::D0Bar};
// setup decays
- model.setHandleDecay(particleTestList);
- for (auto& pCode : particleTestList) CHECK(model.isDecayHandled(pCode));
+ decay.setHandleDecay(particleTestList);
+ for (auto& pCode : particleTestList) CHECK(decay.isDecayHandled(pCode));
// individually
- model.setHandleDecay(Code::KMinus);
+ decay.setHandleDecay(Code::KMinus);
+
+ // impossible
+ CHECK_THROWS(decay.setHandleDecay(Code::Photon));
+
+ CHECK(decay.isDecayHandled(Code::PiPlus));
+ CHECK_FALSE(decay.isDecayHandled(Code::Photon));
// possible decays
- CHECK_FALSE(model.canHandleDecay(Code::Proton));
- CHECK_FALSE(model.canHandleDecay(Code::Electron));
- CHECK(model.canHandleDecay(Code::PiPlus));
- CHECK(model.canHandleDecay(Code::MuPlus));
+ CHECK_FALSE(decay.canHandleDecay(Code::Photon));
+ CHECK_FALSE(decay.canHandleDecay(Code::Proton));
+ CHECK_FALSE(decay.canHandleDecay(Code::Electron));
+ CHECK(decay.canHandleDecay(Code::PiPlus));
+ CHECK(decay.canHandleDecay(Code::MuPlus));
}
SECTION("pythia interaction") {
@@ -161,10 +165,81 @@ TEST_CASE("PythiaInterface", "[processes]") {
auto& view = *secViewPtr;
auto particle = stackPtr->first();
- corsika::pythia8::Interaction model;
- model.doInteraction(view);
- [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
+ corsika::pythia8::Interaction collision;
+
+ CHECK(collision.canInteract(Code::Proton));
+ CHECK(collision.canInteract(Code::AntiProton));
+ CHECK(collision.canInteract(Code::Neutron));
+ CHECK(collision.canInteract(Code::AntiNeutron));
+ CHECK(collision.canInteract(Code::PiMinus));
+ CHECK(collision.canInteract(Code::PiPlus));
+ CHECK_FALSE(collision.canInteract(Code::Electron));
+
+ // nuclei not supported
+ CHECK_THROWS(collision.getCrossSection(Code::Proton, Code::Helium, 1_TeV));
+ std::tuple<CrossSectionType, CrossSectionType> xs_test =
+ collision.getCrossSection(Code::Iron, Code::Hydrogen, 1_GeV);
+ CHECK(std::get<0>(xs_test) == std::numeric_limits<double>::infinity() * 1_mb);
+ CHECK(std::get<1>(xs_test) == std::numeric_limits<double>::infinity() * 1_mb);
+
+ collision.getInteractionLength(particle);
+
+ collision.doInteraction(view);
+ [[maybe_unused]] const GrammageType length = collision.getInteractionLength(particle);
CHECK(length / 1_kg * square(1_m) == Approx(43.04).margin(5e-1));
CHECK(view.getSize() == 38);
}
+
+ SECTION("pythia nucleus projectile") {
+
+ // this is a projectile nucleus with very little energy
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Oxygen, 0, 0, 17_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ auto& view = *secViewPtr;
+ auto particle = stackPtr->first();
+
+ corsika::pythia8::Interaction collision;
+
+ GrammageType lambda_test = collision.getInteractionLength(particle);
+ CHECK(lambda_test == std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm));
+
+ CHECK_THROWS(collision.doInteraction(view));
+ }
+
+ SECTION("pythia too low energy") {
+
+ // this is a projectile neutron with very little energy
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Neutron, 0, 0, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ auto& view = *secViewPtr;
+ auto particle = stackPtr->first();
+
+ corsika::pythia8::Interaction collision;
+
+ GrammageType lambda_test = collision.getInteractionLength(particle);
+ CHECK(lambda_test == std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm));
+
+ CHECK_THROWS(collision.doInteraction(view));
+ }
+
+ SECTION("pythia wrong target") {
+
+ // incompatible target
+ auto [env_Fe, csPtr_Fe, nodePtr_Fe] = setup::testing::setup_environment(Code::Iron);
+ [[maybe_unused]] auto const& cs_Fe = *csPtr_Fe;
+ [[maybe_unused]] auto const& env_dummy_Fe = env_Fe;
+ [[maybe_unused]] auto const& node_dummy_Fe = nodePtr_Fe;
+
+ // resonable projectile, but tool low energy
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Proton, 0, 0, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr_Fe,
+ *csPtr_Fe);
+ auto& view = *secViewPtr;
+
+ corsika::pythia8::Interaction collision;
+
+ CHECK_THROWS(collision.doInteraction(view));
+ }
}