diff --git a/CMakeLists.txt b/CMakeLists.txt
index 344257e5b308456242d15444b4853dc97b3fc96e..5432ca03c7c6ee14036ac64ffc3d09d57cad51a8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -258,6 +258,7 @@ add_subdirectory (modules/qgsjetII)
add_subdirectory (modules/urqmd)
add_subdirectory (modules/conex)
add_subdirectory (modules/epos)
+add_subdirectory (modules/fluka)
#+++++++++++++++++++++++++++++++
# unit testing
diff --git a/README.md b/README.md
index 13e840d6a54083e55d3268612164d24b84aa427f..a41b18bf83481f705c1b863723e96f4dc683a05a 100644
--- a/README.md
+++ b/README.md
@@ -63,6 +63,7 @@ You will also need:
- cmake
- git
- g++, gfortran, binutils, make
+- optional: FLUKA (see below)
On a bare Ubuntu 20.04, just add:
``` shell
@@ -98,6 +99,17 @@ make -j8
make install
```
+### FLUKA support
+For legal reasons we do not distribute/bundle FLUKA together with CORSIKA 8.
+If you want to use FLUKA as low-energy hadronic interaction model, you have to download
+it separately from (http://www.fluka.org/), which requires registering there as FLUKA user.
+You need to download binaries suitable for your system (check compiler and glibc versions)
+and unpack the binaries into a directory of your choice. You also need the data archive
+(fluka20xy.z-data.tar.gz) unpacked in the same directory. To compile CORSIKA 8 with FLUKA,
+you need to specify the path of libflukahp.a when invoking cmake with the option
+`-DC8_FLUKALIB=<path>`. CMake will print a status message indicating whether FLUKA support
+is enabled or disabled when the library is (not) found.
+
## Installation (using docker containers)
diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 06bca17cf981b5f8b96300660c3aad713c7b0dc1..b5b4410aa60c1a0c5e8aed7ad4b4b0a8f159ef11 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -132,6 +132,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);
}
}
diff --git a/corsika/detail/framework/geometry/CoordinateSystem.inl b/corsika/detail/framework/geometry/CoordinateSystem.inl
index 12e914ba9d4131209d4be8fa7eec066e9961a97a..f8297585599d2b448a2daa2e00fa3e3f4132409b 100644
--- a/corsika/detail/framework/geometry/CoordinateSystem.inl
+++ b/corsika/detail/framework/geometry/CoordinateSystem.inl
@@ -81,7 +81,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;
diff --git a/corsika/detail/modules/fluka/InteractionModel.inl b/corsika/detail/modules/fluka/InteractionModel.inl
new file mode 100644
index 0000000000000000000000000000000000000000..eb3f2473bac1f244953312629fa7bcd0bb0726a7
--- /dev/null
+++ b/corsika/detail/modules/fluka/InteractionModel.inl
@@ -0,0 +1,238 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <algorithm>
+#include <cstdlib>
+#include <vector>
+#include <iterator>
+#include <set>
+#include <utility>
+
+#include <boost/iterator/zip_iterator.hpp>
+#include <Eigen/Dense>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/fluka/ParticleConversion.hpp>
+
+#include <FLUKA.hpp>
+
+namespace corsika::fluka {
+ template <typename TEnvironment>
+ inline InteractionModel::InteractionModel(TEnvironment const& env)
+ : materials_{genFlukaMaterials(env)}
+ , cumsgx_{std::make_unique<double[]>(materials_.size() * 3)} {
+ for (auto const& [code, matno] : materials_) {
+ CORSIKA_LOGGER_DEBUG(logger_, "FLUKA material initialization: {} -> {}",
+ get_name(code, full_name{}), matno);
+ }
+
+ if (int const ndmhep = ::fluka::ndmhep_(); ::fluka::nmxhep != ndmhep) {
+ CORSIKA_LOGGER_CRITICAL(logger_, "HEPEVT dimension mismatch. FLUKA reports %d",
+ ndmhep);
+ throw std::runtime_error{"FLUKA HEPEVT dimension mismatch"};
+ }
+ }
+
+ inline bool InteractionModel::isValid(Code projectileID, int material,
+ HEPEnergyType /*sqrtS*/) const {
+ if (!fluka::canInteract(projectileID)) {
+ // invalid projectile
+ return false;
+ }
+
+ if (material < 0) {
+ // invalid/uninitialized target
+ return false;
+ }
+
+ // TODO: check validity range of sqrtS
+ return true;
+ }
+
+ inline bool InteractionModel::isValid(Code projectileID, Code targetID,
+ HEPEnergyType sqrtS) const {
+ return isValid(projectileID, getMaterialIndex(targetID), sqrtS);
+ }
+
+ inline int InteractionModel::getMaterialIndex(Code targetID) const {
+ auto compare = [=](std::pair<Code, int> const& p) { return p.first == targetID; };
+ if (auto it = std::find_if(materials_.cbegin(), materials_.cend(), compare);
+ it == materials_.cend()) {
+ return -1;
+ } else {
+ return it->second;
+ }
+ }
+
+ inline CrossSectionType InteractionModel::getCrossSection(
+ Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ auto const flukaMaterial = getMaterialIndex(targetId);
+
+ HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+ if (!isValid(projectileId, flukaMaterial, sqrtS)) { return CrossSectionType::zero(); }
+
+ COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+ FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+ HEPEnergyType const Elab = projectileLab4mom.getTimeLikeComponent();
+ auto constexpr invGeV = 1 / 1_GeV;
+ double const labMomentum =
+ projectileLab4mom.getSpaceLikeComponents().getNorm() * invGeV;
+
+ auto const plab = projectileLab4mom.getSpaceLikeComponents();
+
+ CORSIKA_LOGGER_DEBUG(logger_, fmt::format("Elab = {} GeV", Elab * invGeV));
+ auto const flukaCodeProj =
+ static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+
+ double const dummyEkin = 0;
+ CrossSectionType const xs = ::fluka::sgmxyz_(&flukaCodeProj, &flukaMaterial,
+ &dummyEkin, &labMomentum, &iflxyz_) *
+ 1_mb;
+ return xs;
+ }
+
+ template <typename TSecondaryView>
+ inline void InteractionModel::doInteraction(TSecondaryView& view,
+ Code const projectileId,
+ Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) {
+
+ auto const flukaCodeProj =
+ static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+ auto const flukaMaterial = getMaterialIndex(targetId);
+
+ HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+ if (!isValid(projectileId, flukaMaterial, sqrtS)) {
+ std::string const errmsg = fmt::format(
+ "Event generation with invalid configuration requested: proj: {}, target: {}",
+ get_name(projectileId, full_name{}), get_name(targetId, full_name{}));
+ CORSIKA_LOGGER_CRITICAL(logger_, errmsg);
+ throw std::runtime_error{errmsg.c_str()};
+ }
+
+ COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+ FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+ auto constexpr invGeV = 1 / 1_GeV;
+
+ auto const plab = projectileLab4mom.getSpaceLikeComponents();
+ auto const& cs = plab.getCoordinateSystem();
+ auto const labMomentum = plab.getNorm();
+ double const labMomentumGeV = labMomentum * invGeV;
+
+ auto const direction = (plab / labMomentum).getComponents().getEigenVector();
+
+ double const dummyEkin = 0;
+ ::fluka::evtxyz_(&flukaCodeProj, &flukaMaterial, &dummyEkin, &labMomentumGeV,
+ &direction[0], &direction[1], &direction[2], &iflxyz_, cumsgx_.get(),
+ cumsgx_.get() + materials_.size(),
+ cumsgx_.get() + materials_.size() * 2);
+
+ for (int i = 0; i < ::fluka::hepevt_.nhep; ++i) {
+ int const status = ::fluka::hepevt_.isthep[i];
+ if (status != 1) // skip non-final-state particles
+ continue;
+
+ auto const pdg = static_cast<corsika::PDGCode>(::fluka::hepevt_.idhep[i]);
+ auto const c8code = corsika::convert_from_PDG(pdg);
+ auto const mom = QuantityVector<hepenergy_d>{
+ Eigen::Map<Eigen::Vector3d>(&::fluka::hepevt_.phep[i][0]) *
+ (1_GeV).magnitude()};
+ auto const pPrime = mom.getNorm();
+ auto const c8mass = corsika::get_mass(c8code);
+
+ auto const fourMomCollisionFrame =
+ FourVector{calculate_total_energy(pPrime, c8mass), MomentumVector{cs, mom}};
+ auto const fourMomOrigFrame = targetRestBoost.fromCoM(fourMomCollisionFrame);
+ auto const momOrigFrame = fourMomOrigFrame.getSpaceLikeComponents();
+ auto const p = momOrigFrame.getNorm();
+
+ view.addSecondary(std::tuple{c8code, corsika::calculate_kinetic_energy(p, c8mass),
+ momOrigFrame / p});
+ }
+ }
+
+ template <typename TEnvironment>
+ inline std::vector<std::pair<Code, int>> InteractionModel::genFlukaMaterials(
+ TEnvironment const& env) {
+ auto const& universe = *(env.getUniverse());
+
+ // generate complete list of all nuclei types in universe
+ auto const allElementsInUniverse = std::invoke([&]() {
+ std::set<Code> allElementsInUniverse;
+ auto collectElements = [&](auto& vtn) {
+ if (vtn.hasModelProperties()) {
+ auto const& comp =
+ vtn.getModelProperties().getNuclearComposition().getComponents();
+ for (auto const c : comp) allElementsInUniverse.insert(c);
+ }
+ };
+ universe.walk(collectElements);
+ return allElementsInUniverse;
+ });
+
+ /*
+ * We define one material per element/isotope we have in C8. Cross-section averaging
+ * and target isotope selection happen in C8.
+ */
+
+ int const nElements = allElementsInUniverse.size();
+ auto nelmfl = std::make_unique<int[]>(nElements);
+ std::vector<int> izelfl;
+ izelfl.reserve(nElements);
+ auto wfelml = std::make_unique<double[]>(nElements);
+ auto const mxelfl = nElements;
+ double const pptmax = 1e11; // GeV
+ double const ef2dp3 = 0; // GeV, 0 means default is used
+ double const df2dp3 = -1; // default
+ bool const lprint = true;
+ auto mtflka = std::make_unique<int[]>(mxelfl);
+ // magic number that FLUKA uses to see if it's the right version
+ char const crvrck[] = "76466879";
+ int const size = 8;
+
+ std::fill(&nelmfl[0], &nelmfl[nElements], 1);
+ std::fill(&wfelml[0], &wfelml[nElements], 1.);
+ std::transform(allElementsInUniverse.cbegin(), allElementsInUniverse.cend(),
+ std::back_inserter(izelfl), get_nucleus_Z);
+
+ // check if FLUPRO is set
+ if (std::getenv("FLUPRO") == nullptr) {
+ throw std::runtime_error{"FLUPRO not set; required to initialize FLUKA"};
+ }
+
+ // call FLUKA
+ ::fluka::stpxyz_(&nElements, nelmfl.get(), izelfl.data(), wfelml.get(), &mxelfl,
+ &pptmax, &ef2dp3, &df2dp3, &iflxyz_, &lprint, mtflka.get(), crvrck,
+ &size);
+
+ // now create & fill vector of (C8 Code, FLUKA mat. no.) pairs
+ std::vector<std::pair<Code, int>> mapping;
+ mapping.reserve(nElements);
+
+ auto it = boost::make_zip_iterator(
+ boost::make_tuple(allElementsInUniverse.begin(), &mtflka[0]));
+ auto const end = boost::make_zip_iterator(
+ boost::make_tuple(allElementsInUniverse.end(), &mtflka[nElements]));
+ for (; it != end; ++it) {
+ boost::tuple<Code const&, int&> tup = *it;
+ mapping.emplace_back(tup.get<0>(), tup.get<1>());
+ }
+
+ return mapping;
+ }
+} // namespace corsika::fluka
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index 76cbe20d7dbb1799c6c39c0d7265ffda2c34a881..f797debd2f4d909951c0684f5028face4d1f847c 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -86,7 +86,7 @@ namespace corsika {
//! internal method
void setBoost(double const coshEta, double const sinhEta);
- private:
+ public:
Eigen::Matrix2d boost_;
Eigen::Matrix2d inverseBoost_;
CoordinateSystemPtr const originalCS_;
diff --git a/corsika/modules/Epos.hpp b/corsika/modules/Epos.hpp
index b4aa9086cf563e35b0ffb23d4aa9c85f6375cf60..0fb3db62407f93b713048182faa6a82e870dab72 100644
--- a/corsika/modules/Epos.hpp
+++ b/corsika/modules/Epos.hpp
@@ -13,7 +13,7 @@
#include <corsika/framework/process/InteractionProcess.hpp>
/**
- * @file Sibyll.hpp
+ * @file Epos.hpp
*
* Includes all the parts of the EPOS model. Defines the InteractionProcess<TModel>
* classes needed for the ProcessSequence.
@@ -27,4 +27,4 @@ namespace corsika::epos {
* to provide all the functions for ProcessSequence.
*/
class Interaction : public InteractionModel, public InteractionProcess<Interaction> {};
-} // namespace corsika::epos
\ No newline at end of file
+} // namespace corsika::epos
diff --git a/corsika/modules/FLUKA.hpp b/corsika/modules/FLUKA.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..43f4e32e30ebf2d05d7f2edaca9b9c6b41bdf700
--- /dev/null
+++ b/corsika/modules/FLUKA.hpp
@@ -0,0 +1,35 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/fluka/InteractionModel.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+/**
+ * @file FLUKA.hpp
+ *
+ * Includes all the parts of the FLUKA model. Defines the InteractionProcess<TModel>
+ * classes needed for the ProcessSequence.
+ */
+
+namespace corsika::fluka {
+ /**
+ * fluka::Interaction is the process for ProcessSequence.
+ *
+ * The fluka::InteractionModel is wrapped as an InteractionProcess here in order
+ * to provide all the functions for ProcessSequence.
+ */
+ class Interaction : public fluka::InteractionModel,
+ public corsika::InteractionProcess<Interaction> {
+ public:
+ template <typename TEnvironment>
+ Interaction(TEnvironment const& env)
+ : fluka::InteractionModel{env} {}
+ };
+} // namespace corsika::fluka
diff --git a/corsika/modules/Random.hpp b/corsika/modules/Random.hpp
index a81bc1ad2d262f0cc049067d58aaa13df372b6b9..60d8b9a9a03a8d35a6655dc24e7ad2894b4e34b5 100644
--- a/corsika/modules/Random.hpp
+++ b/corsika/modules/Random.hpp
@@ -22,3 +22,4 @@
#include <corsika/modules/urqmd/Random.hpp>
#include <corsika/modules/qgsjetII/Random.hpp>
#include <corsika/modules/conex/Random.hpp>
+#include <corsika/modules/fluka/Random.hpp>
diff --git a/corsika/modules/epos/ParticleConversion.hpp b/corsika/modules/epos/ParticleConversion.hpp
index 8e501d67f544c3ac6fc0bd23c247d21e9f5e1680..7faf53da628a4a72c0abfd42432cf23f3e16f5fb 100644
--- a/corsika/modules/epos/ParticleConversion.hpp
+++ b/corsika/modules/epos/ParticleConversion.hpp
@@ -8,6 +8,10 @@
#pragma once
+#include <array>
+#include <cstdint>
+#include <type_traits>
+
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
diff --git a/corsika/modules/fluka/InteractionModel.hpp b/corsika/modules/fluka/InteractionModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8796e48b8e8c4fbb018d6002ce944dd4ea970a2b
--- /dev/null
+++ b/corsika/modules/fluka/InteractionModel.hpp
@@ -0,0 +1,73 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <vector>
+#include <utility>
+#include <memory>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+
+namespace corsika::fluka {
+ /**
+ * This class exposes the (hadronic) interactions of FLUKA. FLUKA needs to be
+ * initialized with a predefined set of target materials and a flag describing the type
+ * of interactions (elastic, inelastic, electromagnetic dissociation). Currently, only
+ * inelastic events are supported.
+ */
+ class InteractionModel {
+ public:
+ /**
+ * Create a new InteractionModel. The FLUKA materials are collected from the elements
+ * present in the environment. Each element is its own FLUKA material, no FLUKA
+ * compounds are used.
+ */
+ template <typename TEnvironment>
+ InteractionModel(TEnvironment const&);
+
+ //! Return the cross-section of a given combination of projectile/target.
+ CrossSectionType getCrossSection(Code projectileId, Code targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
+
+ bool isValid(Code projectileID, Code targetID, HEPEnergyType sqrtS) const;
+ bool isValid(Code projectileID, int material, HEPEnergyType sqrtS) const;
+
+ //! convert target Code to FLUKA material number
+ int getMaterialIndex(Code targetID) const;
+
+ /**
+ * Perform an interaction. Since FLUKA expects a fixed-target configuration, we
+ * perform a Lorentz transform into the rest frame of the target.
+ */
+ template <typename TSecondaryView>
+ void doInteraction(TSecondaryView& view, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+ private:
+ default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("fluka");
+ std::vector<std::pair<Code, int>> const
+ materials_; //!< map target Code to FLUKA material no.
+ std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_FLUKA_Interaction");
+ std::unique_ptr<double[]> cumsgx_; //!< dump for evtxyz cumsg*, never read again
+
+ template <typename TEnvironment>
+ static std::vector<std::pair<Code, int>> genFlukaMaterials(TEnvironment const&);
+ };
+
+ inline static int const iflxyz_ = 1;
+} // namespace corsika::fluka
+
+#include <corsika/detail/modules/fluka/InteractionModel.inl>
diff --git a/corsika/modules/fluka/ParticleConversion.hpp b/corsika/modules/fluka/ParticleConversion.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..228d568ac9a889e08e2f9aa5137b3571ffc8de99
--- /dev/null
+++ b/corsika/modules/fluka/ParticleConversion.hpp
@@ -0,0 +1,50 @@
+/*
+ * (c) Copyright 2023 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 <array>
+#include <cstdint>
+#include <type_traits>
+#include <vector>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <fmt/core.h>
+
+#include <string>
+
+namespace corsika::fluka {
+
+ enum class FLUKACode : int;
+ using FLUKACodeIntType = std::underlying_type<FLUKACode>::type;
+
+#include <corsika/modules/fluka/Generated.inc>
+
+ FLUKACode constexpr convertToFluka(Code const c8id) {
+ if (is_nucleus(c8id)) {
+ throw std::runtime_error{"nucleus conversion to FLUKA not implemented"};
+ }
+
+ FLUKACode const flukaID = corsika2fluka[static_cast<corsika::CodeIntType>(c8id)];
+ if (flukaID == FLUKACode::Unknown) {
+ throw std::runtime_error{
+ fmt::format("no correspondig FLUKA id for {}", get_name(c8id)).c_str()};
+ }
+
+ return flukaID;
+ }
+
+ int constexpr convertToFlukaRaw(Code const code) {
+ return static_cast<FLUKACodeIntType>(convertToFluka(code));
+ }
+
+ bool const canInteract(Code const code) {
+ if (is_nucleus(code)) return false; // nuclei support not yet implemented
+ return flukaCanInteract[static_cast<CodeIntType>(code)];
+ }
+} // namespace corsika::fluka
diff --git a/corsika/modules/fluka/Random.hpp b/corsika/modules/fluka/Random.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..73cf5377a37d0453248b7fdfff0c04f9607ab981
--- /dev/null
+++ b/corsika/modules/fluka/Random.hpp
@@ -0,0 +1,26 @@
+/*
+ * (c) Copyright 2023 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 <random>
+
+#include <corsika/framework/random/RNGManager.hpp>
+
+namespace fluka {
+ /**
+ * the random number generator function for FLUKA
+ */
+ double rndm_interface() {
+ static corsika::default_prng_type& rng =
+ corsika::RNGManager<>::getInstance().getRandomStream("fluka");
+ static std::uniform_real_distribution<double> dist;
+ return dist(rng);
+ }
+
+} // namespace fluka
diff --git a/modules/fluka/CMakeLists.txt b/modules/fluka/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b9cad711bba440c7f9a7d5c6dea34711872768d7
--- /dev/null
+++ b/modules/fluka/CMakeLists.txt
@@ -0,0 +1,55 @@
+set(C8_FLUKALIB CACHE STRING "path to libflukahp.a")
+
+if (EXISTS "${C8_FLUKALIB}")
+ message("libflukahp.a fount at ${C8_FLUKALIB}; FLUKA support is enabled.")
+ set (input_dir ${PROJECT_SOURCE_DIR}/src/modules/fluka)
+
+ # we remove flrndm_.o from the original libflukahp.a and save the result as libflukahp-norndm.a
+ add_custom_command(
+ OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/libflukahp-norndm.a
+ COMMAND ${input_dir}/strip_flukahp.sh ${C8_FLUKALIB} ${CMAKE_CURRENT_BINARY_DIR}
+ )
+ add_custom_target(generate_libfluka-norndm DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/libflukahp-norndm.a)
+
+ add_library(flukahp-norndm STATIC IMPORTED)
+ add_dependencies(flukahp-norndm generate_libfluka-norndm)
+ set_property(TARGET flukahp-norndm PROPERTY IMPORTED_LOCATION ${CMAKE_CURRENT_BINARY_DIR}/libflukahp-norndm.a)
+
+ find_library(MATH_LIBRARY m)
+ target_link_libraries(flukahp-norndm INTERFACE gfortran ${MATH_LIBRARY})
+ target_include_directories (flukahp-norndm INTERFACE
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+ $<INSTALL_INTERFACE:include/corsika_modules/fluka>
+ )
+
+ # build shared library containing FLUKA + our own implementation of flrndm_()
+ # Note: This apparently works only when built as shared lib, not possible so easily as static lib
+ # Note also that this lib is not self-contained as it doesn't contain fluka::rndm_interface(), so
+ # you can't load it at runtime
+ add_library(fluka SHARED fluka.cpp)
+ target_link_libraries(fluka PRIVATE flukahp-norndm)
+
+ target_include_directories (fluka PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+ $<INSTALL_INTERFACE:include/corsika_modules/fluka>
+ )
+
+ install (
+ FILES
+ fluka.hpp
+ DESTINATION include/corsika_modules/fluka
+ )
+
+ install (
+ TARGETS fluka
+ EXPORT CORSIKA8PublicTargets
+ ARCHIVE DESTINATION lib/corsika
+ INCLUDES DESTINATION include/corsika_modules/fluka
+ )
+
+ # add fluka to corsika8 build
+ add_dependencies (CORSIKA8 fluka)
+ target_link_libraries (CORSIKA8 INTERFACE fluka)
+else()
+ message("C8_FLUKALIB not set/invalid, building without FLUKA support")
+endif()
diff --git a/modules/fluka/FLUKA.hpp b/modules/fluka/FLUKA.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..12b8a17562cee175314c8ba6c6ff0a5154eb683b
--- /dev/null
+++ b/modules/fluka/FLUKA.hpp
@@ -0,0 +1,229 @@
+/*
+ * (c) Copyright 2022 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 <cstddef>
+#include <array>
+
+namespace fluka {
+ /**
+ * \fluka fluka::rndm_interface
+ *
+ * this is the random number hook to external packages.
+ *
+ * CORSIKA8, for example, has to provide an implementation of this.
+ **/
+ extern double rndm_interface();
+
+ size_t constexpr nmxhep = 10000;
+ template <typename T>
+ using hepmc_array = std::array<T, nmxhep>;
+
+ extern "C" {
+ /**
+ *---------------------------------------------------------------------------
+ * hep standard event commonblock.
+ *---------------------------------------------------------------------------
+ *---------------------------------------------------------------------------
+ *
+ * nevhep - event number
+ * nhep - number of entries in the event record
+ *
+ * isthep(i) - status code
+ * idhep(i) - particle id (particle data group standard)
+ *
+ * jmohep(1,i) - position of mother particle in list
+ * jmohep(2,i) - position of second mother particle in list
+ * jdahep(1,i) - position of first daughter in list
+ * jdahep(2,i) - position of first daughter in list
+ *
+ * phep(1,i) - p_x momentum in gev/c
+ * phep(2,i) - p_y momentum in gev/c
+ * phep(3,i) - p_z momentum in gev/c
+ * phep(4,i) - energy in gev
+ * phep(5,i) - mass in gev/c**2
+ *
+ * vhep(1,i) - x position of production vertex in mm
+ * vhep(2,i) - y position of production vertex in mm
+ * vhep(3,i) - z position of production vertex in mm
+ * vhep(4,i) - time of production in mm/c
+ */
+
+ extern struct {
+ int nevhep; // event number
+ int nhep; // number of entries
+ hepmc_array<int> isthep; // status code
+ hepmc_array<int> idhep; // PDG particle id
+ hepmc_array<int[2]> jmohep; // position of first, second mother
+ hepmc_array<int[2]> jdahep; // position of first, last daughter
+ hepmc_array<double[5]> phep; // 4-momemtum, mass (GeV)
+ hepmc_array<double[4]> vhep; // vertex, production time in mm
+ } hepevt_;
+
+ /* Iflxyz must be consistent in all calls, to Stpxyz, Sgmxyz, Evtxyz
+ * Iflxyz = 1 -> only inelastic
+ * Iflxyz = 10 -> only elastic
+ * Iflxyz = 11 -> inelastic + elastic
+ * Iflxyz =100 -> only emd
+ * Iflxyz =101 -> inelastic + emd
+ * Iflxyz =110 -> elastic + emd
+ * Iflxyz =111 -> inelastic + elastic + emd
+ */
+
+ /**---------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * *
+ * SeTuP for X/Y/Z interaction types: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 24 October 2022 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 01-Nov-22 by Alfredo Ferrari *
+ * Private *
+ * *
+ * Input variables: *
+ * *
+ * Nmatfl = Number of requested Fluka materials *
+ * Nelmfl(m) = Number of elements in the m_th requested Fluka *
+ * material, Nelmfl(m)=1 means simple element, no *
+ * compound *
+ * Izelfl(n) = Cumulative array of the Z of the elements con- *
+ * stituting the requested Fluka materials *
+ * (the Z for the elements of the m_th materials *
+ * start at n = Sum_i=1^m-1 [Nelmfl(i)] + 1) *
+ * Wfelml(n) = Cumulative array of the weight fractions of the*
+ * elements constituting the requested Fluka *
+ * materials (the weight fractions for the elem- *
+ * ents of the m_th materials start at *
+ * n = Sum_i=1^m-1 [Nelmfl(i)] + 1) *
+ * Mxelfl = Dimension of the Iz/Wfelml arrays, it must be *
+ * Mxelfl >= Sum_i=1^nmatfl [Nelmfl(i)] *
+ * Pptmax = Maximum momentum (GeV/c) to be used in initia- *
+ * lization (optional) *
+ * Ef2dp3 = Transition energy (GeV) from Fluka (Peanut) to *
+ * Dpmjet3 for hA interactions (optional) *
+ * Df2dp3 = Smearing (+/-Df2dp3) (GeV) of the transition *
+ * energy from Fluka (Peanut) to Dpmjet3 for hA *
+ * interactions (optional) *
+ * Lprint = Material printout flag *
+ * *
+ * Output variables: *
+ * *
+ * Mtflka(m) = Fluka material number corresponding to the m_th*
+ * requested material *
+ * *
+ *----------------------------------------------------------------------*/
+ void stpxyz_(int const* NMATFL, int const NELMFL[], int const IZELFL[],
+ double const WFELFL[], int const* MXELFL, double const* PPTMAX,
+ double const* EF2DP3, double const* DF2DP3, int const* IFLXYZ,
+ bool const* LPRINT, int* MTFLKA, char const* CRVRCK, int const*);
+
+ /**---------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * *
+ * EVenT for X/Y/Z interaction types: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 24 October 2022 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 28-Oct-22 by Alfredo Ferrari *
+ * Private *
+ * *
+ * Output variables: *
+ * *
+ * Cumsgi(i) = cumulative macroscopic inelastic cross section *
+ * (cm^2/g x rho) for the i_th element of the mmat*
+ * material *
+ * Cumsge(i) = cumulative macroscopic elastic cross section *
+ * (cm^2/g x rho) for the i_th element of the mmat*
+ * material *
+ * Cumsgm(i) = cumulative macroscopic emd cross section *
+ * (cm^2/g x rho) for the i_th element of the mmat*
+ * material *
+ * *
+ * *
+ *----------------------------------------------------------------------*/
+ void evtxyz_(int const* KPROJ, int const* MMAT, double const* EKIN, double const* PPROJ,
+ double const* TXX, double const* TYY, double const* TZZ, int const* IFLXYZ,
+ double CUMSGI[], double CUMSGE[], double CUMSGM[]);
+
+ /**---------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * SiGMa (mb) for X/Y/Z interaction types: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 24 October 2022 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 27-Oct-22 by Alfredo Ferrari *
+ * Private *
+ * *
+ *----------------------------------------------------------------------*/
+ double sgmxyz_(int const* KPROJ, int const* MMAT, double const* EKIN,
+ double const* PPROJ, int const* IFLXYZ);
+
+ /**---------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2023-2023 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * FiLL HEP common: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 06 February 2023 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 07-Feb-23 by Alfredo Ferrari *
+ * Private *
+ * *
+ *----------------------------------------------------------------------*/
+ void fllhep_();
+
+ /**----------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2023-2023 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * N DiMension of the HEP common: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 09 March 2023 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 09-Mar-23 by Alfredo Ferrari *
+ * Private *
+ * *
+ *----------------------------------------------------------------------*/
+ int ndmhep_();
+
+ //! random-number generator called from within FLUKA
+ double flrndm_();
+ }
+} // namespace fluka
diff --git a/modules/fluka/fluka.cpp b/modules/fluka/fluka.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..985bde38bd414cae76dd2113c4df266820578f16
--- /dev/null
+++ b/modules/fluka/fluka.cpp
@@ -0,0 +1,32 @@
+/*
+ * (c) Copyright 2023 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 <FLUKA.hpp>
+
+namespace fluka {
+ double (*rndmPtr)() = &rndm_interface;
+
+ extern "C" {
+ double flrndm_() { return ::fluka::rndmPtr(); }
+
+ //! overwrite function pointer to be used as FLUKA RNG (flrndm_())
+ void setFlukaRNG(double (*func)()) { ::fluka::rndmPtr = func; }
+
+ /**
+ * The following (function) pointers make sure the corresponding objects in libflukahp.a
+ * won't get dropped from the final file during linking.
+ */
+
+ [[maybe_unused]] extern auto* const hepevt_ptr = &hepevt_;
+ [[maybe_unused]] extern auto* const stpxyc_ptr = &stpxyz_;
+ [[maybe_unused]] extern auto* const evtxyz_ptr = &evtxyz_;
+ [[maybe_unused]] extern auto* const sgmxyz_ptr = &sgmxyz_;
+ [[maybe_unused]] extern auto* const ndmhep_ptr = &ndmhep_;
+ }
+
+} // namespace fluka
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index f8d5df55a95bd8e4547d58322fecb6cc23e5f262..34fe8bdd3040ce796ec5a6c435c4fa5680f2c0e0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,6 +1,7 @@
-add_subdirectory (framework/core)
+add_subdirectory (framework/core)
add_subdirectory (media)
-add_subdirectory (modules/sibyll)
-add_subdirectory (modules/sophia)
-add_subdirectory (modules/qgsjetII)
-add_subdirectory (modules/epos)
+add_subdirectory (modules/sibyll)
+add_subdirectory (modules/sophia)
+add_subdirectory (modules/qgsjetII)
+add_subdirectory (modules/epos)
+add_subdirectory (modules/fluka)
diff --git a/src/modules/fluka/CMakeLists.txt b/src/modules/fluka/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2d2d26b447c277b20da5760a7236cab69c334f5e
--- /dev/null
+++ b/src/modules/fluka/CMakeLists.txt
@@ -0,0 +1,26 @@
+set (input_dir ${PROJECT_SOURCE_DIR}/src/modules/fluka)
+set (output_dir ${PROJECT_BINARY_DIR}/corsika/modules/fluka)
+
+file (MAKE_DIRECTORY ${output_dir})
+
+add_custom_command (
+ OUTPUT ${output_dir}/Generated.inc
+ COMMAND ${input_dir}/code_generator.py
+ ${PROJECT_BINARY_DIR}/corsika/framework/core/particle_db.pkl
+ ${input_dir}/fluka_codes.dat
+ DEPENDS ${input_dir}/code_generator.py
+ ${input_dir}/fluka_codes.dat
+ GenParticlesHeaders # for particle_db.pkl
+ WORKING_DIRECTORY
+ ${output_dir}/
+ COMMENT "Generate conversion tables for particle codes FLUKA <-> CORSIKA"
+ VERBATIM
+ )
+
+add_custom_target (SourceDirLinkFLUKA DEPENDS ${output_dir}/Generated.inc)
+add_dependencies (CORSIKA8 SourceDirLinkFLUKA)
+
+install (
+ FILES ${output_dir}/Generated.inc
+ DESTINATION include/corsika/modules/fluka
+)
diff --git a/src/modules/fluka/code_generator.py b/src/modules/fluka/code_generator.py
new file mode 100755
index 0000000000000000000000000000000000000000..3574fa25756193f42f21166088e70506a457ede2
--- /dev/null
+++ b/src/modules/fluka/code_generator.py
@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+
+# (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+#
+# See file AUTHORS for a list of contributors.
+#
+# This software is distributed under the terms of the GNU General Public
+# Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+# the license.
+
+
+import pickle, sys, itertools
+
+
+def load_particledb(filename):
+ '''
+ loads the pickled particle_db (which is an OrderedDict)
+ '''
+ with open(filename, "rb") as f:
+ particle_db = pickle.load(f)
+ return particle_db
+
+
+def read_fluka_codes(filename, particle_db):
+ '''
+ reads to FLUKA codes data file
+
+ For particles that can interact in FLUKA, add can-interact flag
+ '''
+ with open(filename) as f:
+ for line in f:
+ line = line.strip()
+ if len(line) == 0 or line[0].startswith('#'):
+ continue
+ line = line.split()
+ if len(line) == 2:
+ canInteractFlag = "no"
+ identifier, fluka_code = line
+ else:
+ identifier, fluka_code, canInteractFlag = line
+
+ try:
+ particle_db[identifier]["fluka_code"] = int(fluka_code)
+ particle_db[identifier]["fluka_canInteract"] = (canInteractFlag == "yes")
+ except KeyError as e:
+ raise Exception("Identifier '{:s}' not found in CORSIKA8 particle_db".format(identifier))
+
+
+def generate_fluka_enum(particle_db):
+ '''
+ generates the enum to access epos particles by readable names
+ '''
+ output = ("enum class FLUKACode : int {\n"
+ " Unknown = 0,\n")
+ for identifier, pData in particle_db.items():
+ if 'fluka_code' in pData:
+ output += " {:s} = {:d},\n".format(identifier, pData['fluka_code'])
+ output += "};\n"
+ return output
+
+
+def generate_corsika2fluka(particle_db):
+ '''
+ generates the look-up table to convert corsika codes to epos codes
+ '''
+ string = "static inline std::array<FLUKACode, {:d}> constexpr corsika2fluka = {{\n".format(len(particle_db))
+ for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
+ if 'fluka_code' in pData:
+ string += " FLUKACode::{:s}, \n".format(identifier)
+ else:
+ string += " FLUKACode::Unknown, // {:s}\n".format(identifier + ' not implemented in FLUKA')
+ string += "};\n"
+ return string
+
+
+def generate_fluka_canInteract(particle_db):
+ '''
+ generates the look-up table to convert corsika codes to epos codes
+ '''
+ string = "static inline std::vector<bool> const flukaCanInteract = {\n"
+ for identifier, pData in particle_db.items():
+ canInteract = pData.get("fluka_canInteract", False)
+ string += " {:s}, // {:s}\n".format(str(canInteract).lower(), identifier)
+ string += "};\n"
+ return string
+
+
+if __name__ == "__main__":
+ if len(sys.argv) != 3:
+ print("usage: {:s} <particle_db.pkl> <fluka_codes.dat>".format(sys.argv[0]), file=sys.stderr)
+ sys.exit(1)
+
+ print("code_generator.py for FLUKA")
+
+ particle_db = load_particledb(sys.argv[1])
+ read_fluka_codes(sys.argv[2], particle_db)
+
+ with open("Generated.inc", "w") as f:
+ print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
+ print(generate_fluka_enum(particle_db), file=f)
+ print(generate_corsika2fluka(particle_db), file=f)
+ print(generate_fluka_canInteract(particle_db), file=f)
diff --git a/src/modules/fluka/fluka_codes.dat b/src/modules/fluka/fluka_codes.dat
new file mode 100644
index 0000000000000000000000000000000000000000..b7f42677d1b5658672e308c01ceb1afb8346179f
--- /dev/null
+++ b/src/modules/fluka/fluka_codes.dat
@@ -0,0 +1,60 @@
+# FLUKA particle code conversion table
+# C8 identifier, FLUKA code, has cross-section/can interact
+Proton 1 yes
+AntiProton 2 yes
+Electron 3
+Positron 4
+NuE 5
+NuEBar 6
+Photon 7
+Neutron 8 yes
+AntiNeutron 9 yes
+MuPlus 10
+MuMinus 11
+K0Long 12 yes
+PiPlus 13 yes
+PiMinus 14 yes
+KPlus 15 yes
+KMinus 16 yes
+Lambda0 17 yes
+Lambda0Bar 18 yes
+K0Short 19 yes
+SigmaMinus 20 yes
+SigmaPlus 21 yes
+Sigma0 22 yes
+Pi0 23 yes
+K0 24
+K0Bar 25
+NuMu 27
+NuMuBar 28
+SigmaMinusBar 31 yes
+Sigma0Bar 32 yes
+SigmaPlusBar 33 yes
+Xi0 34 yes
+Xi0Bar 35 yes
+XiMinus 36 yes
+XiPlusBar 37 yes
+OmegaMinus 38 yes
+OmegaPlusBar 39 yes
+TauPlus 41
+TauMinus 42
+NuTau 43
+NuTauBar 44
+DPlus 45
+DMinus 46
+D0 47
+D0Bar 48
+DsPlus 49
+DsMinus 50
+LambdaCPlus 51
+XiCPlus 52
+XiC0 53
+XiPrimeCPlus 54
+XiPrimeC0 55
+OmegaC0 56
+LambdaCMinusBar 57
+XiCMinusBar 58
+XiC0Bar 59
+XiPrimeCMinusBar 60
+XiPrimeC0Bar 61
+OmegaC0Bar 62
diff --git a/src/modules/fluka/strip_flukahp.sh b/src/modules/fluka/strip_flukahp.sh
new file mode 100755
index 0000000000000000000000000000000000000000..0ee9693bfcb16c61d90d48929fa1b431031e4127
--- /dev/null
+++ b/src/modules/fluka/strip_flukahp.sh
@@ -0,0 +1,36 @@
+#!/bin/sh
+
+# (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+#
+# See file AUTHORS for a list of contributors.
+#
+# This software is distributed under the terms of the GNU General Public
+# Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+# the license.
+
+# This script strips off flrndm() from the libflukahp.a so that we can provide our own
+# implementation.
+
+flukalibOrig=`realpath $1`
+target="$2"
+
+if [ ! -r "$flukalibOrig" ]; then
+ echo "\"$flukalibOrig\" not readable" 1>&2
+ exit 1
+fi
+
+tmpdir=`mktemp -d fluka_objectsXXXXXX`
+workdir=`pwd`
+
+echo "extracting objects from $flukalibOrig into `realpath $tmpdir`..."
+cd "$tmpdir"
+ar x "$flukalibOrig"
+cd "$workdir"
+rm "$tmpdir/flrndm.o"
+
+[ -f "libflukahp-norndm.a" ] && rm "libflukahp-norndm.a"
+
+echo "creating libflukahp-norndm.a..."
+ar -rcs "$target/libflukahp-norndm.a" "$tmpdir"/*.o
+
+rm -r "$tmpdir"
diff --git a/tests/common/SetupStack.hpp b/tests/common/SetupStack.hpp
index 4a7ceff3ff4590f1bd20a8736aec0ddc4581d704..a598aa30af96517181f75b0afb3e72d6089ac1c6 100644
--- a/tests/common/SetupStack.hpp
+++ b/tests/common/SetupStack.hpp
@@ -24,7 +24,7 @@ namespace corsika::test {
/*
* the version with history
*/
- using Stack = detaill::StackWithHistory;
+ using Stack = detail::StackWithHistory;
#else // WITH_HISTORY
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 015093588f7ec42d0496ba5bb0d323badc425f4d..16a0e66f0d34824f0ed173a5a66a759999dab17f 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -15,7 +15,11 @@ set (test_modules_sources
testRadio.cpp
testEMThinning.cpp
testSophia.cpp
- )
+)
+
+if (EXISTS "${C8_FLUKALIB}")
+ list(APPEND test_modules_sources "testFluka.cpp")
+endif()
CORSIKA_ADD_TEST (testModules SOURCES ${test_modules_sources})
diff --git a/tests/modules/testFluka.cpp b/tests/modules/testFluka.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1978c3f0478c365f682b75b7d0f035dc6745f23b
--- /dev/null
+++ b/tests/modules/testFluka.cpp
@@ -0,0 +1,191 @@
+/*
+ * (c) Copyright 2023 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/modules/FLUKA.hpp>
+#include <corsika/modules/fluka/Random.hpp>
+
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+
+#include <catch2/catch.hpp>
+
+#include <fstream>
+#include <iomanip>
+
+#include <SetupTestStack.hpp>
+#include <SetupTestEnvironment.hpp>
+
+using namespace corsika;
+
+template <typename TStackView>
+auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
+ MomentumVector sum{vCS, 0_eV, 0_eV, 0_eV};
+ for (auto const& p : view) { sum += p.getMomentum(); }
+ return sum;
+}
+
+TEST_CASE("FLUKACodeConversion") {
+ REQUIRE(corsika::fluka::convertToFluka(Code::PiPlus) ==
+ corsika::fluka::FLUKACode::PiPlus);
+ REQUIRE(corsika::fluka::convertToFlukaRaw(Code::PiPlus) == 13);
+ REQUIRE(corsika::fluka::convertToFlukaRaw(Code::Proton) == 1);
+ REQUIRE(corsika::fluka::convertToFlukaRaw(Code::Lambda0) == 17);
+
+ REQUIRE_THROWS(corsika::fluka::convertToFluka(Code::WPlus));
+}
+
+auto setupEnvironment() {
+ using DummyEnvironmentInterface =
+ IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+ using DummyEnvironment = Environment<DummyEnvironmentInterface>;
+ using MyHomogeneousModel = MediumPropertyModel<
+ UniformMagneticField<HomogeneousMedium<DummyEnvironmentInterface>>>;
+ RNGManager<>::getInstance().registerRandomStream("fluka");
+ DummyEnvironment env;
+ auto& universe = *env.getUniverse();
+ CoordinateSystemPtr const& cs = env.getCoordinateSystem();
+ universe.setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, 0_T), 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition{
+ std::vector<Code>{Code::Hydrogen, Code::Oxygen, Code::Nitrogen, Code::Argon},
+ std::vector<double>{.25, .25, .25, .25}});
+
+ return env;
+}
+
+static auto const env = setupEnvironment();
+static corsika::fluka::InteractionModel flukaModel{env};
+static auto const& cs = env.getCoordinateSystem();
+
+TEST_CASE("FLUKA") {
+ //~ auto tup = setupFluka();
+ //~ corsika::fluka::InteractionModel& flukaModel = std::get<0>(tup);
+ //~ auto const env = std::get<1>(tup);
+
+ SECTION("getMaterialIndex") {
+ REQUIRE(flukaModel.getMaterialIndex(Code::Hydrogen) > 0);
+ REQUIRE(flukaModel.getMaterialIndex(Code::Oxygen) > 0);
+ REQUIRE(flukaModel.getMaterialIndex(Code::Nitrogen) > 0);
+ REQUIRE(flukaModel.getMaterialIndex(Code::Argon) > 0);
+
+ // not initialized
+ REQUIRE(flukaModel.getMaterialIndex(Code::Uranium) < 0);
+ }
+
+ SECTION("getCrossSection") {
+ auto const projectileCode =
+ GENERATE(Code::PiMinus, Code::PiMinus, Code::PiMinus, Code::KMinus, Code::K0Long,
+ Code::K0Short, Code::Lambda0, Code::SigmaPlus, Code::Proton,
+ Code::AntiProton, Code::KMinus, Code::K0Long, Code::K0Short,
+ Code::Lambda0, Code::SigmaPlus, Code::Proton, Code::AntiProton);
+
+ auto const targetCode = GENERATE(Code::Oxygen, Code::Hydrogen);
+
+ HEPEnergyType const p = 100_GeV;
+ auto const projectile4mom =
+ FourVector{calculate_total_energy(p, get_mass(projectileCode)),
+ MomentumVector{cs, 0_eV, 0_eV, p}};
+ auto const target4mom =
+ FourVector{get_mass(targetCode), MomentumVector{cs, 0_eV, 0_eV, 0_eV}};
+
+ CHECK(flukaModel.getCrossSection(projectileCode, targetCode, projectile4mom,
+ target4mom) > 0_mb);
+ }
+
+ SECTION("getCrossSection invalid") {
+ auto const projectileCode = GENERATE(Code::Electron, Code::Photon);
+ auto const targetCode = GENERATE(Code::Oxygen, Code::Hydrogen);
+
+ HEPEnergyType const p = 100_GeV;
+ auto const projectile4mom =
+ FourVector{calculate_total_energy(p, get_mass(projectileCode)),
+ MomentumVector{cs, 0_eV, 0_eV, p}};
+ auto const target4mom =
+ FourVector{get_mass(targetCode), MomentumVector{cs, 0_eV, 0_eV, 0_eV}};
+
+ CHECK(flukaModel.getCrossSection(projectileCode, targetCode, projectile4mom,
+ target4mom) == CrossSectionType::zero());
+ }
+
+ SECTION("doInteraction") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ auto const& cs = *csPtr;
+
+ auto const projectileCode = GENERATE(Code::PiPlus, Code::PiMinus, Code::KPlus,
+ Code::K0Long, Code::Lambda0, Code::SigmaPlus);
+ auto const p = GENERATE(1_GeV, 20_GeV, 100_GeV, 1_TeV);
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Hydrogen, 1_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
+ { [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
+
+ auto const targetCode = Code::Oxygen;
+ auto const projectile4mom =
+ FourVector{calculate_total_energy(p, get_mass(projectileCode)),
+ MomentumVector{cs, 0_eV, 0_eV, p}};
+ auto const target4mom =
+ FourVector{get_mass(targetCode), MomentumVector{cs, 0_eV, 0_eV, 0_eV}};
+
+ flukaModel.doInteraction(*secViewPtr, projectileCode, targetCode, projectile4mom,
+ target4mom);
+ auto const pSum = sumMomentum(*secViewPtr, cs);
+
+ CHECK((pSum - projectile4mom.getSpaceLikeComponents()).getNorm() / p ==
+ Approx(0).margin(1e-4));
+ CHECK((pSum.getNorm() - p) / p == Approx(0).margin(1e-4));
+ CHECK(secViewPtr->getSize() > 1);
+ }
+
+ SECTION("doInteraction-invalid-projectile") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ auto const& cs = *csPtr;
+
+ auto const projectileCode = GENERATE(Code::Electron, Code::MuPlus);
+ auto const p = 50_GeV;
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Hydrogen, 1_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
+ { [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
+
+ auto const targetCode = Code::Oxygen;
+ auto const projectile4mom =
+ FourVector{calculate_total_energy(p, get_mass(projectileCode)),
+ MomentumVector{cs, 0_eV, 0_eV, p}};
+ auto const target4mom =
+ FourVector{get_mass(targetCode), MomentumVector{cs, 0_eV, 0_eV, 0_eV}};
+
+ REQUIRE_THROWS(flukaModel.doInteraction(*secViewPtr, projectileCode, targetCode,
+ projectile4mom, target4mom));
+ }
+
+ SECTION("doInteraction-invalid-target") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ auto const& cs = *csPtr;
+
+ auto const projectileCode = Code::AntiNeutron;
+ auto const p = 50_GeV;
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Hydrogen, 1_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
+ { [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
+
+ auto const targetCode = Code::Uranium;
+ auto const projectile4mom =
+ FourVector{calculate_total_energy(p, get_mass(projectileCode)),
+ MomentumVector{cs, 0_eV, 0_eV, p}};
+ auto const target4mom =
+ FourVector{get_mass(targetCode), MomentumVector{cs, 0_eV, 0_eV, 0_eV}};
+
+ REQUIRE_THROWS(flukaModel.doInteraction(*secViewPtr, projectileCode, targetCode,
+ projectile4mom, target4mom));
+ }
+}