lots of things corrected, cross-sections working now

corsika/detail/modules/fluka/InteractionModel.inl

@@ -5,10 +5,11 @@
  * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
  * the license.
  */
- 
+
 #pragma once
 
 #include <algorithm>
+#include <cstdlib>
 #include <vector>
 #include <iterator>
 #include <set>
@@ -23,66 +24,94 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
 #include <FLUKA.hpp>
-#include <ParticleConversion.hpp>
+#include <corsika/modules/fluka/ParticleConversion.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 InteractionModel::InteractionModel(TEnvironment const& env)
+      : materials_{genFlukaMaterials(env)} {
+    for (auto const& [code, matno] : materials_) {
+      std::cout << get_name(code, full_name{}) << "    " << matno << std::endl;
+    }
+  }
+
+  inline bool InteractionModel::isValid(Code projectileID, int material,
+                                        HEPEnergyType /*sqrtS*/) const {
+    if (!fluka::canInteract(projectileID)) {
+      // invalid projectile
+      //~ std::cout << "invalid projecile (cannot interact)";
+      return false;
     }
-    
-  inline bool InteractionModel::isValid(Code projectileID, Code targetID, HEPEnergyType /*sqrtS*/) const {
-      if (fluka::canInteract(projectileID)) {
-          // invalid projectile
-        return false;
-      }
-      
-      if (getMaterialIndex(targetID) == -1) {
-        // unknown material
-        return false;
-      }
-        
-      // TODO: check validity range of sqrtS
-      return true;
+
+    if (material < 0) {
+      // invalid/uninitialized target
+      return false;
+    }
+
+    // TODO: check validity range of sqrtS
+    return true;
+  }
+
+  inline bool InteractionModel::isValid(Code projectileID, Code targetID,
+                                        HEPEnergyType sqrtS) const {
+    return isValid(projectileID, getMaterialIndex(targetID), sqrtS);
   }
-  
+
   inline int InteractionModel::getMaterialIndex(Code targetID) const {
-    if (auto it = std::find(materials_.cbegin(), materials_.cend(), [=](std::pair<Code, int> const& p) {return p.first == targetID;});
+    auto compare = [=](std::pair<Code, int> const& p) { return p.first == targetID; };
+    if (auto it = std::find_if(materials_.cbegin(), materials_.cend(), compare);
         it == materials_.cend()) {
-        return -1;
+      return -1;
     } else {
-        return it->second;
+      return it->second;
     }
   }
-        
+
   inline CrossSectionType InteractionModel::getCrossSection(
       Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
       FourMomentum const& targetP4) const {
+    auto const flukaCodeProj =
+        static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+    auto const flukaMaterial = getMaterialIndex(targetId);
+
     HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
-    if (!isValid(projectileId, targetId, sqrtS)) { return CrossSectionType::zero(); }
-    
-    COMBoost const targetRestBoost{projectileP4.getSpaceLikeComponents(), get_mass(targetId)};
+    if (!isValid(projectileId, flukaMaterial, sqrtS)) { return CrossSectionType::zero(); }
+
+    COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
     FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
-    //~ projec
-    
-    int const iflxyz = 1;
-    
-    // SIGREA = SGMXYZ ( KPROJ, MMAT, EKIN, PPROJ, iflxyz );
+    HEPEnergyType const Elab = projectileLab4mom.getTimeLikeComponent();
+    auto constexpr invGeV = 1 / 1_GeV;
+    double const EkinLab = (Elab - get_mass(projectileId)) * invGeV;
+    double const labMomentum =
+        projectileLab4mom.getSpaceLikeComponents().getNorm() * invGeV;
+
+    auto const plab = projectileLab4mom.getSpaceLikeComponents();
+
+    std::cout << targetRestBoost.boost_ << '\n';
+    std::cout << targetRestBoost.rotatedCS_->getTransform().matrix() << '\n';
+    std::cout << "Elab / GeV = " << Elab * invGeV << '\n';
+    std::cout << "plab = " << plab << '\n';
+    std::cout << "EkinLab / GeV = " << EkinLab << '\n';
+    std::cout << "labMomentum / GeV = "
+              << projectileLab4mom.getSpaceLikeComponents() * invGeV << '\n';
+
+    CrossSectionType const xs = ::fluka::sgmxyz_(&flukaCodeProj, &flukaMaterial, &EkinLab,
+                                                 &labMomentum, &iflxyz_) *
+                                1_mb;
+    return xs;
   }
-      
-      template <typename TSecondaryView>
-      inline void InteractionModel::doInteraction(
-      TSecondaryView& view, Code const projectileId, Code const targetId,
-      FourMomentum const& projectileP4, FourMomentum const& targetP4) {
-          
-      }
-  
-  
+
+  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());
+  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([&]() {
@@ -97,26 +126,26 @@ namespace corsika::fluka {
       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.
+
+    /*
+     * 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;    
+    auto const mxelfl = nElements;
     double const pptmax = 1e11; // GeV
-    double const ef2dp3 = 0; // GeV, 0 means default is used
-    double const df2dp3 = -1; // default
+    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
-    
-    
+    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
@@ -126,31 +155,34 @@ namespace corsika::fluka {
      *    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);
-    
+                   std::back_inserter(izelfl), get_nucleus_Z);
+
+    // check if FLUPRO is set
+    if (std::getenv("FLUPRO") == nullptr) {
+      throw std::runtime_error{"FLUPRO not set; required to initialize 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);
+
     // 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]));
+
+    auto it = boost::make_zip_iterator(
+        boost::make_tuple(allElementsInUniverse.begin(), &mtflka[0]));
+    auto const 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>());
+      boost::tuple<Code const&, int&> tup = *it;
+      mapping.emplace_back(tup.get<0>(), tup.get<1>());
     }
-    
-    return mapping;
 
-    }
-}
+    return mapping;
+  }
+} // namespace corsika::fluka

corsika/modules/FLUKA.hpp

+/*
+ * (c) Copyright 2020 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/modules/fluka/ParticleConversion.hpp>
+#include <corsika/modules/fluka/InteractionModel.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+/**
+ * @file Epos.hpp
+ *
+ * Includes all the parts of the EPOS model. Defines the InteractionProcess<TModel>
+ * classes needed for the ProcessSequence.
+ */
+
+namespace corsika::fluka {
+  /**
+   * fluka::Interaction is the process for ProcessSequence.
+   *
+   * The fluka::InteractionModel is wrapped as an InteractionProcess here in order
+   * to provide all the functions for ProcessSequence.
+   */
+  class Interaction : public fluka::InteractionModel,
+                      public corsika::InteractionProcess<Interaction> {};
+} // namespace corsika::fluka

corsika/modules/fluka/InteractionModel.hpp

@@ -8,6 +8,9 @@
 
 #pragma once
 
+#include <vector>
+#include <utility>
+
 #include <corsika/media/Environment.hpp>
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/framework/geometry/FourVector.hpp>
@@ -16,31 +19,33 @@
 
 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&);
+  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;
+    bool isValid(Code projectileID, int material, 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);
+
+  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
   };
-}
+
+  inline static int const iflxyz_ = 1;
+} // namespace corsika::fluka
 
 #include <corsika/detail/modules/fluka/InteractionModel.inl>

modules/fluka/FLUKA.hpp

@@ -9,68 +9,123 @@
 #pragma once
 
 namespace fluka {
-extern "C" {
-/*----------------------------------------------------------------------*
-*                                                                      *
-*     Copyright (C) 2022-2022      by    Alfredo Ferrari & Paola Sala  *
-*     All Rights Reserved.                                             *
-*                                                                      *
-*                                                                      *
-*     SeTuP for X/Y/Z interaction types:                               *
-*                                                                      *
-*     Authors:                           Alfredo Ferrari & Paola Sala  *
-*                                                                      *
-*                                                                      *
-*     Created on  24 October 2022  by    Alfredo Ferrari & Paola Sala  *
-*                                            Private        Private    *
-*                                                                      *
-*     Last change on  01-Nov-22    by             Alfredo Ferrari      *
-*                                                     Private          *
-*                                                                      *
-*     Input variables:                                                 *
-*                                                                      *
-*           Nmatfl    = Number of requested Fluka materials            *
-*           Nelmfl(m) = Number of elements in the m_th requested Fluka *
-*                       material, Nelmfl(m)=1 means simple element, no *
-*                       compound                                       *
-*           Izelfl(n) = Cumulative array of the Z of the elements con- *
-*                       stituting the requested Fluka materials        *
-*                      (the Z for the elements of the m_th materials   *
-*                       start at n = Sum_i=1^m-1 [Nelmfl(i)] + 1)      *
-*           Wfelml(n) = Cumulative array of the weight fractions of the*
-*                       elements constituting the requested Fluka      *
-*                       materials (the weight fractions for the elem-  *
-*                       ents of the m_th materials start at            *
-*                        n = Sum_i=1^m-1 [Nelmfl(i)] + 1)              *
-*           Mxelfl    = Dimension of the Iz/Wfelml arrays, it must be  *
-*                       Mxelfl >= Sum_i=1^nmatfl [Nelmfl(i)]           *
-*           Pptmax    = Maximum momentum (GeV/c) to be used in initia- *
-*                       lization (optional)                            *
-*           Ef2dp3    = Transition energy (GeV) from Fluka (Peanut) to *
-*                       Dpmjet3 for hA interactions (optional)         *
-*           Df2dp3    = Smearing (+/-Df2dp3) (GeV) of the transition   *
-*                       energy from Fluka (Peanut) to Dpmjet3 for hA   *
-*                       interactions (optional)                        *
-*           Lprint    = Material printout flag                         *
-*                                                                      *
-*     Output variables:                                                *
-*                                                                      *
-*           Mtflka(m) = Fluka material number corresponding to the m_th*
-*                       requested material                             *
-*                                                                      *
-*----------------------------------------------------------------------*
-*/
-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]);
-}
-}
+  extern "C" {
+  /*  Iflxyz must be consistent in all calls, to Stpxyz, Sgmxyz, Evtxyz
+   *    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
+   */
+
+  /*----------------------------------------------------------------------*
+   *                                                                      *
+   *     Copyright (C) 2022-2022      by    Alfredo Ferrari & Paola Sala  *
+   *     All Rights Reserved.                                             *
+   *                                                                      *
+   *                                                                      *
+   *     SeTuP for X/Y/Z interaction types:                               *
+   *                                                                      *
+   *     Authors:                           Alfredo Ferrari & Paola Sala  *
+   *                                                                      *
+   *                                                                      *
+   *     Created on  24 October 2022  by    Alfredo Ferrari & Paola Sala  *
+   *                                            Private        Private    *
+   *                                                                      *
+   *     Last change on  01-Nov-22    by             Alfredo Ferrari      *
+   *                                                     Private          *
+   *                                                                      *
+   *     Input variables:                                                 *
+   *                                                                      *
+   *           Nmatfl    = Number of requested Fluka materials            *
+   *           Nelmfl(m) = Number of elements in the m_th requested Fluka *
+   *                       material, Nelmfl(m)=1 means simple element, no *
+   *                       compound                                       *
+   *           Izelfl(n) = Cumulative array of the Z of the elements con- *
+   *                       stituting the requested Fluka materials        *
+   *                      (the Z for the elements of the m_th materials   *
+   *                       start at n = Sum_i=1^m-1 [Nelmfl(i)] + 1)      *
+   *           Wfelml(n) = Cumulative array of the weight fractions of the*
+   *                       elements constituting the requested Fluka      *
+   *                       materials (the weight fractions for the elem-  *
+   *                       ents of the m_th materials start at            *
+   *                        n = Sum_i=1^m-1 [Nelmfl(i)] + 1)              *
+   *           Mxelfl    = Dimension of the Iz/Wfelml arrays, it must be  *
+   *                       Mxelfl >= Sum_i=1^nmatfl [Nelmfl(i)]           *
+   *           Pptmax    = Maximum momentum (GeV/c) to be used in initia- *
+   *                       lization (optional)                            *
+   *           Ef2dp3    = Transition energy (GeV) from Fluka (Peanut) to *
+   *                       Dpmjet3 for hA interactions (optional)         *
+   *           Df2dp3    = Smearing (+/-Df2dp3) (GeV) of the transition   *
+   *                       energy from Fluka (Peanut) to Dpmjet3 for hA   *
+   *                       interactions (optional)                        *
+   *           Lprint    = Material printout flag                         *
+   *                                                                      *
+   *     Output variables:                                                *
+   *                                                                      *
+   *           Mtflka(m) = Fluka material number corresponding to the m_th*
+   *                       requested material                             *
+   *                                                                      *
+   *----------------------------------------------------------------------*/
+  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]);
+
+  /*----------------------------------------------------------------------*
+   *                                                                      *
+   *     Copyright (C) 2022-2022      by    Alfredo Ferrari & Paola Sala  *
+   *     All Rights Reserved.                                             *
+   *                                                                      *
+   *                                                                      *
+   *     EVenT for X/Y/Z interaction types:                               *
+   *                                                                      *
+   *     Authors:                           Alfredo Ferrari & Paola Sala  *
+   *                                                                      *
+   *                                                                      *
+   *     Created on  24 October 2022  by    Alfredo Ferrari & Paola Sala  *
+   *                                            Private        Private    *
+   *                                                                      *
+   *     Last change on  28-Oct-22    by             Alfredo Ferrari      *
+   *                                                     Private          *
+   *                                                                      *
+   *     Output variables:                                                *
+   *                                                                      *
+   *           Cumsgi(i) = cumulative macroscopic inelastic cross section *
+   *                       (cm^2/g x rho) for the i_th element of the mmat*
+   *                       material                                       *
+   *           Cumsge(i) = cumulative macroscopic elastic   cross section *
+   *                       (cm^2/g x rho) for the i_th element of the mmat*
+   *                       material                                       *
+   *           Cumsgm(i) = cumulative macroscopic emd       cross section *
+   *                       (cm^2/g x rho) for the i_th element of the mmat*
+   *                       material                                       *
+   *                                                                      *
+   *                                                                      *
+   *----------------------------------------------------------------------*/
+  //~ void evtxyz_ (int& , MMAT  , EKIN  , PPROJ , TXX, TYY, TZZ,
+  //~ &                    IFLXYZ, CUMSGI, CUMSGE, CUMSGM )
+
+  /*----------------------------------------------------------------------*
+   *                                                                      *
+   *     Copyright (C) 2022-2022      by    Alfredo Ferrari & Paola Sala  *
+   *     All Rights Reserved.                                             *
+   *                                                                      *
+   *     SiGMa (mb) for X/Y/Z interaction types:                          *
+   *                                                                      *
+   *     Authors:                           Alfredo Ferrari & Paola Sala  *
+   *                                                                      *
+   *                                                                      *
+   *     Created on  24 October 2022  by    Alfredo Ferrari & Paola Sala  *
+   *                                            Private        Private    *
+   *                                                                      *
+   *     Last change on  27-Oct-22    by             Alfredo Ferrari      *
+   *                                                     Private          *
+   *                                                                      *
+   *----------------------------------------------------------------------*/
+  double sgmxyz_(int const* KPROJ, int const* MMAT, double const* EKIN,
+                 double const* PPROJ, int const* IFLXYZ);
+  }
+} // namespace fluka

src/modules/fluka/fluka_codes.dat

+# FLUKA particle code conversion table
+# C8 identifier, FLUKA code, has cross-section/can interact
+Proton	1	yes
+AntiProton	2	yes
+Electron	3
+Positron	4
+NuE	5
+NuEBar	6
 Photon	7
-Electron	4
-Positron	3
-
+Neutron	8	yes
+AntiNeutron	9	yes
 MuPlus	10
 MuMinus	11
-Pi0	23
+K0Long	12	yes
 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
-# 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
+K0Short	19	yes
+SigmaMinus	20 yes
+SigmaPlus	21  yes
+Sigma0	22  yes
+Pi0	23  yes
+K0 24
+K0Bar 25
 NuMu	27
 NuMuBar	28
-
-D0	16
-DPlus	25
-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
+SigmaMinusBar	31  yes
+Sigma0Bar	32  yes
+SigmaPlusBar	33  yes
+Xi0	34  yes
+Xi0Bar	35  yes
+XiMinus	36  yes
+XiPlusBar	37  yes
+OmegaMinus	38  yes
+OmegaPlusBar	39  yes
 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
+DPlus	45
+DMinus	46
+D0	47
+D0Bar	48
+DsPlus	49
+DsMinus	50
+LambdaCPlus	51
+XiCPlus	52
+XiC0	53
+XiPrimeCPlus	54
+XiPrimeC0	55
+OmegaC0	56
+LambdaCMinusBar	57
+XiCMinusBar	58
+XiC0Bar	59
+XiPrimeCMinusBar	60
+XiPrimeC0Bar	61
+OmegaC0Bar	62

tests/modules/testFluka.cpp

@@ -6,10 +6,13 @@
  * the license.
  */
 
-#include <corsika/modules/fluka/InteractionModel.hpp>
+#include <corsika/modules/FLUKA.hpp>
 //~ #include <SetupTestEnvironment.hpp>
 
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
 #include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
 #include <corsika/framework/geometry/CoordinateSystem.hpp>
 
 #include <corsika/media/Environment.hpp>
@@ -18,21 +21,73 @@
 
 #include <catch2/catch.hpp>
 
+#include <fstream>
+#include <iomanip>
+
 using namespace corsika;
 
+TEST_CASE("FLUKACodeConversion") {
+  REQUIRE(corsika::fluka::convertToFluka(Code::PiPlus) ==
+          corsika::fluka::FLUKACode::PiPlus);
+  REQUIRE(corsika::fluka::convertToFlukaRaw(Code::PiPlus) == 13);
+  REQUIRE(corsika::fluka::convertToFlukaRaw(Code::Proton) == 1);
+  REQUIRE(corsika::fluka::convertToFlukaRaw(Code::Lambda0) == 17);
+}
+
 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};
-    
+  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>{.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},
+    };
+
+    HEPEnergyType const p = 100_GeV;
+    for (auto const& [projectileCode, targetCode] : combinations) {
+      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);
+    }
+  }
 }