Processes/LongitudinalProfile/CMakeLists.txt

-set (
-  MODEL_SOURCES
-  LongitudinalProfile.cc
-  )
-
-set (
-  MODEL_HEADERS
-  LongitudinalProfile.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/longitudinal_profile
-  )
-
-add_library (ProcessLongitudinalProfile STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessLongitudinalProfile ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessLongitudinalProfile
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessLongitudinalProfile
-  CORSIKAenvironment
-  CORSIKAunits
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAsetup
-  )
-
-target_include_directories (
-  ProcessLongitudinalProfile 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessLongitudinalProfile
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-# CORSIKA_ADD_TEST(testNullModel)
-#target_link_libraries (
-#  testNullModel  ProcessNullModel
-#  CORSIKAsetup
-#  CORSIKAgeometry
-#  CORSIKAunits
-#  CORSIKAthirdparty # for catch2
-#  )
-

Processes/LongitudinalProfile/LongitudinalProfile.cc

-/*
- * (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.
- */
-
-#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/logging/Logging.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <cmath>
-#include <iomanip>
-#include <limits>
-
-using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-using Track = Trajectory;
-
-using namespace corsika::process::longitudinal_profile;
-using namespace corsika::units::si;
-
-LongitudinalProfile::LongitudinalProfile(environment::ShowerAxis const& shower_axis,
-                                         units::si::GrammageType dX)
-    : dX_(dX)
-    , shower_axis_{shower_axis}
-    , profiles_{static_cast<unsigned int>(shower_axis.maximumX() / dX_) + 1} {}
-
-template <>
-corsika::process::EProcessReturn LongitudinalProfile::DoContinuous(Particle const& vP,
-                                                                   Track const& vTrack) {
-  auto const pid = vP.GetPID();
-
-  GrammageType const grammageStart = shower_axis_.projectedX(vTrack.GetPosition(0));
-  GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1));
-
-  C8LOG_INFO(
-      "pos1={} m, pos2={}, X={} g/cm2", vTrack.GetPosition(0).GetCoordinates() / 1_m,
-      vTrack.GetPosition(1).GetCoordinates() / 1_m, grammageStart / 1_g * square(1_cm));
-
-  const int binStart = std::ceil(grammageStart / dX_);
-  const int binEnd = std::floor(grammageEnd / dX_);
-
-  for (int b = binStart; b <= binEnd; ++b) {
-    if (pid == particles::Code::Gamma) {
-      profiles_.at(b)[ProfileIndex::Gamma]++;
-    } else if (pid == particles::Code::Positron) {
-      profiles_.at(b)[ProfileIndex::Positron]++;
-    } else if (pid == particles::Code::Electron) {
-      profiles_.at(b)[ProfileIndex::Electron]++;
-    } else if (pid == particles::Code::MuPlus) {
-      profiles_.at(b)[ProfileIndex::MuPlus]++;
-    } else if (pid == particles::Code::MuMinus) {
-      profiles_.at(b)[ProfileIndex::MuMinus]++;
-    } else if (particles::IsHadron(pid)) {
-      profiles_.at(b)[ProfileIndex::Hadron]++;
-    }
-  }
-
-  return corsika::process::EProcessReturn::eOk;
-}
-
-void LongitudinalProfile::save(std::string const& filename, const int width,
-                               const int precision) {
-  std::ofstream f{filename};
-  f << "# X / g·cm¯², gamma, e+, e-, mu+, mu-, all hadrons" << std::endl;
-  for (size_t b = 0; b < profiles_.size(); ++b) {
-    f << std::setprecision(5) << std::setw(11) << b * (dX_ / (1_g / 1_cm / 1_cm));
-    for (auto const& N : profiles_.at(b)) {
-      f << std::setw(width) << std::setprecision(precision) << std::scientific << N;
-    }
-    f << std::endl;
-  }
-}

Processes/LongitudinalProfile/LongitudinalProfile.h

-/*
- * (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/environment/ShowerAxis.h>
-#include <corsika/process/ContinuousProcess.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <array>
-#include <fstream>
-#include <limits>
-#include <string>
-
-namespace corsika::process::longitudinal_profile {
-
-  /**
-   * \class LongitudinalProfile
-   *
-   * is a ContinuousProcess, which is constructed from an environment::ShowerAxis
-   * object, and a dX in units of g/cm2
-   * (corsika::units::si::GrammageType).
-   *
-   * LongitudinalProfile does then convert each single Track of the
-   * simulation into a projected grammage range and counts for
-   * different particle species when they cross dX (default: 10g/cm2)
-   * boundaries.
-   *
-   **/
-
-  class LongitudinalProfile
-      : public corsika::process::ContinuousProcess<LongitudinalProfile> {
-
-  public:
-    LongitudinalProfile(environment::ShowerAxis const&,
-                        units::si::GrammageType dX = std::invoke([]() {
-                          using namespace units::si;
-                          return 10_g / square(1_cm);
-                        })); // profile binning);
-
-    template <typename TParticle, typename TTrack>
-    corsika::process::EProcessReturn DoContinuous(TParticle const&, TTrack const&);
-
-    template <typename TParticle, typename TTrack>
-    corsika::units::si::LengthType MaxStepLength(TParticle const&, TTrack const&) {
-      return units::si::meter * std::numeric_limits<double>::infinity();
-    }
-
-    void save(std::string const&, int const width = 14, int const precision = 6);
-
-  private:
-    units::si::GrammageType const dX_;
-
-    environment::ShowerAxis const& shower_axis_;
-    using ProfileEntry = std::array<uint32_t, 6>;
-    enum ProfileIndex {
-      Gamma = 0,
-      Positron = 1,
-      Electron = 2,
-      MuPlus = 3,
-      MuMinus = 4,
-      Hadron = 5
-    };
-    std::vector<ProfileEntry> profiles_; // longitudinal profile
-  };
-
-} // namespace corsika::process::longitudinal_profile

Processes/ObservationPlane/ObservationPlane_interpolation.cc

-/*
- * (c) Copyright 2019 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.
- */
-
-#include <corsika/process/observation_plane/ObservationPlane.h>
-
-#include <fstream>
-
-using namespace corsika::process::observation_plane;
-using namespace corsika::units::si;
-
-ObservationPlane::ObservationPlane(
-    geometry::Plane const& obsPlane,
-    geometry::Vector<units::si::dimensionless_d> const& x_axis,
-    std::string const& filename, bool deleteOnHit)
-    : plane_(obsPlane)
-    , outputStream_(filename)
-    , deleteOnHit_(deleteOnHit)
-    , xAxis_(x_axis.normalized())
-    , yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
-  outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;
-}
-
-corsika::process::EProcessReturn ObservationPlane::DoContinuous(
-    setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
-  TimeType const timeOfIntersection =
-      (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
-      trajectory.GetV0().dot(plane_.GetNormal());
-
-  if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
-
-  if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
-    return process::EProcessReturn::eOk;
-  }
-
-  auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
-
-  outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
-                << particle.GetEnergy() * (1 / 1_eV) << ' '
-                << displacement.dot(xAxis_) / 1_m << ' ' << displacement.dot(yAxis_) / 1_m
-                << ' ' << std::endl;
-
-  if (deleteOnHit_) {
-    return process::EProcessReturn::eParticleAbsorbed;
-  } else {
-    return process::EProcessReturn::eOk;
-  }
-}
-
-LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
-                                           setup::Trajectory const& trajectory) {
-  TimeType const timeOfIntersection =
-      (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
-      trajectory.GetV0().dot(plane_.GetNormal());
-
-  if (timeOfIntersection < TimeType::zero()) {
-    return std::numeric_limits<double>::infinity() * 1_m;
-  }
-
-  auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
-  return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
-}

Processes/OnShellCheck/CMakeLists.txt

-set (
-  MODEL_SOURCES
-  OnShellCheck.cc
-)
-
-set (
-  MODEL_HEADERS
-  OnShellCheck.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/on_shell_check
-  )
-
-add_library (ProcessOnShellCheck STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessOnShellCheck ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessOnShellCheck
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessOnShellCheck
-  CORSIKAunits
-  CORSIKAparticles
-  CORSIKAprocesssequence
-  CORSIKAsetup
-  )
-
-target_include_directories (
-  ProcessOnShellCheck 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessOnShellCheck
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
-  )
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST(testOnShellCheck SOURCES
-  testOnShellCheck.cc
-  ${MODEL_HEADERS}
-)
-
-target_link_libraries (
-  testOnShellCheck
-  ProcessOnShellCheck
-  CORSIKAunits
-  CORSIKAstackinterface
-  CORSIKAprocesssequence
-  CORSIKAsetup
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAtesting
-  )

Processes/OnShellCheck/OnShellCheck.cc

-/*
- * (c) Copyright 2018 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/geometry/FourVector.h>
-#include <corsika/process/on_shell_check/OnShellCheck.h>
-
-using namespace std;
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units::si;
-using namespace corsika::particles;
-using namespace corsika::setup;
-
-namespace corsika::process {
-  namespace on_shell_check {
-
-    void OnShellCheck::Init() {
-      std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100
-                << "%" << endl
-                << "              energy tolerance is set to " << energy_tolerance_ * 100
-                << "%" << std::endl;
-    }
-
-    EProcessReturn OnShellCheck::DoSecondaries(corsika::setup::StackView& vS) {
-      for (auto& p : vS) {
-        auto const pid = p.GetPID();
-        if (!particles::IsHadron(pid) || particles::IsNucleus(pid)) continue;
-        auto const e_original = p.GetEnergy();
-        auto const p_original = p.GetMomentum();
-        auto const Plab = corsika::geometry::FourVector(e_original, p_original);
-        auto const m_kinetic = Plab.GetNorm();
-        auto const m_corsika = particles::GetMass(pid);
-        auto const m_err_abs = abs(m_kinetic - m_corsika);
-        if (m_err_abs >= mass_tolerance_ * m_corsika) {
-          const HEPEnergyType e_shifted =
-              sqrt(p_original.GetSquaredNorm() + m_corsika * m_corsika);
-          auto const e_shift_relative = (e_shifted / e_original - 1);
-          count_ = count_ + 1;
-          average_shift_ += abs(e_shift_relative);
-          if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
-          std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
-                    << std::setw(40) << std::setfill(' ')
-                    << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
-                    << std::setw(40) << std::setfill(' ')
-                    << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
-                    << std::setw(40) << std::setfill(' ')
-                    << "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
-                    << std::setw(40) << std::setfill(' ')
-                    << "mass tolerance (GeV): " << (m_corsika * mass_tolerance_) / 1_GeV
-                    << std::endl;
-          /*
-            For now we warn if the necessary shift is larger than 1%.
-            we could promote this to an error.
-          */
-          if (abs(e_shift_relative) > energy_tolerance_) {
-            std::cout << "OnShellCheck: warning! shifted particle energy by "
-                      << e_shift_relative * 100 << " %" << std::endl;
-            if (throw_error_)
-              throw std::runtime_error(
-                  "OnShellCheck: error! shifted energy by large amount!");
-          }
-
-          // reset energy
-          p.SetEnergy(e_shifted);
-        } else
-          std::cout << "OnShellCheck: particle mass for " << pid << " OK" << std::endl;
-      }
-      return EProcessReturn::eOk;
-    }
-  } // namespace on_shell_check
-} // namespace corsika::process

Processes/OnShellCheck/OnShellCheck.h

-/*
- * (c) Copyright 2018 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.
- */
-
-#ifndef _corsika_process_on_shell_check_OnShellCheck_h_
-#define _corsika_process_on_shell_check_OnShellCheck_h_
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/SecondariesProcess.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-  namespace on_shell_check {
-    class OnShellCheck : public process::SecondariesProcess<OnShellCheck> {
-      double average_shift_ = 0;
-      double max_shift_ = 0;
-      double count_ = 0;
-
-    public:
-      OnShellCheck(const double vMassTolerance, const double vEnergyTolerance,
-                   const bool vError)
-          : mass_tolerance_(vMassTolerance)
-          , energy_tolerance_(vEnergyTolerance)
-          , throw_error_(vError) {}
-
-      ~OnShellCheck() {
-        std::cout << "OnShellCheck: summary" << std::endl
-                  << " particles shifted: " << int(count_) << std::endl;
-        if (count_)
-          std::cout << " average energy shift (%): " << average_shift_ / count_ * 100.
-                    << std::endl
-                    << " max. energy shift (%): " << max_shift_ * 100. << std::endl;
-      };
-
-      EProcessReturn DoSecondaries(corsika::setup::StackView&);
-
-      void Init();
-
-    private:
-      double mass_tolerance_;
-      double energy_tolerance_;
-      bool throw_error_;
-    };
-  } // namespace on_shell_check
-} // namespace corsika::process
-
-#endif

Processes/OnShellCheck/testOnShellCheck.cc

-/*
- * (c) Copyright 2018 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/process/on_shell_check/OnShellCheck.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/geometry/FourVector.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <corsika/setup/SetupStack.h>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::process::on_shell_check;
-using namespace corsika::units;
-using namespace corsika::units::si;
-
-TEST_CASE("OnShellCheck", "[processes]") {
-  feenableexcept(FE_INVALID);
-  using EnvType = setup::Environment;
-  EnvType env;
-  const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  // setup empty particle stack
-  setup::Stack stack;
-  stack.Clear();
-  // two energies
-  const HEPEnergyType E = 10_GeV;
-  // list of arbitrary particles
-  std::array const particleList{particles::Code::PiPlus, particles::Code::PiMinus,
-                                particles::Code::Helium, particles::Code::Gamma};
-
-  std::array const mass_shifts{1.1, 1.001, 1.0, 1.0};
-
-  SECTION("check particle masses") {
-
-    OnShellCheck check(1.e-2, 0.01, false);
-
-    check.Init();
-
-    // add primary particle to stack
-    auto particle = stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            particles::Code::Proton, E,
-            corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
-            geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-    // view on secondary particles
-    setup::StackView view{particle};
-    // ref. to primary particle through the secondary view.
-    // only this way the secondary view is populated
-    auto projectile = view.GetProjectile();
-    // add secondaries, all with energies above the threshold
-    // only cut is by species
-    int count = -1;
-    for (auto proType : particleList) {
-      count++;
-      const auto pz = sqrt((E - particles::GetMass(proType) * mass_shifts[count]) *
-                           (E + particles::GetMass(proType) * mass_shifts[count]));
-      auto const momentum = corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
-      projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                         corsika::stack::MomentumVector, geometry::Point,
-                                         units::si::TimeType>{
-          proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-    }
-    check.DoSecondaries(view);
-    int i = -1;
-    for (auto& p : view) {
-      i++;
-      auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
-      auto const m_kinetic = Plab.GetNorm();
-      if (i == 0)
-        CHECK(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
-      else if (i == 1)
-        CHECK_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
-    }
-  }
-}

Processes/Proposal/CMakeLists.txt

-file(READ CMakeLists_PROPOSAL.txt overwrite_CMakeLists_PROPOSAL_txt)
-file(WRITE PROPOSAL/CMakeLists.txt ${overwrite_CMakeLists_PROPOSAL_txt})
-add_subdirectory (PROPOSAL)
-
-FILE (GLOB MODEL_SOURCES *.cc)
-FILE (GLOB MODEL_HEADERS *.h)
-SET (MODEL_NAMESPACE corsika/process/proposal)
-
-ADD_LIBRARY (ProcessProposal STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessProposal ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-SET_TARGET_PROPERTIES (
-  ProcessProposal
-  PROPERTIES VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  )
-
-TARGET_LINK_LIBRARIES (
-  ProcessProposal
-  CORSIKAparticles
-  CORSIKAutilities
-  CORSIKAunits
-  CORSIKAthirdparty
-  CORSIKAgeometry
-  CORSIKAlogging
-  CORSIKAenvironment
-  ProcessParticleCut
-  PROPOSAL::PROPOSAL
-  )
-
-TARGET_INCLUDE_DIRECTORIES (
-  ProcessProposal
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-INSTALL (
-  TARGETS ProcessProposal
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  )

Processes/Proposal/CMakeLists_PROPOSAL.txt

-
-cmake_minimum_required(VERSION 3.8)
-
-if(CMAKE_PROJECT_NAME STREQUAL corsika)
-  message(STATUS "Including PROPOSAL as part of CORSIKA8")
-  set (IN_CORSIKA8 ON)
-else(CMAKE_PROJECT_NAME STREQUAL corsika)
-  set (IN_CORSIKA8 OFF)
-endif(CMAKE_PROJECT_NAME STREQUAL corsika)
-
-project(PROPOSAL VERSION 6.1.2 LANGUAGES CXX)
-
-if (NOT IN_CORSIKA8)
-IF(APPLE)
-    # In newer version of cmake this will be the default
-    SET(CMAKE_MACOSX_RPATH 1)
-ENDIF(APPLE)
-
-# sets standard installtion paths
-include(GNUInstallDirs)
-
-### full RPATH
-### copied from https://cmake.org/Wiki/CMake_RPATH_handling
-### set the RPATH so that for using PROPOSAL in python
-### the DYLD_LIBRARY_PATH must not be set in the bashrc
-### But for using PROPOSAL as c-Library, this path still
-### has to be set
-
-# use, i.e. don't skip the full RPATH for the build tree
-SET(CMAKE_SKIP_BUILD_RPATH  FALSE)
-
-# when building, don't use the install RPATH already
-# (but later on when installing)
-OPTION(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
-list(APPEND CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib;${CMAKE_INSTALL_PREFIX}/lib64")
-message(STATUS "CMAKE_INSTALL_RPATH=${CMAKE_INSTALL_RPATH}")
-
-# add the automatically determined parts of the RPATH
-# which point to directories outside the build tree to the install RPATH
-OPTION(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
-### end full RPATH
-
-# Some additional options
-OPTION(ADD_PYTHON "Choose to compile the python wrapper library" ON)
-endif (NOT IN_CORSIKA8)
-
-#################################################################
-####################       PROPOSAL      ########################
-#################################################################
-file(GLOB_RECURSE SRC_FILES ${PROJECT_SOURCE_DIR}/src/PROPOSAL/*)
-if (IN_CORSIKA8)
-  add_library(PROPOSAL STATIC ${SRC_FILES})
-else (IN_CORSIKA8)
-  add_library(PROPOSAL SHARED ${SRC_FILES})
-endif (IN_CORSIKA8)
-
-add_library(PROPOSAL::PROPOSAL ALIAS PROPOSAL)
-#target_compile_features(PROPOSAL PUBLIC cxx_std_11)
-#set_target_properties(PROPOSAL PROPERTIES CXX_EXTENSIONS OFF)
-target_include_directories(
-    PROPOSAL PUBLIC
-    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
-    $<INSTALL_INTERFACE:include>
-)
-target_compile_options(PROPOSAL PRIVATE -Wall -Wextra -Wnarrowing -Wpedantic -fdiagnostics-show-option -Wno-format-security)
-install(
-    TARGETS PROPOSAL
-    EXPORT PROPOSALTargets
-    LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
-    ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
-)
-
-# input version from the project call, so the library knows its own version
-configure_file(
-    "${PROJECT_SOURCE_DIR}/include/PROPOSAL/version.h.in"
-    "${PROJECT_BINARY_DIR}/include/PROPOSAL/version.h"
-)
-install(
-    FILES ${PROJECT_BINARY_DIR}/include/PROPOSAL/version.h
-    DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/include/PROPOSAL
-)
-target_include_directories(
-    PROPOSAL PUBLIC
-    $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-    $<INSTALL_INTERFACE:include>
-)
-
-# install header files
-file(GLOB_RECURSE INC_FILES ${PROJECT_SOURCE_DIR}/include/*)
-foreach(INC_FILE ${INC_FILES})
-    file(RELATIVE_PATH REL_FILE ${PROJECT_SOURCE_DIR}/include ${INC_FILE})
-    get_filename_component(DIR ${REL_FILE} DIRECTORY)
-    install(FILES include/${REL_FILE} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${DIR})
-endforeach()
-
-
-if(NOT IS_SYMLINK ${CMAKE_BINARY_DIR}/resources)
-    execute_process(COMMAND ln -sv ${CMAKE_CURRENT_SOURCE_DIR}/resources ${CMAKE_BINARY_DIR}/resources OUTPUT_VARIABLE link_msg OUTPUT_STRIP_TRAILING_WHITESPACE)
-    message(STATUS "Symlink to resources created:")
-    message(STATUS "  ${link_msg}")
-endif()
-
-
-#################################################################
-#################           spdlog        #######################
-#################################################################
-
-if (NOT IN_CORSIKA8)
-set(CMAKE_POSITION_INDEPENDENT_CODE ON)
-add_subdirectory("vendor/spdlog")
-endif (NOT IN_CORSIKA8)
-target_link_libraries(PROPOSAL PRIVATE spdlog::spdlog)
-
-
-#################################################################
-#################           Tests        ########################
-#################################################################
-
-if (NOT IN_CORSIKA8)
-if(CMAKE_PROJECT_NAME STREQUAL PROPOSAL)
-    include(CTest)
-endif()
-if(CMAKE_PROJECT_NAME STREQUAL PROPOSAL AND BUILD_TESTING)
-    add_subdirectory(tests)
-else()
-    MESSAGE(STATUS "No tests will be build.")
-endif()
-endif (NOT IN_CORSIKA8)
-
-#################################################################
-#################      Documentation     ########################
-#################################################################
-if (NOT IN_CORSIKA8)
-add_subdirectory(doc)
-endif (NOT IN_CORSIKA8)
-
-
-#################################################################
-#################           python      #########################
-#################################################################
-if (NOT IN_CORSIKA8)
-IF(ADD_PYTHON)
-    message(STATUS "Building the python wrapper library.")
-    find_package(PythonLibs REQUIRED)
-    add_subdirectory("vendor/pybind11")
-    file(GLOB_RECURSE PYTHON_SRC_FILES ${PROJECT_SOURCE_DIR}/interface/python/*)
-
-    pybind11_add_module(pyproposal SHARED ${PYTHON_SRC_FILES})
-    set_target_properties(pyproposal PROPERTIES OUTPUT_NAME proposal)
-    target_include_directories(pyproposal PRIVATE ${PYTHON_INCLUDE_DIRS})
-    target_include_directories(pyproposal PRIVATE ${PROJECT_SOURCE_DIR}/interface/python)
-    target_link_libraries(pyproposal PRIVATE PROPOSAL)
-    target_compile_options(pyproposal PRIVATE -fvisibility=hidden)
-    install(TARGETS pyproposal EXPORT PROPOSALTargets DESTINATION ${CMAKE_INSTALL_LIBDIR})
-ELSE(ADD_PYTHON)
-    MESSAGE(STATUS "No python wrapper library will be build.")
-ENDIF(ADD_PYTHON)
-endif (NOT IN_CORSIKA8)
-
-
-
-#################################################################
-#################         INSTALLATION        ###################
-#################################################################
-if (NOT IN_CORSIKA8)
-install(
-    EXPORT PROPOSALTargets
-    FILE PROPOSALConfig.cmake
-    NAMESPACE PROPOSAL::
-    DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/PROPOSAL
-)
-
-include(CMakePackageConfigHelpers)
-write_basic_package_version_file(
-    "PROPOSALConfigVersion.cmake"
-    VERSION ${PACKAGE_VERSION}
-    COMPATIBILITY SameMajorVersion
-)
-install(
-    FILES "${CMAKE_CURRENT_BINARY_DIR}/PROPOSALConfigVersion.cmake"
-    DESTINATION lib/cmake/PROPOSAL
-)
-endif (NOT IN_CORSIKA8)

Processes/Proposal/ContinuousProcess.cc

-/*
- * (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.
- */
-
-#include <PROPOSAL/PROPOSAL.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/process/proposal/ContinuousProcess.h>
-#include <corsika/process/proposal/Interaction.h>
-#include <corsika/setup/SetupEnvironment.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/COMBoost.h>
-#include <corsika/logging/Logging.h>
-
-#include <iostream>
-
-namespace corsika::process::proposal {
-
-  void ContinuousProcess::BuildCalculator(particles::Code code,
-                                          environment::NuclearComposition const& comp) {
-    // search crosssection builder for given particle
-    auto p_cross = cross.find(code);
-    if (p_cross == cross.end())
-      throw std::runtime_error("PROPOSAL could not find corresponding builder");
-
-    // interpolate the crosssection for given media and energy cut. These may
-    // take some minutes if you have to build the tables and cannot read the
-    // from disk
-    auto c = p_cross->second(media.at(comp.hash()), emCut_);
-
-    // Build displacement integral and scattering object and interpolate them too and
-    // saved in the calc map by a key build out of a hash of composed of the component and
-    // particle code.
-    auto disp = PROPOSAL::make_displacement(c, true);
-    auto scatter =
-        PROPOSAL::make_scattering("highland", particle[code], media.at(comp.hash()));
-    calc[std::make_pair(comp.hash(), code)] =
-        std::make_tuple(std::move(disp), std::move(scatter));
-  }
-
-  template <>
-  ContinuousProcess::ContinuousProcess(setup::Environment const& _env,
-                                       corsika::units::si::HEPEnergyType _emCut)
-      : ProposalProcessBase(_env, _emCut) {}
-
-  template <>
-  void ContinuousProcess::Scatter(setup::Stack::ParticleType& vP,
-                                  corsika::units::si::HEPEnergyType const& loss,
-                                  corsika::units::si::GrammageType const& grammage) {
-    using namespace corsika::units::si; // required for operator::_MeV
-
-    // Get or build corresponding calculators
-    auto c = GetCalculator(vP, calc);
-
-    // Cast corsika vector to proposal vector
-    auto vP_dir = vP.GetDirection();
-    auto d = vP_dir.GetComponents();
-    auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
-                                        d.GetZ().magnitude());
-
-    auto E_f = vP.GetEnergy() - loss;
-
-    // draw random numbers required for scattering process
-    std::uniform_real_distribution<double> distr(0., 1.);
-    auto rnd = array<double, 4>();
-    for (auto& it : rnd) it = distr(fRNG);
-
-    // calculate deflection based on particle energy, loss
-    auto [mean_dir, final_dir] = get<eSCATTERING>(c->second)->Scatter(
-        grammage / 1_g * square(1_cm), vP.GetEnergy() / 1_MeV, E_f / 1_MeV, direction,
-        rnd);
-
-    // TODO: neglect mean direction deflection because Trajectory is a const ref
-    (void)mean_dir;
-
-    // update particle direction after continuous loss caused by multiple
-    // scattering
-    auto vec = corsika::geometry::QuantityVector(
-        final_dir.GetX() * E_f, final_dir.GetY() * E_f, final_dir.GetZ() * E_f);
-    vP.SetMomentum(corsika::stack::MomentumVector(vP_dir.GetCoordinateSystem(), vec));
-  }
-
-  template <>
-  EProcessReturn ContinuousProcess::DoContinuous(setup::Stack::ParticleType& vP,
-                                                 setup::Trajectory const& vT) {
-    using namespace corsika::units::si; // required for operator::_MeV
-
-    if (!CanInteract(vP.GetPID())) return process::EProcessReturn::eOk;
-    if (vT.GetLength() == 0_m) return process::EProcessReturn::eOk;
-
-    // calculate passed grammage
-    auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
-
-    // Get or build corresponding track integral calculator and solve the
-    // integral
-    auto c = GetCalculator(vP, calc);
-    auto final_energy = get<eDISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
-                            vP.GetEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
-                        1_MeV;
-    auto dE = vP.GetEnergy() - final_energy;
-    energy_lost_ += dE;
-
-    // if the particle has a charge take multiple scattering into account
-    if (vP.GetChargeNumber() != 0) Scatter(vP, dE, dX);
-    vP.SetEnergy(final_energy);
-    vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
-    return process::EProcessReturn::eOk;
-  }
-
-  template <>
-  units::si::LengthType ContinuousProcess::MaxStepLength(
-      setup::Stack::ParticleType const& vP, setup::Trajectory const& vT) {
-    using namespace corsika::units::si; // required for operator::_MeV
-
-    if (!CanInteract(vP.GetPID()))
-      return units::si::meter * std::numeric_limits<double>::infinity();
-
-    // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
-    // hyper parameter which must be adjusted.
-    //
-    // in any case: never go below 0.99*emCut_ This needs to be
-    // slightly smaller than emCut_ since, either this Step is limited
-    // by energy_lim, then the particle is stopped in a very short
-    // range (before doing anythin else) and is then removed
-    // instantly. The exact position where it reaches emCut is not
-    // important, the important fact is that its E_kin is zero
-    // afterwards.
-    //
-    auto energy_lim = std::max(0.9 * vP.GetEnergy(), 0.99 * emCut_);
-
-    // solving the track integral for giving energy lim
-    auto c = GetCalculator(vP, calc);
-    auto grammage = get<eDISPLACEMENT>(c->second)->SolveTrackIntegral(
-                        vP.GetEnergy() / 1_MeV, energy_lim / 1_MeV) *
-                    1_g / square(1_cm);
-
-    // return it in distance aequivalent
-    auto dist = vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage);
-    C8LOG_TRACE("PROPOSAL::MaxStepLength X={} g/cm2, l={} m ",
-                grammage / 1_g * square(1_cm), dist / 1_m);
-    return dist;
-  }
-
-  void ContinuousProcess::showResults() const {
-    using namespace corsika::units::si; // required for operator::_MeV
-    std::cout << " ******************************" << std::endl
-              << " PROCESS::ContinuousProcess: " << std::endl;
-    std::cout << " energy lost dE (GeV)      :  " << energy_lost_ / 1_GeV << std::endl;
-  }
-
-  void ContinuousProcess::reset() {
-    using namespace corsika::units::si; // required for operator::_MeV
-    energy_lost_ = 0_GeV;
-  }
-
-} // namespace corsika::process::proposal

Processes/Proposal/Interaction.cc

-/*
- * (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.
- */
-
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/process/proposal/Interaction.h>
-#include <corsika/setup/SetupEnvironment.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/COMBoost.h>
-#include <limits>
-#include <memory>
-#include <random>
-#include <tuple>
-
-namespace corsika::process::proposal {
-
-  template <>
-  Interaction::Interaction(setup::Environment const& _env,
-                           corsika::units::si::HEPEnergyType _emCut)
-      : ProposalProcessBase(_env, _emCut) {}
-
-  void Interaction::BuildCalculator(particles::Code code,
-                                    environment::NuclearComposition const& comp) {
-    // search crosssection builder for given particle
-    auto p_cross = cross.find(code);
-    if (p_cross == cross.end())
-      throw std::runtime_error("PROPOSAL could not find corresponding builder");
-
-    // interpolate the crosssection for given media and energy cut. These may
-    // take some minutes if you have to build the tables and cannot read the
-    // from disk
-    auto c = p_cross->second(media.at(comp.hash()), emCut_);
-
-    // Look which interactions take place and build the corresponding
-    // interaction and secondarie builder. The interaction integral will
-    // interpolated too and saved in the calc map by a key build out of a hash
-    // of composed of the component and particle code.
-    auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
-    calc[std::make_pair(comp.hash(), code)] = std::make_tuple(
-        PROPOSAL::make_secondaries(inter_types, particle[code], media.at(comp.hash())),
-        PROPOSAL::make_interaction(c, true));
-  }
-
-  template <>
-  corsika::process::EProcessReturn Interaction::DoInteraction(setup::StackView& view) {
-    using namespace corsika::units::si; // required for operator::_MeV
-
-    auto const projectile = view.GetProjectile();
-
-    if (CanInteract(projectile.GetPID())) {
-      // Get or build corresponding calculators
-      auto c = GetCalculator(projectile, calc);
-
-      // Get the rates of the interaction types for every component.
-      std::uniform_real_distribution<double> distr(0., 1.);
-
-      // sample a interaction-type, loss and component
-      auto rates = get<eINTERACTION>(c->second)->Rates(projectile.GetEnergy() / 1_MeV);
-      auto [type, comp_ptr, v] = get<eINTERACTION>(c->second)->SampleLoss(
-          projectile.GetEnergy() / 1_MeV, rates, distr(fRNG));
-
-      // Read how much random numbers are required to calculate the secondaries.
-      // Calculate the secondaries and deploy them on the corsika stack.
-      auto rnd =
-          vector<double>(get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type));
-      for (auto& it : rnd) it = distr(fRNG);
-      auto point = PROPOSAL::Vector3D(projectile.GetPosition().GetX() / 1_cm,
-                                      projectile.GetPosition().GetY() / 1_cm,
-                                      projectile.GetPosition().GetZ() / 1_cm);
-      auto projectile_dir = projectile.GetDirection();
-      auto d = projectile_dir.GetComponents();
-      auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
-                                          d.GetZ().magnitude());
-      auto loss = make_tuple(static_cast<int>(type), point, direction,
-                             v * projectile.GetEnergy() / 1_MeV, 0.);
-      auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries(
-          projectile.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd);
-      for (auto& s : sec) {
-        auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
-        auto vec = corsika::geometry::QuantityVector(
-            get<PROPOSAL::Loss::DIRECTION>(s).GetX() * E,
-            get<PROPOSAL::Loss::DIRECTION>(s).GetY() * E,
-            get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * E);
-        auto p =
-            corsika::stack::MomentumVector(projectile_dir.GetCoordinateSystem(), vec);
-        auto sec_code = corsika::particles::ConvertFromPDG(
-            static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
-        view.AddSecondary(
-            make_tuple(sec_code, E, p, projectile.GetPosition(), projectile.GetTime()));
-      }
-    }
-    return process::EProcessReturn::eOk;
-  }
-
-  template <>
-  corsika::units::si::GrammageType Interaction::GetInteractionLength(
-      setup::Stack::StackIterator const& projectile) {
-    using namespace corsika::units::si; // required for operator::_MeV
-
-    if (CanInteract(projectile.GetPID())) {
-      auto c = GetCalculator(projectile, calc);
-      return get<eINTERACTION>(c->second)->MeanFreePath(projectile.GetEnergy() / 1_MeV) *
-             1_g / (1_cm * 1_cm);
-    }
-    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-  }
-} // namespace corsika::process::proposal

Processes/Proposal/Interaction.h

-/*
- * (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 <PROPOSAL/PROPOSAL.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/InteractionProcess.h>
-#include <corsika/process/proposal/ProposalProcessBase.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/random/UniformRealDistribution.h>
-#include <array>
-
-namespace corsika::process::proposal {
-
-  //!
-  //! Electro-magnetic and gamma stochastic losses produced by proposal. It makes
-  //! use of interpolation tables which are runtime intensive calculation, but can be
-  //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
-  //!
-  class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
-
-    enum { eSECONDARIES, eINTERACTION };
-    using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
-                               unique_ptr<PROPOSAL::Interaction>>;
-
-    std::unordered_map<calc_key_t, calculator_t, hash>
-        calc; //!< Stores the secondaries and interaction calculators.
-
-    //!
-    //! Build the secondaries and interaction calculators and add it to calc.
-    //!
-    void BuildCalculator(particles::Code, environment::NuclearComposition const&) final;
-
-  public:
-    //!
-    //! Produces the stoachastic loss calculator for leptons based on nuclear
-    //! compositions and stochastic description limited by the particle cut.
-    //!
-    template <typename TEnvironment>
-    Interaction(TEnvironment const& env, corsika::units::si::HEPEnergyType emCut);
-
-    //!
-    //! Calculate the rates for the different targets and interactions. Sample a
-    //! pair of interaction-type, component and rate, followed by sampling a loss and
-    //! produce the corresponding secondaries and store them on the particle stack.
-    //!
-    template <typename TSecondaryView>
-    corsika::process::EProcessReturn DoInteraction(TSecondaryView&);
-
-    //!
-    //! Calculates the  mean free path length
-    //!
-    template <typename TParticle>
-    corsika::units::si::GrammageType GetInteractionLength(TParticle const& p);
-  };
-} // namespace corsika::process::proposal

PROPOSAL @ 92b9b9ca

-Subproject commit 92b9b9ca20826725cb805d135a1a5d5b29a7e06a

Processes/Proposal/ProposalProcessBase.cc

-/*
- * (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.
- */
-
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/process/proposal/ProposalProcessBase.h>
-#include <corsika/setup/SetupEnvironment.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/COMBoost.h>
-#include <cstdlib>
-#include <iostream>
-#include <limits>
-#include <memory>
-#include <random>
-#include <tuple>
-
-namespace corsika::process::proposal {
-  bool ProposalProcessBase::CanInteract(particles::Code pcode) const {
-    if (std::find(begin(tracked), end(tracked), pcode) != end(tracked)) return true;
-    return false;
-  }
-
-  ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env,
-                                           corsika::units::si::HEPEnergyType _emCut)
-      : emCut_(_emCut)
-      , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
-    using namespace corsika::units::si; // required for operator::_MeV
-    _env.GetUniverse()->walk([&](auto& vtn) {
-      if (vtn.HasModelProperties()) {
-        const auto& prop = vtn.GetModelProperties();
-        const auto& medium = mediumData(prop.medium(corsika::geometry::Point(
-            geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(), 0_cm,
-            0_cm, 0_cm)));
-
-        auto comp_vec = std::vector<PROPOSAL::Components::Component>();
-        const auto& comp = prop.GetNuclearComposition();
-        auto frac_iter = comp.GetFractions().cbegin();
-        for (auto& pcode : comp.GetComponents()) {
-          comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
-                                *frac_iter);
-          ++frac_iter;
-        }
-
-        media[comp.hash()] =
-            PROPOSAL::Medium(medium.name(), medium.Ieff(), -medium.Cbar(), medium.aa(),
-                             medium.sk(), medium.x0(), medium.x1(), medium.dlt0(),
-                             medium.corrected_density(), comp_vec);
-      }
-    });
-
-    PROPOSAL::InterpolationDef::order_of_interpolation = 2;
-    PROPOSAL::InterpolationDef::nodes_cross_section = 100;
-    PROPOSAL::InterpolationDef::nodes_propagate = 1000;
-
-    //! If corsika data exist store interpolation tables to the corresponding
-    //! path, otherwise interpolation tables would only stored in main memory if
-    //! no explicit intrpolation def is specified.
-    if (auto data_path = std::getenv("CORSIKA_DATA")) {
-      PROPOSAL::InterpolationDef::path_to_tables = std::string(data_path) + "/PROPOSAL";
-    } else {
-      throw std::runtime_error(
-          "It is not recommended to run PROPOSAL without its tables in "
-          "$CORSIKA_DATA/PROPOSAL. This would be extremely slow. Please provide the "
-          "table directory. ");
-    }
-  }
-
-  size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept {
-    return p.first ^ std::hash<particles::Code>{}(p.second);
-  }
-
-} // namespace corsika::process::proposal

Processes/Proposal/ProposalProcessBase.h

-/*
- * (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 <PROPOSAL/PROPOSAL.h>
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/setup/SetupEnvironment.h>
-#include <array>
-
-namespace corsika::process::proposal {
-
-  //!
-  //! Particles which can be handled by proposal. That means they can be
-  //! propagated and decayed if they decays.
-  //!
-  static constexpr std::array<particles::Code, 7> tracked{
-      particles::Code::Gamma,    particles::Code::Electron, particles::Code::Positron,
-      particles::Code::MuMinus,  particles::Code::MuPlus,   particles::Code::TauPlus,
-      particles::Code::TauMinus,
-  };
-
-  //!
-  //! Internal map from particle codes to particle properties required for
-  //! crosssections, decay and scattering algorithms. In the future the
-  //! particles may be created by reading out the Corsica constants.
-  //!
-  static std::map<particles::Code, PROPOSAL::ParticleDef> particle = {
-      {particles::Code::Gamma, PROPOSAL::GammaDef()},
-      {particles::Code::Electron, PROPOSAL::EMinusDef()},
-      {particles::Code::Positron, PROPOSAL::EPlusDef()},
-      {particles::Code::MuMinus, PROPOSAL::MuMinusDef()},
-      {particles::Code::MuPlus, PROPOSAL::MuPlusDef()},
-      {particles::Code::TauMinus, PROPOSAL::TauMinusDef()},
-      {particles::Code::TauPlus, PROPOSAL::TauPlusDef()}};
-
-  //!
-  //! Crosssection factories for different particle types.
-  //!
-  template <typename T>
-  static auto cross_builder =
-      [](PROPOSAL::Medium& m, corsika::units::si::HEPEnergyType emCut) {
-        using namespace corsika::units::si;
-        auto p_cut =
-            std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, true);
-        return PROPOSAL::DefaultCrossSections<T>::template Get<std::false_type>(
-            T(), m, p_cut, true);
-      };
-
-  //!
-  //! PROPOSAL default crosssections are maped to corresponding corsika particle
-  //! code.
-  //!
-  static std::map<particles::Code,
-                  std::function<PROPOSAL::crosssection_list_t<PROPOSAL::ParticleDef,
-                                                              PROPOSAL::Medium>(
-                      PROPOSAL::Medium&, corsika::units::si::HEPEnergyType)>>
-      cross = {{particles::Code::Gamma, cross_builder<PROPOSAL::GammaDef>},
-               {particles::Code::Electron, cross_builder<PROPOSAL::EMinusDef>},
-               {particles::Code::Positron, cross_builder<PROPOSAL::EPlusDef>},
-               {particles::Code::MuMinus, cross_builder<PROPOSAL::MuMinusDef>},
-               {particles::Code::MuPlus, cross_builder<PROPOSAL::MuPlusDef>},
-               {particles::Code::TauMinus, cross_builder<PROPOSAL::TauMinusDef>},
-               {particles::Code::TauPlus, cross_builder<PROPOSAL::TauPlusDef>}};
-
-  //!
-  //! PROPOSAL base process which handels mapping of particle codes to
-  //! stored interpolation tables.
-  //!
-  class ProposalProcessBase {
-  protected:
-    corsika::units::si::HEPEnergyType
-        emCut_;                 //!< Stochastic losses smaller than the given cut
-                                //!< will be handeled continuously.
-    corsika::random::RNG& fRNG; //!< random number generator used by proposal
-
-    std::unordered_map<std::size_t, PROPOSAL::Medium>
-        media; //!< maps nuclear composition from univers to media to produce
-               //!< crosssections, which requires further ionization constants.
-
-    //!
-    //! Store cut and  nuclear composition of the whole universe in media which are
-    //! required for creating crosssections by proposal.
-    //!
-    ProposalProcessBase(corsika::setup::Environment const& _env,
-                        corsika::units::si::HEPEnergyType _emCut);
-
-    //!
-    //! Checks if a particle can be processed by proposal
-    //!
-    bool CanInteract(particles::Code pcode) const;
-
-    using calc_key_t = std::pair<std::size_t, particles::Code>;
-
-    //!
-    //! Hash to store interpolation tables related to a pair of particle and nuclear
-    //! composition.
-    //!
-    struct hash {
-      size_t operator()(const calc_key_t& p) const noexcept;
-    };
-
-    //!
-    //! Builds the calculator to the corresponding class
-    //!
-    virtual void BuildCalculator(particles::Code,
-                                 environment::NuclearComposition const&) = 0;
-
-    //!
-    //! Searches the particle dependet calculator dependent of actuall medium composition
-    //! and particle type. If no calculator is found, the corresponding new calculator is
-    //! built and then returned.
-    //!
-    template <typename Particle, typename Calculators>
-    auto GetCalculator(Particle& vP, Calculators& calc) {
-      const auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
-      auto calc_it = calc.find(std::make_pair(comp.hash(), vP.GetPID()));
-      if (calc_it != calc.end()) return calc_it;
-      BuildCalculator(vP.GetPID(), comp);
-      return GetCalculator(vP, calc);
-    }
-  };
-} // namespace corsika::process::proposal

Processes/Pythia/testPythia8.cc

-/*
- * (c) Copyright 2018 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 <Pythia8/Pythia.h>
-#include <corsika/process/pythia/Decay.h>
-#include <corsika/process/pythia/Interaction.h>
-
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/geometry/Point.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <catch2/catch.hpp>
-
-TEST_CASE("Pythia", "[processes]") {
-
-  SECTION("linking pythia") {
-    using namespace Pythia8;
-    using std::cout;
-    using std::endl;
-
-    // Generator. Process selection. LHC initialization. Histogram.
-    Pythia pythia;
-
-    pythia.readString("Next:numberShowInfo = 0");
-    pythia.readString("Next:numberShowProcess = 0");
-    pythia.readString("Next:numberShowEvent = 0");
-
-    pythia.readString("ProcessLevel:all = off");
-
-    pythia.init();
-
-    Event& event = pythia.event;
-    event.reset();
-
-    pythia.particleData.mayDecay(321, true);
-    double pz = 100.;
-    double m = 0.49368;
-    event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
-
-    if (!pythia.next())
-      cout << "decay failed!" << endl;
-    else
-      cout << "particles after decay: " << event.size() << endl;
-    event.list();
-
-    // loop over final state
-    for (int i = 0; i < pythia.event.size(); ++i)
-      if (pythia.event[i].isFinal()) {
-        cout << "particle: id=" << pythia.event[i].id() << endl;
-      }
-  }
-
-  SECTION("pythia interface") {
-    using namespace corsika;
-
-    random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-
-    process::pythia::Decay model;
-  }
-}
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/setup/SetupEnvironment.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-using namespace corsika;
-using namespace corsika::units::si;
-
-template <typename TStackView>
-auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
-  geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
-
-  for (auto const& p : view) { sum += p.GetMomentum(); }
-
-  return sum;
-}
-
-TEST_CASE("pythia process") {
-
-  auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Proton);
-  auto const& cs = *csPtr;
-  [[maybe_unused]] auto const& env_dummy = env;
-  [[maybe_unused]] auto const& node_dummy = nodePtr;
-
-  SECTION("pythia decay") {
-    feenableexcept(FE_INVALID);
-    const HEPEnergyType P0 = 10_GeV;
-    auto [stackPtr, secViewPtr] =
-        setup::testing::setupStack(particles::Code::PiPlus, 0, 0, P0, nodePtr, *csPtr);
-    const auto plab = corsika::stack::MomentumVector(
-        cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
-    auto& stack = *stackPtr;
-    auto& view = *secViewPtr;
-    auto particle = stackPtr->first();
-
-    random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-
-    process::pythia::Decay model;
-
-    [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
-    model.DoDecay(view);
-    CHECK(stack.getEntries() == 3);
-    auto const pSum = sumMomentum(view, cs);
-    CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(1e-4));
-    CHECK((pSum.norm() - plab.norm()) / 1_GeV == Approx(0).margin(1e-4));
-  }
-
-  SECTION("pythia decay config") {
-    process::pythia::Decay model({particles::Code::PiPlus, particles::Code::PiMinus});
-    REQUIRE(model.IsDecayHandled(particles::Code::PiPlus));
-    REQUIRE(model.IsDecayHandled(particles::Code::PiMinus));
-    REQUIRE_FALSE(model.IsDecayHandled(particles::Code::KPlus));
-
-    const std::vector<particles::Code> particleTestList = {
-        particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-        particles::Code::Lambda0Bar, particles::Code::D0Bar};
-
-    // setup decays
-    model.SetHandleDecay(particleTestList);
-    for (auto& pCode : particleTestList) REQUIRE(model.IsDecayHandled(pCode));
-
-    // individually
-    model.SetHandleDecay(particles::Code::KMinus);
-
-    // possible decays
-    REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Proton));
-    REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Electron));
-    REQUIRE(model.CanHandleDecay(particles::Code::PiPlus));
-    REQUIRE(model.CanHandleDecay(particles::Code::MuPlus));
-  }
-
-  SECTION("pythia interaction") {
-
-    feenableexcept(FE_INVALID);
-    auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::PiPlus, 0,
-                                                             0, 100_GeV, nodePtr, *csPtr);
-    auto& view = *secViewPtr;
-    auto particle = stackPtr->first();
-
-    process::pythia::Interaction model;
-
-    [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
-    [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
-  }
-}

Processes/QGSJetII/code_generator.py

-#!/usr/bin/env python3
-
-# (c) Copyright 2020 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)
-    definition of particle_db dict is: "name", "antiName", "pdg", "mass", "electric_charge", "lifetime", "ngc_code", "isNucleus", "isHadron"
-    '''
-    with open(filename, "rb") as f:
-        particle_db = pickle.load(f)
-    return particle_db
-
-
-def set_default_qgsjetII_definition(particle_db):
-    '''
-    Also particles not explicitly known by QGSJetII may in fact interact via mapping 
-    to cross section types (xsType) and hadron type (hadronType)
-
-    This is achieved here.
-
-    The function return nothing, but modified the input particle_db by adding the 
-    fields 'xsType' and 'hadronType'
-    '''
-    for identifier, pData in particle_db.items():
-        # the cross-section types
-        xsType = "CannotInteract"
-        hadronType = "UndefinedType"
-        if (pData['isNucleus']):
-            xsType = "Baryons"
-            hadronType = "NucleusType"
-        elif (pData['isHadron']):
-            pdg = abs(pData['pdg'])
-            anti = pData['pdg'] < 0
-            isBaryon = (1000 <= pdg < 4000)
-            charge = pData['electric_charge']
-            if (pdg>=100 and pdg<300 and pdg!=130): # light mesons
-                xsType = "LightMesons"
-                if (charge==0):
-                    hadronType = "NeutralLightMesonType"
-                else:
-                    if (charge>0):
-                        hadronType = "PiPlusType"
-                    else:
-                        hadronType = "PiMinusType"               
-            elif ((pdg>=300 and pdg<400) or pdg in [130, 10313, 10323]): # kaons
-                xsType = "Kaons"
-                if (charge>0):
-                    hadronType = "KaonPlusType"
-                else:
-                    hadronType = "KaonMinusType"
-                if (charge==0):
-                    hadronType = "Kaon0SType"
-                    if (pdg == 130):
-                        hadronType = "Kaon0LType"
-                    elif (pdg == 310):
-                        hadronType = "Kaon0SType"
-            elif (isBaryon or pData['isNucleus']): # baryons
-                xsType = "Baryons"
-                if (charge==0):
-                    if (anti):
-                        hadronType = "AntiNeutronType"
-                    else: 
-                        hadronType = "NeutronType"
-                else:
-                    if (charge>0):
-                        hadronType = "ProtonType"
-                    else:
-                        hadronType = "AntiProtonType"
-            # all othe not-captured cased are hopefully irrelevant
-            
-        pData['qgsjetII_xsType'] = xsType
-        pData['qgsjetII_hadronType'] = hadronType
-
-            
-def read_qgsjetII_codes(filename, particle_db):
-    '''
-    reads the qgsjet-codes data file. For particles known to QGSJetII the 'qgsjetII_code' is set in the particle_db, as
-    well as the 'xsType' is updated in case it is different from its default value set above. 
-    '''
-    with open(filename) as f:
-        for line in f:
-            line = line.strip()
-            if len(line)==0 or line[0] == '#':
-                continue
-            line = line.split('#')[0]
-            print (line)
-            identifier, model_code, xsType = line.split()
-            try:
-                particle_db[identifier]["qgsjetII_code"] = int(model_code)
-                particle_db[identifier]["qgsjetII_xsType"] = xsType
-            except KeyError as e:
-                raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
-
-            
-def generate_qgsjetII_enum(particle_db):
-    '''
-    generates the enum to access qgsjetII particles by readable names
-    '''
-    output = "enum class QgsjetIICode : int8_t {\n"
-    for identifier, pData in particle_db.items():
-        if 'qgsjetII_code' in pData:
-            output += "  {:s} = {:d},\n".format(identifier, pData['qgsjetII_code'])
-    output += "};\n"
-    return output
-
-
-def generate_corsika2qgsjetII(particle_db):    
-    '''
-    generates the look-up table to convert corsika codes to qgsjetII codes
-    '''
-    string = "std::array<QgsjetIICode, {:d}> constexpr corsika2qgsjetII = {{\n".format(len(particle_db))
-    for identifier, pData in particle_db.items():
-        if 'qgsjetII_code' in pData:
-            string += "  QgsjetIICode::{:s}, \n".format(identifier)
-        else:
-            string += "  QgsjetIICode::Unknown, // {:s}\n".format(identifier + ' not implemented in QGSJetII')
-    string += "};\n"
-    return string
-    
-
-def generate_corsika2qgsjetII_xsType(particle_db):    
-    '''
-    generates the look-up table to convert corsika codes to qgsjetII codes
-    '''
-    string = "std::array<QgsjetIIXSClass, {:d}> constexpr corsika2qgsjetIIXStype = {{\n".format(len(particle_db))
-    for identifier, pData in particle_db.items():
-        modelCodeXS = pData.get("qgsjetII_xsType", "CannotInteract")
-        string += "  QgsjetIIXSClass::{:s}, // {:s}\n".format(modelCodeXS, identifier if modelCodeXS else identifier + " (not implemented in QGSJETII)")
-    string += "};\n"
-    return string
-
-
-def generate_corsika2qgsjetII_hadronType(particle_db):    
-    '''
-    generates the look-up table to convert corsika codes to qgsjetII codes
-    '''
-    string = "std::array<QgsjetIIHadronType, {:d}> constexpr corsika2qgsjetIIHadronType = {{\n".format(len(particle_db))
-    for identifier, pData in particle_db.items():
-        modelCode = pData.get("qgsjetII_hadronType", "UndefinedType")
-        string += "  QgsjetIIHadronType::{:s}, // {:s}\n".format(modelCode, identifier if modelCode else identifier + " (not implemented in QGSJETII)")
-    string += "};\n"
-    return string
-
-
-def generate_qgsjetII2corsika(particle_db) :
-    '''
-    generates the look-up table to convert qgsjetII codes to corsika codes    
-    '''
-    minID = 0
-    for identifier, pData in particle_db.items() :
-        if 'qgsjetII_code' in pData:
-            minID = min(minID, pData['qgsjetII_code'])
-
-    string = "QgsjetIICodeIntType constexpr minQgsjetII = {:d};\n\n".format(minID)
-
-    pDict = {}
-    for identifier, pData in particle_db.items() :
-        if 'qgsjetII_code' in pData:
-            model_code = pData['qgsjetII_code'] - minID
-            pDict[model_code] = identifier
-    
-    nPart = max(pDict.keys()) - min(pDict.keys()) + 1
-    string += "std::array<corsika::particles::Code, {:d}> constexpr qgsjetII2corsika = {{\n".format(nPart)
-    
-    for iPart in range(nPart) :
-        identifier = pDict.get(iPart, "Unknown")
-        qgsID = iPart + minID
-        string += "  corsika::particles::Code::{:s}, // {:d} \n".format(identifier, qgsID)
-    
-    string += "};\n"
-    return string
-
-
-if __name__ == "__main__":
-    if len(sys.argv) != 3:
-        print("usage: {:s} <particle_db.pkl> <qgsjetII_codes.dat>".format(sys.argv[0]), file=sys.stderr)
-        sys.exit(1)
-        
-    print("code_generator.py for QGSJETII")
-    
-    particle_db = load_particledb(sys.argv[1])
-    read_qgsjetII_codes(sys.argv[2], particle_db)
-    set_default_qgsjetII_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_qgsjetII_enum(particle_db), file=f)
-        print(generate_corsika2qgsjetII(particle_db), file=f)
-        print(generate_qgsjetII2corsika(particle_db), file=f)
-        print(generate_corsika2qgsjetII_xsType(particle_db), file=f)
-        print(generate_corsika2qgsjetII_hadronType(particle_db), file=f)

Processes/Sibyll/sibyll2.3d.cc

-/*
- * (c) Copyright 2018 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/process/sibyll/sibyll2.3d.h>
-
-#include <corsika/random/RNGManager.h>
-#include <random>
-
-int get_nwounded() { return s_chist_.nwd; }
-double get_sibyll_mass2(int& id) { return s_mass1_.am2[abs(id) - 1]; }
-
-double s_rndm_(int&) {
-  static corsika::random::RNG& rng =
-      corsika::random::RNGManager::GetInstance().GetRandomStream("sibyll");
-
-  std::uniform_real_distribution<double> dist;
-  return dist(rng);
-}

Processes/Sibyll/sibyll2.3d.h

-/*
- * (c) Copyright 2018 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
-//----------------------------------------------
-//  C++ interface for the SIBYLL event generator
-//----------------------------------------------
-// wrapper
-
-extern "C" {
-
-typedef char s_name[6];
-
-// SIBYLL particle stack (FORTRAN COMMON)
-// variables are: np : numer of particles on stack
-//                 p : 4momentum + mass of particles on stack
-//             llist : id of particles on stack
-extern struct {
-  double p[5][8000];
-  int llist[8000];
-  int np;
-} s_plist_;
-
-// additional information about interactions.
-// number of wounded nucleons, number of hard and soft scatterings etc.
-extern struct { int nnsof[20], nnjet[20], jdif[20], nwd, njet, nsof; } s_chist_;
-
-extern struct {
-  double cbr[223 + 16 + 12 + 8];
-  int kdec[1338 + 6 * (16 + 12 + 8)];
-  int lbarp[99];
-  int idb[99];
-} s_csydec_;
-
-// additional particle stack for the mother particles of unstable particles
-// stable particles have entry zero
-extern struct { int llist1[8000]; } s_plist1_;
-
-// tables with particle properties
-// charge, strangeness and baryon number
-extern struct {
-  int ichp[99];
-  int istr[99];
-  int ibar[99];
-} s_chp_;
-
-// tables with particle properties
-// mass and mass squared
-extern struct {
-  double am[99];
-  double am2[99];
-} s_mass1_;
-
-// table with particle names
-extern struct { char namp[6][99]; } s_cnam_;
-
-// debug info
-extern struct {
-  int ncall;
-  int ndebug;
-  int lun;
-} s_debug_;
-
-// lund random generator setup
-// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
-
-// sibyll main subroutine
-void sibyll_(const int&, const int&, const double&);
-
-// subroutine to initiate sibyll
-void sibyll_ini_();
-
-// subroutine to SET DECAYS
-void dec_ini_();
-
-// subroutine to initiate random number generator
-// void rnd_ini_();
-
-// print event
-void sib_list_(int&);
-
-// decay routine
-void decsib_();
-
-// interaction length
-// double fpni_(double&, int&);
-
-void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&, double&);
-void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
-                   double&);
-
-double s_rndm_(int&);
-
-int get_nwounded();
-double get_sibyll_mass2(int&);
-
-// phojet random generator setup
-void pho_rndin_(int&, int&, int&, int&);
-}