From dd14940d325604f6122519074e6212546f64f57c Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 2 Feb 2023 14:56:24 +0100
Subject: [PATCH] first steps
---
.../detail/modules/fluka/InteractionModel.inl | 157 ++++++++++++++++++
corsika/modules/fluka/InteractionModel.hpp | 46 +++++
src/modules/fluka/CMakeLists.txt | 26 +++
src/modules/fluka/code_generator.py | 102 ++++++++++++
src/modules/fluka/fluka_codes.dat | 123 ++++++++++++++
tests/modules/testFluka.cpp | 38 +++++
6 files changed, 492 insertions(+)
create mode 100644 corsika/detail/modules/fluka/InteractionModel.inl
create mode 100644 corsika/modules/fluka/InteractionModel.hpp
create mode 100644 src/modules/fluka/CMakeLists.txt
create mode 100755 src/modules/fluka/code_generator.py
create mode 100644 src/modules/fluka/fluka_codes.dat
create mode 100644 tests/modules/testFluka.cpp
diff --git a/corsika/detail/modules/fluka/InteractionModel.inl b/corsika/detail/modules/fluka/InteractionModel.inl
new file mode 100644
index 000000000..f7ae6fd32
--- /dev/null
+++ b/corsika/detail/modules/fluka/InteractionModel.inl
@@ -0,0 +1,157 @@
+/*
+ * (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 <vector>
+#include <iterator>
+#include <set>
+#include <utility>
+
+#include <boost/iterator/zip_iterator.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <FLUKA.hpp>
+
+namespace corsika::fluka {
+ template <typename TEnvironment>
+ inline InteractionModel::InteractionModel(TEnvironment const& env) :materials_{genFlukaMaterials(env)} {
+ for (auto const& [code, matno]: materials_) {
+ std::cout << code << " " << matno << std::endl;
+ }
+ }
+
+ inline bool InteractionModel::isValid(Code projectileID, Code targetID /*, HEPEnergyType sqrtS*/) {
+ static std::array<Code> constexpr validProjectiles{Code::Proton, Code::AntiProton, Code::Neutron, Code::AntiNeutron,
+ Code::PiPlus, Code::PiMinus, Code::KPlus, Code::KMinus, Code::K0Long, Code::K0Short}
+
+ if (std::find(validProjectiles.cbegin(), validProjectiles.cend(), projectileID) == validProjectiles.cend()) {
+ // invalid projectile
+ return false;
+ }
+
+ if (getMaterialIndex(targetID) == -1) {
+ // unknown material
+ return false;
+ }
+
+ // TODO: check validity range of sqrtS
+ }
+
+ inline int InteractionModel::getMaterialIndex(Code targetID) const {
+ if (auto it = std::find(materials_.cbegin(), materials_.cend(), [=targetID](std::pair<Code, int>& p) {return p.first == targetID;});
+ 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 {
+
+ if (!isValid(projectileId, targetId)) { return CrossSectionType::zero(); }
+
+ COMBoost const targetRestBoost{projecileP4.getSpaceLikeComponents(), get_mass(targetId)};
+ FourMomentum const projectileLab4mom = targetRestBoost.toCOM(projectileP4);
+ projec
+
+ int const iflxyz = 1;
+
+ // SIGREA = SGMXYZ ( KPROJ, MMAT, EKIN, PPROJ, iflxyz );
+ }
+
+ template <typename TSecondaryView>
+ inline void InteractionModel::doInteraction(
+ TSecondaryView& view, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4) {
+
+ }
+
+
+ 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);
+ char crvrck[8+1] = "76466879"; // magic number that FLUKA uses to see if it's the right version
+
+
+ /*
+ * 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
+ */
+ int const iflxyz_ = 1;
+
+ 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);
+
+ // call FLUKA
+ ::fluka::stpxyz_(&nElements, nelmfl.get(), izelfl.data(), wfelml.get(), &mxelfl,
+ &pptmax, &ef2dp3, &df2dp3, &iflxyz_, &lprint, mtflka.get(), crvrck);
+
+ // 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 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;
+
+ }
+}
diff --git a/corsika/modules/fluka/InteractionModel.hpp b/corsika/modules/fluka/InteractionModel.hpp
new file mode 100644
index 000000000..836ec2105
--- /dev/null
+++ b/corsika/modules/fluka/InteractionModel.hpp
@@ -0,0 +1,46 @@
+/*
+ * (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 <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/PhysicalUnits.hpp>
+
+namespace corsika::fluka {
+ class InteractionModel {
+ public:
+ template <typename TEnvironment>
+ InteractionModel(TEnvironment const&);
+
+ CrossSectionType getCrossSection(
+ Code projectileId, Code targetId, FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
+
+ bool isValid(Code projectileID, Code targetID, HEPEnergyType sqrtS) const;
+
+ int getMaterialIndex(Code targetID) const;
+
+ template <typename TSecondaryView>
+ void doInteraction(
+ TSecondaryView& view, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+ //~ static int const iflxyz_ = 1; //!< select interaction types (see fluka.h); hardcoded to inel. for now
+ private:
+ std::vector<std::pair<Code, int>> const materials_; //!< map target Code to FLUKA material no.
+
+ template <typename TEnvironment>
+ static std::vector<std::pair<Code, int>> genFlukaMaterials(TEnvironment const&);
+ // TODO: random number stream
+ };
+}
+
+#include <corsika/detail/modules/fluka/InteractionModel.inl>
diff --git a/src/modules/fluka/CMakeLists.txt b/src/modules/fluka/CMakeLists.txt
new file mode 100644
index 000000000..8ce30ed94
--- /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 000000000..690d1ca32
--- /dev/null
+++ b/src/modules/fluka/code_generator.py
@@ -0,0 +1,102 @@
+#!/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_epos_codes(filename, particle_db):
+ '''
+ reads to epos codes data file
+
+ For particls known to epos, add 'epos_code' and 'epos_xsType' to particle_db
+ '''
+ with open(filename) as f:
+ for line in f:
+ line = line.strip()
+ if len(line) == 0 or line[0].startswith('#'):
+ continue
+ if len(line) == 2:
+ canInteractFlag = "no"
+ identifier, fluka_code = line.split()
+ else:
+ identifier, fluka_code, canInteractFlag = line.split()
+
+ 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"
+ 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 = "std::array<FLUKACode, {:d}> constexpr corsika2epos = {{\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 = "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 EPOS")
+
+ particle_db = load_particledb(sys.argv[1])
+ read_fluka_codes(sys.argv[2], particle_db)
+ # ~ set_default_epos_definition(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 000000000..e2c1c319c
--- /dev/null
+++ b/src/modules/fluka/fluka_codes.dat
@@ -0,0 +1,123 @@
+Photon 7
+Electron 4
+Positron 3
+
+MuPlus 10
+MuMinus 11
+Pi0 23
+PiPlus 13 yes
+PiMinus 14 yes
+K0Long 12 yes
+KPlus 15 yes
+KMinus 16 yes
+Neutron 8 yes
+Proton 1 yes
+AntiProton 2 yes
+K0Short 19 yes
+Eta 0
+Lambda0 17 yes
+SigmaPlus 21
+Sigma0 22
+SigmaMinus 20
+Xi0 34
+XiMinus 36
+OmegaMinus 38
+AntiNeutron 9 yes
+Lambda0Bar 18 yes
+SigmaMinusBar 31
+Sigma0Bar 32
+SigmaPlusBar 33
+Xi0Bar 34
+XiPlusBar 37
+OmegaPlusBar 39
+
+Omega 0
+Rho0 0
+RhoPlus 0
+RhoRhoMinus 0
+DeltaPlusPlus 0
+DeltaPlus 0
+Delta0 0
+DeltaMinus 0
+DeltaMinusMinusBar 0
+DeltaMinusBar 0
+Delta0Bar 0
+DeltaPlusBar 0
+KStar0 0
+KStarPlus 0
+KStarMinusBar 0
+KStar0Bar 0
+NuE 5
+NuEBar 6
+NuMu 27
+NuMubar 28
+
+D0 16
+DPlus 25
+DMinusBar 24
+D0Bar 15
+DSPlus 0
+DSMinusBar 0
+EtaC 0
+DStar0 0
+DStarPlus 0
+DStarMinusBar 0
+DStar0Bar 0
+DStarSPlus 0
+DStarsMinusBar 0
+ 0
+JPsi 0
+TauPlus 41
+TauMinus 42
+NuTau 43
+NuTauBar 44
+ 0
+ 0
+LambdaCPlus 17
+XiCPlus 34
+XiC0 36
+SigmaCPlus 21
+SigmaCPlus 22
+SigmaC0 20
+XiPrimeCPlus 34
+XiPrimeC0 36
+OmegaC0 38
+
+LambdaCMinusBar 18
+XiCMinusBar 35
+XiC0Bar 37
+SigmaCMinusMinusBar 31
+SigmaCMinusBar 32
+SigmaC0Bar 33
+XiPrimeCMinusBar 35
+XiPrimeC0Bar 37
+OmegaC0Bar 39
+
+SigmaStarCPlusPlus 0
+SigmaStarCPlus 0
+SigmaStarC0 0
+
+SigmaStarCBar-- 0
+SigmaStarCMinusBar 0
+SigmaStarC0Bar 0
+
+B0 0
+BPlus 0
+BMinusBar 0
+B0Bar 0
+BS0 0
+BS0Bar 0
+BCPlus 0
+BCMinusBar 0
+LambdaB0 0
+SigmaBMinus 0
+SigmaBPlus 0
+XiB0 0
+XiBMinus 0
+OmegaBMinus 0
+LambdaB0Bar 0
+SigmaBPlusBar 0
+SigmaBMinusBar 0
+XiB0Bar 0
+XiBPlusBar 0
+OmegaBPlusBar 0
diff --git a/tests/modules/testFluka.cpp b/tests/modules/testFluka.cpp
new file mode 100644
index 000000000..96f265bf6
--- /dev/null
+++ b/tests/modules/testFluka.cpp
@@ -0,0 +1,38 @@
+/*
+ * (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/InteractionModel.hpp>
+//~ #include <SetupTestEnvironment.hpp>
+
+#include <corsika/framework/geometry/Point.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>
+
+using namespace corsika;
+
+TEST_CASE("FLUKA") {
+ using DummyEnvironmentInterface = IMediumModel;
+ using DummyEnvironment = Environment<DummyEnvironmentInterface>;
+ using MyHomogeneousModel = HomogeneousMedium<DummyEnvironmentInterface>;
+
+ DummyEnvironment env;
+ auto& universe = *env.getUniverse();
+ CoordinateSystemPtr const& cs = env.getCoordinateSystem();
+ universe.setModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Hydrogen, Code::Oxygen, Code::Nitrogen, Code::Argon},
+ std::vector<double>{1./4., 1./4., 1./4., 1./4.}));
+
+ corsika::fluka::InteractionModel flukaModel{env};
+
+}
--
GitLab