event generation, test, linking

corsika/detail/modules/fluka/InteractionModel.inl

@@ -21,6 +21,8 @@
 #include <corsika/media/Environment.hpp>
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
 
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
@@ -90,8 +92,6 @@ namespace corsika::fluka {
 
     auto const plab = projectileLab4mom.getSpaceLikeComponents();
 
-    std::cout << targetRestBoost.boost_ << '\n';
-    std::cout << targetRestBoost.rotatedCS_->getTransform().matrix() << '\n';
     CORSIKA_LOGGER_DEBUG(logger_, fmt::format("Elab = {} GeV", Elab * invGeV));
     CORSIKA_LOGGER_DEBUG(logger_, fmt::format("EkinLab = {} GeV", EkinLab * invGeV));
 
@@ -139,24 +139,32 @@ namespace corsika::fluka {
                      cumsgx_.get() + materials_.size(),
                      cumsgx_.get() + materials_.size() * 2);
 
-    //~ extern struct {
-    //~ int nevhep;                  // event number
-    //~ int nhep;                    // number of entries
-    //~ hepmc_array<int> isthep;     // status code
-    //~ hepmc_array<int> idhep;      // PDG particle id
-    //~ hepmc_array<int[2]> jmohep;  // position of first, second mother
-    //~ hepmc_array<int[2]> jdahep;  // position of first, last daughter
-    //~ hepmc_array<double[5]> phep; // 4-momemtum, mass (GeV)
-    //~ hepmc_array<double[4]> vhep; // vertex, production time in mm
-    //~ } hepevt_;
-
     for (int i = 0; i < ::fluka::hepevt_.nhep; ++i) {
-      int const pdg = ::fluka::hepevt_.idhep[i];
       int const status = ::fluka::hepevt_.isthep[i];
-      auto const mom = QuantityVector<hepenergy_d>{
-          Eigen::Map<Eigen::Vector3d>(&::fluka::hepevt_.phep[i][0]) * invGeV.magnitude()};
+      // TODO: enable this when FLUKA writes status code
+      //~ if (status != 1) // skip non-final-state particles
+      //~ continue;
 
-      std::cout << pdg << '\t' << status << '\t' << mom << std::endl;
+      auto const pdg = static_cast<corsika::PDGCode>(::fluka::hepevt_.idhep[i]);
+      auto const c8code = corsika::convert_from_PDG(pdg);
+      auto const mom = QuantityVector<hepenergy_d>{
+          Eigen::Map<Eigen::Vector3d>(&::fluka::hepevt_.phep[i][0]) *
+          (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}};
+      auto const fourMomOrigFrame = targetRestBoost.fromCoM(fourMomCollisionFrame);
+      auto const momOrigFrame = fourMomOrigFrame.getSpaceLikeComponents();
+      auto const p = momOrigFrame.getNorm();
+
+      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;
     }
   }
 
@@ -195,8 +203,9 @@ namespace corsika::fluka {
     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
+    // magic number that FLUKA uses to see if it's the right version
+    char crvrck[] = "76466879";
+    int const size = 8;
 
     std::fill(&nelmfl[0], &nelmfl[nElements], 1);
     std::fill(&wfelml[0], &wfelml[nElements], 1.);
@@ -210,7 +219,8 @@ namespace corsika::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,
+                     &size);
 
     // now create & fill vector of (C8 Code, FLUKA mat. no.) pairs
     std::vector<std::pair<Code, int>> mapping;

corsika/modules/FLUKA.hpp

@@ -13,9 +13,9 @@
 #include <corsika/framework/process/InteractionProcess.hpp>
 
 /**
- * @file Epos.hpp
+ * @file FLUKA.hpp
  *
- * Includes all the parts of the EPOS model. Defines the InteractionProcess<TModel>
+ * Includes all the parts of the FLUKA model. Defines the InteractionProcess<TModel>
  * classes needed for the ProcessSequence.
  */
 

corsika/modules/Random.hpp

@@ -22,3 +22,4 @@
 #include <corsika/modules/urqmd/Random.hpp>
 #include <corsika/modules/qgsjetII/Random.hpp>
 #include <corsika/modules/conex/Random.hpp>
+#include <corsika/modules/fluka/Random.hpp>

corsika/modules/fluka/InteractionModel.hpp

@@ -18,6 +18,7 @@
 #include <corsika/framework/utility/COMBoost.hpp>
 #include <corsika/framework/core/Logging.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
 
 namespace corsika::fluka {
   class InteractionModel {
@@ -39,14 +40,14 @@ namespace corsika::fluka {
                        FourMomentum const& projectileP4, FourMomentum const& targetP4);
 
   private:
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("fluka");
     std::vector<std::pair<Code, int>> const
         materials_; //!< map target Code to FLUKA material no.
     std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_FLUKA_Interaction");
-    std::unique_ptr<double[]> cumsgx_; //!< dump for evtxyz cumsg*
+    std::unique_ptr<double[]> cumsgx_; //!< dump for evtxyz cumsg*, never read again
 
     template <typename TEnvironment>
     static std::vector<std::pair<Code, int>> genFlukaMaterials(TEnvironment const&);
-    // TODO: random number stream
   };
 
   inline static int const iflxyz_ = 1;

corsika/modules/fluka/ParticleConversion.hpp

@@ -27,13 +27,14 @@ namespace corsika::fluka {
 
   FLUKACode constexpr convertToFluka(Code const c8id) {
     if (is_nucleus(c8id)) {
-		throw std::runtime_error{"nucleus conversion to FLUKA not implemented"};
-	}
+      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()};
-	}
+      throw std::runtime_error{
+          fmt::format("no correspondig FLUKA id for {}", get_name(c8id)).c_str()};
+    }
 
     return flukaID;
   }

corsika/modules/fluka/Random.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 <random>
+
+#include <corsika/framework/random/RNGManager.hpp>
+
+namespace fluka {
+  /**
+   * the random number generator function for FLUKA
+   */
+  double rndm_interface() {
+    static corsika::default_prng_type& rng =
+        corsika::RNGManager<>::getInstance().getRandomStream("fluka");
+    static std::uniform_real_distribution<double> dist;
+    return dist(rng);
+  }
+
+} // namespace fluka

modules/fluka/CMakeLists.txt

 set(C8_FLUKALIB CACHE STRING "path to libflukahp.a")
 
 if (DEFINED C8_FLUKALIB)
-    add_library(flukahp STATIC IMPORTED)
-    set_property(TARGET flukahp PROPERTY IMPORTED_LOCATION ${C8_FLUKALIB})
+    set (input_dir ${PROJECT_SOURCE_DIR}/src/modules/fluka)
+    
+    # we remove flrndm_.o from the original libflukahp.a and save the result as libflukahp-norndm.a
+    add_custom_command(OUTPUT ${CMAKE_BINARY_DIR}/libflukahp-norndm.a COMMAND ${input_dir}/strip_flukahp.sh ${C8_FLUKALIB} ${CMAKE_BINARY_DIR})
+    add_custom_target(generate_libfluka-norndm DEPENDS ${CMAKE_BINARY_DIR}/libflukahp-norndm.a)
+
+    add_library(flukahp-norndm STATIC IMPORTED)
+    add_dependencies(flukahp-norndm generate_libfluka-norndm)
+    set_property(TARGET flukahp-norndm PROPERTY IMPORTED_LOCATION ${CMAKE_BINARY_DIR}/libflukahp-norndm.a)
     find_library(MATH_LIBRARY m)
-    target_link_libraries(flukahp INTERFACE gfortran ${MATH_LIBRARY})
-    target_include_directories (flukahp INTERFACE
+    target_link_libraries(flukahp-norndm INTERFACE gfortran ${MATH_LIBRARY})
+    target_include_directories (flukahp-norndm INTERFACE
         $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
         $<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)
+
+    target_include_directories (flukaRndm PUBLIC
+      $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+      $<INSTALL_INTERFACE:include/corsika_modules/fluka>
+      )
+    
+    set_target_properties (
+      flukaRndm
+      PROPERTIES
+      POSITION_INDEPENDENT_CODE 1
+      )
+    
+    install (
+      FILES
+      fluka.hpp
+      DESTINATION include/corsika_modules/urqmd
+      )
+    
+    install (
+      TARGETS flukaRndm
+      EXPORT CORSIKA8PublicTargets
+      ARCHIVE DESTINATION lib/corsika
+      INCLUDES DESTINATION include/corsika_modules/fluka
+      )
+    
     # add fluka to corsika8 build
-    add_dependencies (CORSIKA8 flukahp)
-    target_link_libraries (CORSIKA8 INTERFACE flukahp)
+    add_dependencies (CORSIKA8 flukaRndm)
+    target_link_libraries (CORSIKA8 INTERFACE flukaRndm)
 endif()

modules/fluka/FLUKA.hpp

@@ -8,7 +8,19 @@
 
 #pragma once
 
+#include <cstddef>
+#include <array>
+
 namespace fluka {
+  /**
+   * \fluka fluka::rndm_interface
+   *
+   * this is the random number hook to external packages.
+   *
+   * CORSIKA8, for example, has to provide an implementation of this.
+   **/
+  extern double rndm_interface();
+
   size_t constexpr nmxhep = 10000;
   template <typename T>
   using hepmc_array = std::array<T, nmxhep>;
@@ -116,7 +128,7 @@ namespace fluka {
   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]);
+               bool const* LPRINT, int* MTFLKA, char const* CRVRCK, int const*);
 
   /*----------------------------------------------------------------------*
    *                                                                      *
@@ -191,5 +203,7 @@ namespace fluka {
    *                                                                      *
    *----------------------------------------------------------------------*/
   void fllhep_();
+  
+  double flrndm();
   }
 } // namespace fluka

modules/fluka/flrndm.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 <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();
+  }
+}
+
+}

src/modules/fluka/CMakeLists.txt

@@ -23,4 +23,4 @@ add_dependencies (CORSIKA8 SourceDirLinkFLUKA)
 install (
   FILES ${output_dir}/Generated.inc
   DESTINATION include/corsika/modules/fluka
-  )
+)

src/modules/fluka/strip_flukahp.sh

+#!/bin/sh
+
+# (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.
+
+# This script strips off flrndm() from the libflukahp.a so that we can provide our own
+# implementation.
+
+flukalibOrig="$1"
+target="$2"
+
+tmpdir=`mktemp -d fluka_objectsXXXXXX`
+
+echo "extracting objects from $1 into `realpath $tmpdir`..."
+ar --output "$tmpdir" x "$flukalibOrig"
+rm "$tmpdir/flrndm.o"
+
+[ -f "libflukahp-norndm.a" ] && rm "libflukahp-norndm.a"
+
+echo "creating libflukahp-norndm.a..."
+ar -rcs "$target/libflukahp-norndm.a" "$tmpdir"/*.o
+
+rm -r "$tmpdir"

tests/modules/testFluka.cpp

@@ -7,6 +7,7 @@
  */
 
 #include <corsika/modules/FLUKA.hpp>
+#include <corsika/modules/fluka/Random.hpp>
 //~ #include <SetupTestEnvironment.hpp>
 
 #include <corsika/framework/core/EnergyMomentumOperations.hpp>
@@ -29,6 +30,13 @@
 
 using namespace corsika;
 
+template <typename TStackView>
+auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
+  MomentumVector sum{vCS, 0_eV, 0_eV, 0_eV};
+  for (auto const& p : view) { sum += p.getMomentum(); }
+  return sum;
+}
+
 TEST_CASE("FLUKACodeConversion") {
   REQUIRE(corsika::fluka::convertToFluka(Code::PiPlus) ==
           corsika::fluka::FLUKACode::PiPlus);
@@ -44,6 +52,8 @@ TEST_CASE("FLUKA") {
   using MyHomogeneousModel = MediumPropertyModel<
       UniformMagneticField<HomogeneousMedium<DummyEnvironmentInterface>>>;
 
+  RNGManager<>::getInstance().registerRandomStream("fluka");
+
   DummyEnvironment env;
   auto& universe = *env.getUniverse();
   CoordinateSystemPtr const& cs = env.getCoordinateSystem();
@@ -99,20 +109,35 @@ TEST_CASE("FLUKA") {
   {
     auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
     auto const& cs = *csPtr;
-    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 projectileCode = Code::PiPlus;
-    auto const targetCode = Code::Nitrogen;
-    auto const p = 20_GeV;
-    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 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));
+      }
+    }
   }
 }