contd.

corsika/detail/modules/fluka/InteractionModel.inl

@@ -23,6 +23,7 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
 #include <FLUKA.hpp>
+#include <ParticleConversion.hpp>
 
 namespace corsika::fluka {
   template <typename TEnvironment>
@@ -32,11 +33,8 @@ namespace corsika::fluka {
       }
     }
     
-  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()) {
+  inline bool InteractionModel::isValid(Code projectileID, Code targetID, HEPEnergyType /*sqrtS*/) const {
+      if (fluka::canInteract(projectileID)) {
           // invalid projectile
         return false;
       }
@@ -47,10 +45,11 @@ namespace corsika::fluka {
       }
         
       // TODO: check validity range of sqrtS
+      return true;
   }
   
   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;});
+    if (auto it = std::find(materials_.cbegin(), materials_.cend(), [=](std::pair<Code, int> const& p) {return p.first == targetID;});
         it == materials_.cend()) {
         return -1;
     } else {
@@ -61,12 +60,12 @@ namespace corsika::fluka {
   inline CrossSectionType InteractionModel::getCrossSection(
       Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
       FourMomentum const& targetP4) const {
-          
-    if (!isValid(projectileId, targetId)) { return CrossSectionType::zero(); }
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+    if (!isValid(projectileId, targetId, sqrtS)) { return CrossSectionType::zero(); }
     
-    COMBoost const targetRestBoost{projecileP4.getSpaceLikeComponents(), get_mass(targetId)};
-    FourMomentum const projectileLab4mom = targetRestBoost.toCOM(projectileP4);
-    projec
+    COMBoost const targetRestBoost{projectileP4.getSpaceLikeComponents(), get_mass(targetId)};
+    FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+    //~ projec
     
     int const iflxyz = 1;
     

modules/fluka/ParticleConversion.hpp

+/*
+ * (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) {
+    return flukaCanInteract[static_cast<CodeIntType>(code)];
+  }
+} // namespace corsika::fluka

src/CMakeLists.txt

-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)

src/modules/fluka/code_generator.py

@@ -21,22 +21,23 @@ def load_particledb(filename):
     return particle_db
 
 
-def read_epos_codes(filename, particle_db):
+def read_fluka_codes(filename, particle_db):
     '''
-    reads to epos codes data file
+    reads to FLUKA codes data file
 
-    For particls known to epos, add 'epos_code' and 'epos_xsType' to particle_db
+    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.split()
+                identifier, fluka_code = line
             else:
-                identifier, fluka_code, canInteractFlag = line.split()
+                identifier, fluka_code, canInteractFlag = line
 
             try:
                 particle_db[identifier]["fluka_code"] = int(fluka_code)
@@ -49,7 +50,8 @@ def generate_fluka_enum(particle_db):
     '''
      generates the enum to access epos particles by readable names
     '''
-    output = "enum class FLUKACode : int {\n"
+    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'])
@@ -61,7 +63,7 @@ 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))
+    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:
@@ -76,10 +78,10 @@ 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"
+    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 += "  {:s}, // {:s}\n".format(str(canInteract).lower(), identifier)
     string += "};\n"
     return string
 
@@ -89,11 +91,10 @@ if __name__ == "__main__":
         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")
+    print("code_generator.py for FLUKA")
     
     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)

src/modules/fluka/fluka_codes.dat

@@ -14,7 +14,7 @@ Neutron	8	yes
 Proton	1	yes
 AntiProton	2	yes
 K0Short	19	yes
-Eta	0
+# Eta	0
 Lambda0	17	yes
 SigmaPlus	21
 Sigma0	22
@@ -31,48 +31,48 @@ 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
+# Omega	0
+# Rho0	0
+# RhoPlus	0
+# RhoMinus	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
+NuMuBar	28
 
 D0	16
 DPlus	25
-DMinusBar	24
+DMinus	24
 D0Bar	15
-DSPlus	0
-DSMinusBar	0
-EtaC	0
-DStar0	0
-DStarPlus	0
-DStarMinusBar	0
-DStar0Bar	0
-DStarSPlus	0
-DStarsMinusBar	0
-	0
-JPsi	0
+# 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
+#	0
+#	0
 LambdaCPlus	17
 XiCPlus	34
 XiC0	36
@@ -93,31 +93,31 @@ XiPrimeCMinusBar	35
 XiPrimeC0Bar	37
 OmegaC0Bar	39
 
-SigmaStarCPlusPlus	0
-SigmaStarCPlus	0
-SigmaStarC0	0
+# SigmaStarCPlusPlus	0
+# SigmaStarCPlus	0
+# SigmaStarC0	0
 
-SigmaStarCBar--	0
-SigmaStarCMinusBar	0
-SigmaStarC0Bar	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
+# 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