diff --git a/corsika/detail/modules/fluka/InteractionModel.inl b/corsika/detail/modules/fluka/InteractionModel.inl
index 627e64ba31a0b784a6aa7cc2068c12277d2f802f..8983d3980dc15655990ecb866b44655e4d395521 100644
--- a/corsika/detail/modules/fluka/InteractionModel.inl
+++ b/corsika/detail/modules/fluka/InteractionModel.inl
@@ -130,7 +130,7 @@ namespace corsika::fluka {
auto const plab = projectileLab4mom.getSpaceLikeComponents();
auto const& cs = plab.getCoordinateSystem();
auto const labMomentum = plab.getNorm();
- double const labMomentumGeV = labMomentum * invGeV;
+ double const labMomentumGeV = labMomentum * invGeV * 0;
auto const direction = (plab / labMomentum).getComponents().getEigenVector();
@@ -141,9 +141,8 @@ namespace corsika::fluka {
for (int i = 0; i < ::fluka::hepevt_.nhep; ++i) {
int const status = ::fluka::hepevt_.isthep[i];
- // TODO: enable this when FLUKA writes status code
- //~ if (status != 1) // skip non-final-state particles
- //~ continue;
+ 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);
@@ -152,7 +151,6 @@ namespace corsika::fluka {
(1_GeV).magnitude()};
auto const pPrime = mom.getNorm();
auto const c8mass = corsika::get_mass(c8code);
- auto const flMass = ::fluka::hepevt_.phep[i][5 - 1];
auto const fourMomCollisionFrame =
FourVector{calculate_total_energy(pPrime, c8mass), MomentumVector{cs, mom}};
@@ -162,9 +160,6 @@ namespace corsika::fluka {
view.addSecondary(std::tuple{c8code, corsika::calculate_kinetic_energy(p, c8mass),
momOrigFrame / p});
- std::cout << static_cast<int>(pdg) << '\t' << get_name(c8code, full_name{}) << '\t'
- << momOrigFrame / 1_GeV << '\t' << flMass << ' ' << c8mass / 1_GeV << " "
- << std::endl;
}
}
@@ -172,8 +167,8 @@ namespace corsika::fluka {
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
+ // generate complete list of all nuclei types in universe
auto const allElementsInUniverse = std::invoke([&]() {
std::set<Code> allElementsInUniverse;
auto collectElements = [&](auto& vtn) {
@@ -204,7 +199,7 @@ namespace corsika::fluka {
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 crvrck[] = "76466879";
+ char const crvrck[] = "76466879";
int const size = 8;
std::fill(&nelmfl[0], &nelmfl[nElements], 1);
diff --git a/corsika/modules/FLUKA.hpp b/corsika/modules/FLUKA.hpp
index 2b75999cff0d1ff66cda97350d1d2f55db221a4a..43f4e32e30ebf2d05d7f2edaca9b9c6b41bdf700 100644
--- a/corsika/modules/FLUKA.hpp
+++ b/corsika/modules/FLUKA.hpp
@@ -8,7 +8,6 @@
#pragma once
-#include <corsika/modules/fluka/ParticleConversion.hpp>
#include <corsika/modules/fluka/InteractionModel.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
@@ -27,5 +26,10 @@ namespace corsika::fluka {
* to provide all the functions for ProcessSequence.
*/
class Interaction : public fluka::InteractionModel,
- public corsika::InteractionProcess<Interaction> {};
+ public corsika::InteractionProcess<Interaction> {
+ public:
+ template <typename TEnvironment>
+ Interaction(TEnvironment const& env)
+ : fluka::InteractionModel{env} {}
+ };
} // namespace corsika::fluka
diff --git a/modules/fluka/CMakeLists.txt b/modules/fluka/CMakeLists.txt
index 2d56bdc85789e1eef9d1f611e80baaa60b9ce875..8c9c17fb44805efad2a301bc822c34a0a02632a8 100644
--- a/modules/fluka/CMakeLists.txt
+++ b/modules/fluka/CMakeLists.txt
@@ -17,17 +17,20 @@ if (DEFINED C8_FLUKALIB)
$<INSTALL_INTERFACE:include/corsika_modules/fluka>
)
- # this contains our own implementation of flrndm_()
- add_library(flukaRndm SHARED flrndm.cpp)
- target_link_libraries(flukaRndm PRIVATE flukahp-norndm)
+ # 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 (flukaRndm PUBLIC
+ target_include_directories (fluka PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
$<INSTALL_INTERFACE:include/corsika_modules/fluka>
)
set_target_properties (
- flukaRndm
+ fluka
PROPERTIES
POSITION_INDEPENDENT_CODE 1
)
@@ -35,17 +38,17 @@ if (DEFINED C8_FLUKALIB)
install (
FILES
fluka.hpp
- DESTINATION include/corsika_modules/urqmd
+ DESTINATION include/corsika_modules/fluka
)
install (
- TARGETS flukaRndm
+ TARGETS fluka
EXPORT CORSIKA8PublicTargets
ARCHIVE DESTINATION lib/corsika
INCLUDES DESTINATION include/corsika_modules/fluka
)
# add fluka to corsika8 build
- add_dependencies (CORSIKA8 flukaRndm)
- target_link_libraries (CORSIKA8 INTERFACE flukaRndm)
+ add_dependencies (CORSIKA8 fluka)
+ target_link_libraries (CORSIKA8 INTERFACE fluka)
endif()
diff --git a/modules/fluka/FLUKA.hpp b/modules/fluka/FLUKA.hpp
index 25b55f9042ef60ee59c502adcab79fa6e421e9ac..7c5ab5f168fba1056cb086877fb4fd733bba1cf5 100644
--- a/modules/fluka/FLUKA.hpp
+++ b/modules/fluka/FLUKA.hpp
@@ -26,34 +26,34 @@ namespace fluka {
using hepmc_array = std::array<T, nmxhep>;
extern "C" {
- /*
- c---------------------------------------------------------------------------
- c hep standard event commonblock.
- c---------------------------------------------------------------------------
- c---------------------------------------------------------------------------
- c
- c nevhep - event number
- c nhep - number of entries in the event record
- c
- c isthep(i) - status code
- c idhep(i) - particle id (particle data group standard)
- c
- c jmohep(1,i) - position of mother particle in list
- c jmohep(2,i) - position of second mother particle in list
- c jdahep(1,i) - position of first daughter in list
- c jdahep(2,i) - position of first daughter in list
- c
- c phep(1,i) - p_x momentum in gev/c
- c phep(2,i) - p_y momentum in gev/c
- c phep(3,i) - p_z momentum in gev/c
- c phep(4,i) - energy in gev
- c phep(5,i) - mass in gev/c**2
- c
- c vhep(1,i) - x position of production vertex in mm
- c vhep(2,i) - y position of production vertex in mm
- c vhep(3,i) - z position of production vertex in mm
- c vhep(4,i) - time of production in mm/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
@@ -76,7 +76,7 @@ namespace fluka {
* Iflxyz =111 -> inelastic + elastic + emd
*/
- /*----------------------------------------------------------------------*
+ /**---------------------------------------------------------------------*
* *
* Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
* All Rights Reserved. *
@@ -130,7 +130,7 @@ namespace fluka {
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. *
@@ -165,7 +165,7 @@ namespace fluka {
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. *
@@ -185,7 +185,7 @@ namespace fluka {
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. *
@@ -203,7 +203,8 @@ namespace fluka {
* *
*----------------------------------------------------------------------*/
void fllhep_();
-
- double flrndm();
+
+ //! random-number generator called from within FLUKA
+ double flrndm_();
}
} // namespace fluka
diff --git a/modules/fluka/flrndm.cpp b/modules/fluka/flrndm.cpp
deleted file mode 100644
index 79f3a82aef6aef9ffb9123e47de158dc493b6257..0000000000000000000000000000000000000000
--- a/modules/fluka/flrndm.cpp
+++ /dev/null
@@ -1,31 +0,0 @@
-/*
- * (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 {
-
- /**
- * The following (function) pointers make sure the corresponding objects in libflukahp.a
- * won't get dropped from the final file during linking.
- */
-
- auto* const hepevt_ptr = &hepevt_;
- auto* const stpxyc_ptr = &stpxyz_;
- auto* const evtxyz_ptr = &evtxyz_;
- auto* const sgmxyz_ptr = &sgmxyz_;
- auto* const fllhep_ptr = &fllhep_;
-
-extern "C" {
- double flrndm_() { return ::fluka::rndm_interface(); }
- void bla() {
- fllhep_ptr();
- }
-}
-
-}
diff --git a/modules/fluka/fluka.cpp b/modules/fluka/fluka.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1843eab1bd7618f2622797157ae3ab7336435b62
--- /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 {
+
+ /**
+ * The following (function) pointers make sure the corresponding objects in libflukahp.a
+ * won't get dropped from the final file during linking.
+ */
+
+ auto* const hepevt_ptr = &hepevt_;
+ auto* const stpxyc_ptr = &stpxyz_;
+ auto* const evtxyz_ptr = &evtxyz_;
+ auto* const sgmxyz_ptr = &sgmxyz_;
+
+ 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; }
+ }
+
+} // namespace fluka
diff --git a/tests/modules/testFluka.cpp b/tests/modules/testFluka.cpp
index 6ddcd96c8ba216c866f3def14cc911a9ae59fca1..573332bb9712b097f5ede9ac1c6c771fd4fa218d 100644
--- a/tests/modules/testFluka.cpp
+++ b/tests/modules/testFluka.cpp
@@ -8,7 +8,6 @@
#include <corsika/modules/FLUKA.hpp>
#include <corsika/modules/fluka/Random.hpp>
-//~ #include <SetupTestEnvironment.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
@@ -45,15 +44,13 @@ TEST_CASE("FLUKACodeConversion") {
REQUIRE(corsika::fluka::convertToFlukaRaw(Code::Lambda0) == 17);
}
-TEST_CASE("FLUKA") {
+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();
@@ -63,81 +60,71 @@ TEST_CASE("FLUKA") {
std::vector<Code>{Code::Hydrogen, Code::Oxygen, Code::Nitrogen, Code::Argon},
std::vector<double>{.25, .25, .25, .25}});
- corsika::fluka::InteractionModel flukaModel{env};
-
- // test 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);
- REQUIRE(flukaModel.getMaterialIndex(Code::Uranium) < 0);
-
- { // test getCrossSection for allowed projectile/target combinations
- auto combinationsOK = std::vector{
- std::tuple{Code::PiMinus, Code::Hydrogen},
- std::tuple{Code::PiMinus, Code::Nitrogen},
- std::tuple{Code::PiMinus, Code::Oxygen},
- std::tuple{Code::KMinus, Code::Oxygen},
- std::tuple{Code::K0Long, Code::Oxygen},
- std::tuple{Code::K0Short, Code::Oxygen},
- std::tuple{Code::Lambda0, Code::Oxygen},
- std::tuple{Code::SigmaPlus, Code::Oxygen},
- std::tuple{Code::Proton, Code::Oxygen},
- std::tuple{Code::AntiProton, Code::Oxygen},
- std::tuple{Code::KMinus, Code::Hydrogen},
- std::tuple{Code::K0Long, Code::Hydrogen},
- std::tuple{Code::K0Short, Code::Hydrogen},
- std::tuple{Code::Lambda0, Code::Hydrogen},
- std::tuple{Code::SigmaPlus, Code::Hydrogen},
- std::tuple{Code::Proton, Code::Hydrogen},
- std::tuple{Code::AntiProton, Code::Hydrogen},
- };
+ 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);
+ 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;
- for (auto const& [projectileCode, targetCode] : combinationsOK) {
- 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);
- }
+ 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("doInteraction") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
auto const& cs = *csPtr;
- auto const projectiles = std::array{Code::PiPlus, Code::PiMinus, Code::KPlus,
- Code::K0Long, Code::Lambda0, Code::SigmaPlus};
- auto const momenta = std::array{100_MeV, 1_GeV, 20_GeV, 100_GeV, 1_TeV};
-
- for (auto const p : momenta) {
- for (auto const projectileCode : projectiles) {
- std::cout << "===== " << projectileCode << " @ " << p << "\n";
-
- 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);
- std::cout << "psum = " << pSum << '\n';
- CHECK((pSum - projectile4mom.getSpaceLikeComponents()).getNorm() / p ==
- Approx(0).margin(1e-4));
- CHECK((pSum.getNorm() - p) / p == Approx(0).margin(1e-4));
- }
- }
+ 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);
}
}