first steps

corsika/detail/modules/fluka/InteractionModel.inl

+/*
+ * (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;
+
+    }
+}

corsika/modules/fluka/InteractionModel.hpp

+/*
+ * (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>

src/modules/fluka/CMakeLists.txt

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

src/modules/fluka/code_generator.py

+#!/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)

src/modules/fluka/fluka_codes.dat

+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

tests/modules/testFluka.cpp

+/*
+ * (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};
+    
+}