From 3b40b8cc9800cafe57e99774fac6ad97d4a78c39 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 9 Mar 2023 16:20:25 +0100
Subject: [PATCH] event generation, test, linking
---
.../detail/modules/fluka/InteractionModel.inl | 50 ++++++++++-------
corsika/modules/FLUKA.hpp | 4 +-
corsika/modules/Random.hpp | 1 +
corsika/modules/fluka/InteractionModel.hpp | 5 +-
corsika/modules/fluka/ParticleConversion.hpp | 9 +--
corsika/modules/fluka/Random.hpp | 26 +++++++++
modules/fluka/CMakeLists.txt | 47 ++++++++++++++--
modules/fluka/FLUKA.hpp | 16 +++++-
modules/fluka/flrndm.cpp | 31 +++++++++++
src/modules/fluka/CMakeLists.txt | 2 +-
src/modules/fluka/strip_flukahp.sh | 28 ++++++++++
tests/modules/testFluka.cpp | 55 ++++++++++++++-----
12 files changed, 223 insertions(+), 51 deletions(-)
create mode 100644 corsika/modules/fluka/Random.hpp
create mode 100644 modules/fluka/flrndm.cpp
create mode 100755 src/modules/fluka/strip_flukahp.sh
diff --git a/corsika/detail/modules/fluka/InteractionModel.inl b/corsika/detail/modules/fluka/InteractionModel.inl
index 2798d9fb5..627e64ba3 100644
--- a/corsika/detail/modules/fluka/InteractionModel.inl
+++ b/corsika/detail/modules/fluka/InteractionModel.inl
@@ -21,6 +21,8 @@
#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>
@@ -90,8 +92,6 @@ namespace corsika::fluka {
auto const plab = projectileLab4mom.getSpaceLikeComponents();
- std::cout << targetRestBoost.boost_ << '\n';
- std::cout << targetRestBoost.rotatedCS_->getTransform().matrix() << '\n';
CORSIKA_LOGGER_DEBUG(logger_, fmt::format("Elab = {} GeV", Elab * invGeV));
CORSIKA_LOGGER_DEBUG(logger_, fmt::format("EkinLab = {} GeV", EkinLab * invGeV));
@@ -139,24 +139,32 @@ namespace corsika::fluka {
cumsgx_.get() + materials_.size(),
cumsgx_.get() + materials_.size() * 2);
- //~ 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_;
-
for (int i = 0; i < ::fluka::hepevt_.nhep; ++i) {
- int const pdg = ::fluka::hepevt_.idhep[i];
int const status = ::fluka::hepevt_.isthep[i];
- auto const mom = QuantityVector<hepenergy_d>{
- Eigen::Map<Eigen::Vector3d>(&::fluka::hepevt_.phep[i][0]) * invGeV.magnitude()};
+ // TODO: enable this when FLUKA writes status code
+ //~ if (status != 1) // skip non-final-state particles
+ //~ continue;
- std::cout << pdg << '\t' << status << '\t' << mom << std::endl;
+ 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 flMass = ::fluka::hepevt_.phep[i][5 - 1];
+
+ 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});
+ std::cout << static_cast<int>(pdg) << '\t' << get_name(c8code, full_name{}) << '\t'
+ << momOrigFrame / 1_GeV << '\t' << flMass << ' ' << c8mass / 1_GeV << " "
+ << std::endl;
}
}
@@ -195,8 +203,9 @@ namespace corsika::fluka {
double const df2dp3 = -1; // default
bool const lprint = true;
auto mtflka = std::make_unique<int[]>(mxelfl);
- char crvrck[8 + 1] =
- "76466879"; // magic number that FLUKA uses to see if it's the right version
+ // magic number that FLUKA uses to see if it's the right version
+ char crvrck[] = "76466879";
+ int const size = 8;
std::fill(&nelmfl[0], &nelmfl[nElements], 1);
std::fill(&wfelml[0], &wfelml[nElements], 1.);
@@ -210,7 +219,8 @@ namespace corsika::fluka {
// call FLUKA
::fluka::stpxyz_(&nElements, nelmfl.get(), izelfl.data(), wfelml.get(), &mxelfl,
- &pptmax, &ef2dp3, &df2dp3, &iflxyz_, &lprint, mtflka.get(), crvrck);
+ &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;
diff --git a/corsika/modules/FLUKA.hpp b/corsika/modules/FLUKA.hpp
index 02555ea5e..2b75999cf 100644
--- a/corsika/modules/FLUKA.hpp
+++ b/corsika/modules/FLUKA.hpp
@@ -13,9 +13,9 @@
#include <corsika/framework/process/InteractionProcess.hpp>
/**
- * @file Epos.hpp
+ * @file FLUKA.hpp
*
- * Includes all the parts of the EPOS model. Defines the InteractionProcess<TModel>
+ * Includes all the parts of the FLUKA model. Defines the InteractionProcess<TModel>
* classes needed for the ProcessSequence.
*/
diff --git a/corsika/modules/Random.hpp b/corsika/modules/Random.hpp
index a81bc1ad2..60d8b9a9a 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/fluka/InteractionModel.hpp b/corsika/modules/fluka/InteractionModel.hpp
index f662040bc..206319f85 100644
--- a/corsika/modules/fluka/InteractionModel.hpp
+++ b/corsika/modules/fluka/InteractionModel.hpp
@@ -18,6 +18,7 @@
#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 {
class InteractionModel {
@@ -39,14 +40,14 @@ namespace corsika::fluka {
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*
+ 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&);
- // TODO: random number stream
};
inline static int const iflxyz_ = 1;
diff --git a/corsika/modules/fluka/ParticleConversion.hpp b/corsika/modules/fluka/ParticleConversion.hpp
index 64e97dea0..770cc8c53 100644
--- a/corsika/modules/fluka/ParticleConversion.hpp
+++ b/corsika/modules/fluka/ParticleConversion.hpp
@@ -27,13 +27,14 @@ namespace corsika::fluka {
FLUKACode constexpr convertToFluka(Code const c8id) {
if (is_nucleus(c8id)) {
- throw std::runtime_error{"nucleus conversion to FLUKA not implemented"};
- }
+ 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()};
- }
+ throw std::runtime_error{
+ fmt::format("no correspondig FLUKA id for {}", get_name(c8id)).c_str()};
+ }
return flukaID;
}
diff --git a/corsika/modules/fluka/Random.hpp b/corsika/modules/fluka/Random.hpp
new file mode 100644
index 000000000..73cf5377a
--- /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
index 626901998..2d56bdc85 100644
--- a/modules/fluka/CMakeLists.txt
+++ b/modules/fluka/CMakeLists.txt
@@ -1,16 +1,51 @@
set(C8_FLUKALIB CACHE STRING "path to libflukahp.a")
if (DEFINED C8_FLUKALIB)
- add_library(flukahp STATIC IMPORTED)
- set_property(TARGET flukahp PROPERTY IMPORTED_LOCATION ${C8_FLUKALIB})
+ 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_BINARY_DIR}/libflukahp-norndm.a COMMAND ${input_dir}/strip_flukahp.sh ${C8_FLUKALIB} ${CMAKE_BINARY_DIR})
+ add_custom_target(generate_libfluka-norndm DEPENDS ${CMAKE_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_BINARY_DIR}/libflukahp-norndm.a)
find_library(MATH_LIBRARY m)
- target_link_libraries(flukahp INTERFACE gfortran ${MATH_LIBRARY})
- target_include_directories (flukahp INTERFACE
+ 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>
)
+ # this contains our own implementation of flrndm_()
+ add_library(flukaRndm SHARED flrndm.cpp)
+ target_link_libraries(flukaRndm PRIVATE flukahp-norndm)
+
+ target_include_directories (flukaRndm PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+ $<INSTALL_INTERFACE:include/corsika_modules/fluka>
+ )
+
+ set_target_properties (
+ flukaRndm
+ PROPERTIES
+ POSITION_INDEPENDENT_CODE 1
+ )
+
+ install (
+ FILES
+ fluka.hpp
+ DESTINATION include/corsika_modules/urqmd
+ )
+
+ install (
+ TARGETS flukaRndm
+ EXPORT CORSIKA8PublicTargets
+ ARCHIVE DESTINATION lib/corsika
+ INCLUDES DESTINATION include/corsika_modules/fluka
+ )
+
# add fluka to corsika8 build
- add_dependencies (CORSIKA8 flukahp)
- target_link_libraries (CORSIKA8 INTERFACE flukahp)
+ add_dependencies (CORSIKA8 flukaRndm)
+ target_link_libraries (CORSIKA8 INTERFACE flukaRndm)
endif()
diff --git a/modules/fluka/FLUKA.hpp b/modules/fluka/FLUKA.hpp
index e24f6d22e..25b55f904 100644
--- a/modules/fluka/FLUKA.hpp
+++ b/modules/fluka/FLUKA.hpp
@@ -8,7 +8,19 @@
#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>;
@@ -116,7 +128,7 @@ namespace fluka {
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 CRVRCK[8]);
+ bool const* LPRINT, int* MTFLKA, char const* CRVRCK, int const*);
/*----------------------------------------------------------------------*
* *
@@ -191,5 +203,7 @@ namespace fluka {
* *
*----------------------------------------------------------------------*/
void fllhep_();
+
+ double flrndm();
}
} // namespace fluka
diff --git a/modules/fluka/flrndm.cpp b/modules/fluka/flrndm.cpp
new file mode 100644
index 000000000..79f3a82ae
--- /dev/null
+++ b/modules/fluka/flrndm.cpp
@@ -0,0 +1,31 @@
+/*
+ * (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/src/modules/fluka/CMakeLists.txt b/src/modules/fluka/CMakeLists.txt
index 8ce30ed94..2d2d26b44 100644
--- a/src/modules/fluka/CMakeLists.txt
+++ b/src/modules/fluka/CMakeLists.txt
@@ -23,4 +23,4 @@ add_dependencies (CORSIKA8 SourceDirLinkFLUKA)
install (
FILES ${output_dir}/Generated.inc
DESTINATION include/corsika/modules/fluka
- )
+)
diff --git a/src/modules/fluka/strip_flukahp.sh b/src/modules/fluka/strip_flukahp.sh
new file mode 100755
index 000000000..0599ca987
--- /dev/null
+++ b/src/modules/fluka/strip_flukahp.sh
@@ -0,0 +1,28 @@
+#!/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="$1"
+target="$2"
+
+tmpdir=`mktemp -d fluka_objectsXXXXXX`
+
+echo "extracting objects from $1 into `realpath $tmpdir`..."
+ar --output "$tmpdir" x "$flukalibOrig"
+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/modules/testFluka.cpp b/tests/modules/testFluka.cpp
index 0770bf06c..6ddcd96c8 100644
--- a/tests/modules/testFluka.cpp
+++ b/tests/modules/testFluka.cpp
@@ -7,6 +7,7 @@
*/
#include <corsika/modules/FLUKA.hpp>
+#include <corsika/modules/fluka/Random.hpp>
//~ #include <SetupTestEnvironment.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
@@ -29,6 +30,13 @@
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);
@@ -44,6 +52,8 @@ TEST_CASE("FLUKA") {
using MyHomogeneousModel = MediumPropertyModel<
UniformMagneticField<HomogeneousMedium<DummyEnvironmentInterface>>>;
+ RNGManager<>::getInstance().registerRandomStream("fluka");
+
DummyEnvironment env;
auto& universe = *env.getUniverse();
CoordinateSystemPtr const& cs = env.getCoordinateSystem();
@@ -99,20 +109,35 @@ TEST_CASE("FLUKA") {
{
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
auto const& cs = *csPtr;
- 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 projectileCode = Code::PiPlus;
- auto const targetCode = Code::Nitrogen;
- auto const p = 20_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}};
-
- flukaModel.doInteraction(*secViewPtr, projectileCode, targetCode, projectile4mom,
- target4mom);
+
+ 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));
+ }
+ }
}
}
--
GitLab