Processes/ExampleProcessors/DummySecondariesProcess.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
-
-#include <corsika/process/SecondariesProcess.h>
-
-#include <chrono>
-#include <thread>
-
-namespace corsika::process {
-  namespace example_processors {
-
-    template <int ISleep>
-    class DummySecondariesProcess
-        : public SecondariesProcess<DummySecondariesProcess<ISleep>> {
-    private:
-    public:
-      template <typename TSecondaries>
-      inline EProcessReturn DoSecondaries(TSecondaries&) {
-        std::this_thread::sleep_for(std::chrono::milliseconds(ISleep));
-        return process::EProcessReturn::eOk;
-      }
-    };
-
-  } // namespace example_processors
-} // namespace corsika::process
\ No newline at end of file

Processes/ExampleProcessors/TestDummy.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/example_processors/DummyBoundaryCrossingProcess.h>
-#include <corsika/process/example_processors/DummyContinuousProcess.h>
-#include <corsika/process/example_processors/DummyDecayProcess.h>
-#include <corsika/process/example_processors/DummyInteractionProcess.h>
-#include <corsika/process/example_processors/DummySecondariesProcess.h>
-
-#include <corsika/process/ProcessReturn.h>
-
-#include <catch2/catch.hpp>
-
-#include <chrono>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::process::example_processors;
-
-using namespace corsika::units::si;
-
-TEST_CASE("Dummy Processes") {
-  DummyBoundaryCrossingProcess<1000> dbc;
-  DummyContinuousProcess<1000> dc;
-  DummyDecayProcess<1000> dd;
-  DummyInteractionProcess<1000> di;
-  DummySecondariesProcess<1000> dse;
-
-  int tmp = 0;
-  SECTION("BoundaryCrossing") {
-    auto start = std::chrono::steady_clock::now();
-    REQUIRE(dbc.DoBoundaryCrossing(tmp, 0, 0) == EProcessReturn::eOk);
-    auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-  }
-
-  SECTION("Continuous") {
-    auto start = std::chrono::steady_clock::now();
-    REQUIRE(dc.DoContinuous(tmp, nullptr) == EProcessReturn::eOk);
-    auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-
-    start = std::chrono::steady_clock::now();
-    REQUIRE(dc.MaxStepLength(nullptr, nullptr) ==
-            units::si::meter * std::numeric_limits<double>::infinity());
-    end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-  }
-
-  SECTION("Decay") {
-    auto start = std::chrono::steady_clock::now();
-    REQUIRE(dd.DoDecay(tmp) == EProcessReturn::eOk);
-    auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-
-    start = std::chrono::steady_clock::now();
-    REQUIRE(dd.GetLifetime(tmp) ==
-            units::si::second * std::numeric_limits<double>::infinity());
-    end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-  }
-
-  SECTION("Interaction") {
-    auto start = std::chrono::steady_clock::now();
-    REQUIRE(di.DoInteraction(tmp) == EProcessReturn::eOk);
-    auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-
-    start = std::chrono::steady_clock::now();
-    REQUIRE(di.GetInteractionLength(tmp) ==
-            (units::si::gram / 1_cm / 1_cm) * std::numeric_limits<double>::infinity());
-    end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-  }
-
-  SECTION("Secondaries") {
-    auto start = std::chrono::steady_clock::now();
-    REQUIRE(dse.DoSecondaries(tmp) == EProcessReturn::eOk);
-    auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(1000).margin(1));
-  }
-}

Processes/InteractionCounter/CMakeLists.txt

-set (
-  MODEL_HEADERS
-  InteractionCounter.hpp
-  InteractionHistogram.hpp
-  )
-
-set (
-  MODEL_SOURCES
-  InteractionHistogram.cc
-)
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/interaction_counter
-  )
-
-add_library (ProcessInteractionCounter STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessInteractionCounter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessInteractionCounter
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessInteractionCounter
-  CORSIKAunits
-  CORSIKAutilities
-  CORSIKAprocesssequence
-  CORSIKAthirdparty
-  C8::ext::boost
-  )
-
-target_include_directories (
-  ProcessInteractionCounter 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST(testInteractionCounter SOURCES
-  testInteractionCounter.cc
-  ${MODEL_HEADERS}
-)
-#
-target_compile_definitions (
-  testInteractionCounter
-  PRIVATE
-  REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}"
-  )
-#
-target_link_libraries (
-  testInteractionCounter
-  ProcessInteractionCounter
-  CORSIKAsetup
-  CORSIKAunits
-  CORSIKAenvironment
-  CORSIKAtesting
-  )

Processes/InteractionCounter/InteractionCounter.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/process/interaction_counter/InteractionHistogram.hpp>
-
-#include <corsika/process/InteractionProcess.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/setup/SetupStack.h>
-
-namespace corsika::process::interaction_counter {
-  /*!
-   * Wrapper around an InteractionProcess that fills histograms of the number
-   * of calls to DoInteraction() binned in projectile energy (both in
-   * lab and center-of-mass frame) and species
-   */
-  template <class TCountedProcess>
-  class InteractionCounter
-      : public InteractionProcess<InteractionCounter<TCountedProcess>> {
-
-    TCountedProcess& process_;
-    InteractionHistogram histogram_;
-
-  public:
-    InteractionCounter(TCountedProcess& process)
-        : process_(process) {}
-
-    template <typename TSecondaryView>
-    auto DoInteraction(TSecondaryView& view) {
-      auto const projectile = view.GetProjectile();
-      auto const massNumber = projectile.GetNode()
-                                  ->GetModelProperties()
-                                  .GetNuclearComposition()
-                                  .GetAverageMassNumber();
-      auto const massTarget = massNumber * units::constants::nucleonMass;
-
-      if (auto const projectile_id = projectile.GetPID();
-          projectile_id == particles::Code::Nucleus) {
-        auto const A = projectile.GetNuclearA();
-        auto const Z = projectile.GetNuclearZ();
-        histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget, A, Z);
-      } else {
-        histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget);
-      }
-      return process_.DoInteraction(view);
-    }
-
-    template <typename TParticle>
-    auto GetInteractionLength(TParticle const& particle) const {
-      return process_.GetInteractionLength(particle);
-    }
-
-    InteractionHistogram const& GetHistogram() const { return histogram_; }
-  };
-
-} // namespace corsika::process::interaction_counter

Processes/InteractionCounter/InteractionHistogram.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/interaction_counter/InteractionHistogram.hpp>
-
-#include <fstream>
-#include <string>
-
-using namespace corsika::process::interaction_counter;
-
-InteractionHistogram::InteractionHistogram()
-    : inthist_cms_{detail::hist_factory(num_bins_cms, lower_edge_cms, upper_edge_cms)}
-    , inthist_lab_{detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab)} {}
-
-void InteractionHistogram::fill(particles::Code projectile_id,
-                                units::si::HEPEnergyType lab_energy,
-                                units::si::HEPEnergyType mass_target, int A, int Z) {
-  using namespace units::si;
-  auto constexpr inv_eV = 1 / 1_eV;
-  if (projectile_id == particles::Code::Nucleus) {
-    auto const sqrtS =
-        sqrt(A * A * (units::constants::nucleonMass * units::constants::nucleonMass) +
-             mass_target * mass_target + 2 * lab_energy * mass_target);
-
-    int32_t const pdg = 1'000'000'000l + Z * 10'000l + A * 10l;
-
-    inthist_lab_(pdg, lab_energy * inv_eV);
-    inthist_cms_(pdg, sqrtS * inv_eV);
-  } else {
-    auto const projectile_mass = particles::GetMass(projectile_id);
-    auto const sqrtS = sqrt(projectile_mass * projectile_mass +
-                            mass_target * mass_target + 2 * lab_energy * mass_target);
-
-    inthist_cms_(static_cast<int>(particles::GetPDG(projectile_id)), sqrtS * inv_eV);
-    inthist_lab_(static_cast<int>(particles::GetPDG(projectile_id)), lab_energy * inv_eV);
-  }
-}
-
-void InteractionHistogram::saveLab(std::string const& filename,
-                                   utl::SaveMode mode) const {
-  corsika::utl::save_hist(inthist_lab_, filename, mode);
-}
-
-void InteractionHistogram::saveCMS(std::string const& filename,
-                                   utl::SaveMode mode) const {
-  corsika::utl::save_hist(inthist_cms_, filename, mode);
-}
-
-InteractionHistogram& InteractionHistogram::operator+=(
-    InteractionHistogram const& other) {
-  inthist_lab_ += other.inthist_lab_;
-  inthist_cms_ += other.inthist_cms_;
-
-  return *this;
-}
-
-InteractionHistogram InteractionHistogram::operator+(InteractionHistogram other) const {
-  other.inthist_lab_ += inthist_lab_;
-  other.inthist_cms_ += inthist_cms_;
-
-  return other;
-}

Processes/InteractionCounter/testInteractionCounter.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/interaction_counter/InteractionCounter.hpp>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/setup/SetupStack.h>
-
-#include <catch2/catch.hpp>
-
-#include <numeric>
-
-using namespace corsika;
-using namespace corsika::process::interaction_counter;
-using namespace corsika::units;
-using namespace corsika::units::si;
-
-const std::string refDataDir = std::string(REFDATADIR); // from cmake
-
-struct DummyProcess {
-  template <typename TParticle>
-  GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
-    return 100_g / 1_cm / 1_cm;
-  }
-
-  template <typename TParticle>
-  auto DoInteraction([[maybe_unused]] TParticle& projectile) {
-    return nullptr;
-  }
-};
-
-TEST_CASE("InteractionCounter", "[process]") {
-
-  logging::SetLevel(logging::level::debug);
-
-  DummyProcess d;
-  InteractionCounter countedProcess(d);
-
-  SECTION("GetInteractionLength") {
-    REQUIRE(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
-  }
-
-  auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
-  [[maybe_unused]] auto& env_dummy = env;
-
-  SECTION("DoInteraction nucleus") {
-    unsigned short constexpr A = 14, Z = 7;
-    auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
-                                                             Z, 105_TeV, nodePtr, *csPtr);
-    REQUIRE(stackPtr->getEntries() == 1);
-    REQUIRE(secViewPtr->getEntries() == 0);
-
-    auto const ret = countedProcess.DoInteraction(*secViewPtr);
-    REQUIRE(ret == nullptr);
-
-    auto const& h = countedProcess.GetHistogram().labHist();
-    REQUIRE(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
-    REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
-    auto const& h2 = countedProcess.GetHistogram().CMSHist();
-    REQUIRE(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
-    REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
-
-    countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
-                                          utl::SaveMode::overwrite);
-    countedProcess.GetHistogram().saveCMS("testInteractionCounter_file2.npz",
-                                          utl::SaveMode::overwrite);
-  }
-
-  SECTION("DoInteraction Lambda") {
-    auto constexpr code = particles::Code::Lambda0;
-    auto [stackPtr, secViewPtr] =
-        setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
-    REQUIRE(stackPtr->getEntries() == 1);
-    REQUIRE(secViewPtr->getEntries() == 0);
-
-    auto const ret = countedProcess.DoInteraction(*secViewPtr);
-    REQUIRE(ret == nullptr);
-
-    auto const& h = countedProcess.GetHistogram().labHist();
-    REQUIRE(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
-    REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
-
-    auto const& h2 = countedProcess.GetHistogram().CMSHist();
-    REQUIRE(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
-    REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
-  }
-}
-
-#include <algorithm>
-#include <iterator>
-#include <string>
-#include <fstream>
-
-TEST_CASE("InteractionCounterOutput", "[output validation]") {
-
-  auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
-                       "testInteractionCounter_file2");
-
-  SECTION(std::string("check saved data, ") + file + ".npz") {
-
-    std::cout << file + ".npz vs " << refDataDir + "/" + file + "_REF.npz" << std::endl;
-
-    // compare to binary reference data
-    std::ifstream file1(file + ".npz");
-    std::ifstream file1ref(refDataDir + "/" + file + "_REF.npz");
-
-    std::istreambuf_iterator<char> begin1(file1);
-    std::istreambuf_iterator<char> begin1ref(file1ref);
-
-    std::istreambuf_iterator<char> end;
-
-    while (begin1 != end && begin1ref != end) {
-      CHECK(*begin1 == *begin1ref);
-      ++begin1;
-      ++begin1ref;
-    }
-    CHECK(begin1 == end);
-    CHECK(begin1ref == end);
-    file1.close();
-    file1ref.close();
-  }
-}

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