IAP GITLAB

Skip to content
Snippets Groups Projects
Commit c5ea3b71 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

linking problems fixed

parent 3b40b8cc
No related branches found
No related tags found
1 merge request!468Resolve "Add FLUKA"
Hide whitespace changes
Inline Side-by-side
corsika/detail/modules/fluka/InteractionModel.inl
+ 5
− 10
View file @ c5ea3b71
......@@ -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);
......
corsika/modules/FLUKA.hpp
+ 6
− 2
View file @ c5ea3b71
......@@ -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
modules/fluka/CMakeLists.txt
+ 12
− 9
View file @ c5ea3b71
......@@ -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()
modules/fluka/FLUKA.hpp
+ 35
− 34
View file @ c5ea3b71
......@@ -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
modules/fluka/flrndm.cpp → modules/fluka/fluka.cpp
+ 32
− 0
View file @ c5ea3b71
......@@ -9,23 +9,24 @@
#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();
/**
* 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
tests/modules/testFluka.cpp
+ 60
− 73
View file @ c5ea3b71
......@@ -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);
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum