Processes/ObservationPlane/CMakeLists.txt

-set (
-  MODEL_HEADERS
-  ObservationPlane.h
-  )
-
-set (
-  MODEL_SOURCES
-  ObservationPlane.cc
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/observation_plane
-  )
-
-add_library (ProcessObservationPlane STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessObservationPlane ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-#target dependencies on other libraries(also the header onlys)
-target_link_libraries (
-  ProcessObservationPlane
-  CORSIKAgeometry
-  CORSIKAprocesssequence
-  )
-
-target_include_directories (
-  ProcessObservationPlane 
-  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(testObservationPlane)
-target_link_libraries (
-  testObservationPlane
-  ProcessObservationPlane
-  CORSIKAstackinterface
-  CORSIKAthirdparty # for catch2
-)

Processes/ObservationPlane/ObservationPlane.cc

-/*
- * (c) Copyright 2019 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/observation_plane/ObservationPlane.h>
-
-#include <fstream>
-
-using namespace corsika::process::observation_plane;
-using namespace corsika::units::si;
-
-ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
-                                   std::string const& filename, bool deleteOnHit)
-    : plane_(obsPlane)
-    , outputStream_(filename)
-    , deleteOnHit_(deleteOnHit) {
-  outputStream_ << "#PDG code, energy / eV, distance to center / 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;
-  }
-
-  outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
-                << particle.GetEnergy() * (1 / 1_eV) << ' '
-                << (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
-                << std::endl;
-
-  if (deleteOnHit_) { return process::EProcessReturn::eParticleAbsorbed; }
-  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/ObservationPlane/ObservationPlane.h

-/*
- * (c) Copyright 2019 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/geometry/Plane.h>
-#include <corsika/process/ContinuousProcess.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <fstream>
-
-namespace corsika::process::observation_plane {
-
-  /**
-   * The ObservationPlane writes PDG codes, energies, and distances of particles to the
-   * central point of the plane into its output file. The particles are considered
-   * "absorbed" afterwards.
-   */
-  class ObservationPlane : public corsika::process::ContinuousProcess<ObservationPlane> {
-
-  public:
-    ObservationPlane(geometry::Plane const&, std::string const&, bool = true);
-
-    corsika::process::EProcessReturn DoContinuous(
-        corsika::setup::Stack::ParticleType const& vParticle,
-        corsika::setup::Trajectory const& vTrajectory);
-
-    corsika::units::si::LengthType MaxStepLength(
-        corsika::setup::Stack::ParticleType const&,
-        corsika::setup::Trajectory const& vTrajectory);
-
-  private:
-    geometry::Plane const plane_;
-    std::ofstream outputStream_;
-    bool const deleteOnHit_;
-  };
-} // namespace corsika::process::observation_plane

Processes/ObservationPlane/testObservationPlane.cc

-/*
- * (c) Copyright 2019 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.
- */
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-#include <corsika/process/observation_plane/ObservationPlane.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-using namespace corsika::units::si;
-using namespace corsika::process::observation_plane;
-using namespace corsika;
-using namespace corsika::geometry;
-using namespace corsika::particles;
-
-TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
-
-  auto const& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  /*
-    Test with downward going 1_GeV neutrino, starting at 0,1_m,10m
-
-    ObservationPlane has origin at 0,0,0
-   */
-
-  Point const start(rootCS, {0_m, 1_m, 10_m});
-  Vector<units::si::SpeedType::dimension_type> vec(rootCS, 0_m / second, 0_m / second,
-                                                   -units::constants::c);
-  Line line(start, vec);
-  Trajectory<Line> track(line, 12_m / units::constants::c);
-
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.Clear();
-  {
-    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-      return sqrt((Elab - m) * (Elab + m));
-    };
-    stack.AddParticle(
-        std::tuple<Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, Point,
-                   units::si::TimeType>{
-            Code::NuMu, 1_GeV,
-            corsika::stack::MomentumVector(
-                rootCS, {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::GetMass())}),
-            Point(rootCS, {1_m, 1_m, 10_m}), 0_ns});
-  }
-  auto particle = stack.GetNextParticle();
-
-  SECTION("horizontal plane") {
-
-    Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
-                         Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-    ObservationPlane obs(obsPlane, "particles.dat", true);
-
-    const LengthType length = obs.MaxStepLength(particle, track);
-    const process::EProcessReturn ret = obs.DoContinuous(particle, track);
-
-    REQUIRE(length / 10_m == Approx(1).margin(1e-4));
-    REQUIRE(ret == process::EProcessReturn::eParticleAbsorbed);
-
-    /*
-    SECTION("horizontal plane") {
-      REQUIRE(true); // todo: we have to check content of output file...
-
-    }
-    */
-  }
-
-  SECTION("inclined plane") {}
-
-  SECTION("transparent plane") {
-    Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
-                         Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-    ObservationPlane obs(obsPlane, "particles.dat", false);
-
-    const LengthType length = obs.MaxStepLength(particle, track);
-    const process::EProcessReturn ret = obs.DoContinuous(particle, track);
-
-    REQUIRE(length / 10_m == Approx(1).margin(1e-4));
-    REQUIRE(ret == process::EProcessReturn::eOk);
-  }
-}

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 = environment::Environment<setup::IEnvironmentModel>;
-  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<particles::Code, 4> particleList = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Helium,
-      particles::Code::Gamma};
-
-  std::array<double, 4> 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
-    corsika::stack::SecondaryView 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/ParticleCut/CMakeLists.txt

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

Processes/ParticleCut/ParticleCut.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/logging/Logging.h>
-#include <corsika/process/particle_cut/ParticleCut.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 particle_cut {
-
-    template <typename TParticle>
-    bool ParticleCut::ParticleIsBelowEnergyCut(TParticle const& vP) const {
-      auto const energyLab = vP.GetEnergy();
-      // nuclei
-      if (vP.GetPID() == particles::Code::Nucleus) {
-        // calculate energy per nucleon
-        auto const ElabNuc = energyLab / vP.GetNuclearA();
-        return (ElabNuc < fECut);
-      } else {
-        return (energyLab < fECut);
-      }
-    }
-
-    bool ParticleCut::ParticleIsEmParticle(Code vCode) const {
-      // FOR NOW: switch
-      switch (vCode) {
-        case Code::Gamma:
-        case Code::Electron:
-        case Code::Positron:
-          return true;
-        default:
-          return false;
-      }
-    }
-
-    bool ParticleCut::ParticleIsInvisible(Code vCode) const {
-      switch (vCode) {
-        case Code::NuE:
-        case Code::NuEBar:
-        case Code::NuMu:
-        case Code::NuMuBar:
-          return true;
-
-        default:
-          return false;
-      }
-    }
-
-    template <typename TParticle>
-    bool ParticleCut::checkCutParticle(const TParticle& particle) {
-
-      const Code pid = particle.GetPID();
-      HEPEnergyType energy = particle.GetEnergy();
-      C8LOG_DEBUG(fmt::format("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV", pid,
-                              energy / 1_GeV,
-                              (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV));
-      if (ParticleIsEmParticle(pid)) {
-        C8LOG_DEBUG("removing em. particle...");
-        fEmEnergy += energy;
-        fEmCount += 1;
-        return true;
-      } else if (ParticleIsInvisible(pid)) {
-        C8LOG_DEBUG("removing inv. particle...");
-        fInvEnergy += energy;
-        fInvCount += 1;
-        return true;
-      } else if (ParticleIsBelowEnergyCut(particle)) {
-        C8LOG_DEBUG("removing low en. particle...");
-        fEnergy += energy;
-        return true;
-      } else if (particle.GetTime() > 10_ms) {
-        C8LOG_DEBUG("removing OLD particle...");
-        fEnergy += energy;
-        return true;
-      }
-      return false; // this particle will not be removed/cut
-    }
-
-    EProcessReturn ParticleCut::DoSecondaries(corsika::setup::StackView& vS) {
-      auto particle = vS.begin();
-      while (particle != vS.end()) {
-        if (checkCutParticle(particle)) {
-          particle.Delete();
-        } else {
-          ++particle; // next entry in SecondaryView
-        }
-      }
-      return EProcessReturn::eOk;
-    }
-
-    process::EProcessReturn ParticleCut::DoContinuous(Particle& particle, Track const&) {
-      C8LOG_TRACE("ParticleCut::DoContinuous");
-      if (checkCutParticle(particle)) {
-        C8LOG_TRACE("removing during continuous");
-        return process::EProcessReturn::eParticleAbsorbed;
-      }
-      return process::EProcessReturn::eOk;
-    }
-
-    ParticleCut::ParticleCut(const units::si::HEPEnergyType vCut)
-        : fECut(vCut) {
-
-      fEmEnergy = 0._GeV;
-      fEmCount = 0;
-      fInvEnergy = 0._GeV;
-      fInvCount = 0;
-      fEnergy = 0._GeV;
-    }
-
-    void ParticleCut::ShowResults() {
-      C8LOG_INFO(fmt::format(
-          " ******************************\n"
-          " ParticleCut: \n"
-          " energy in em.  component (GeV):  {}\n"
-          " no. of em.  particles injected:  {}\n"
-          " energy in inv. component (GeV):  {}\n"
-          " no. of inv. particles injected:  {}\n"
-          " energy below particle cut (GeV): {}\n"
-          " ******************************",
-          fEmEnergy / 1_GeV, fEmCount, fInvEnergy / 1_GeV, fInvCount, fEnergy / 1_GeV));
-    }
-  } // namespace particle_cut
-} // namespace corsika::process

Processes/ParticleCut/ParticleCut.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/particles/ParticleProperties.h>
-#include <corsika/process/ContinuousProcess.h>
-#include <corsika/process/SecondariesProcess.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-  namespace particle_cut {
-    class ParticleCut : public process::SecondariesProcess<ParticleCut>,
-                        public corsika::process::ContinuousProcess<ParticleCut> {
-
-      using Particle = corsika::setup::Stack::ParticleType;
-      using Track = corsika::setup::Trajectory;
-
-      units::si::HEPEnergyType const fECut;
-
-      units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
-      units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
-      unsigned int fEmCount = 0;
-      units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
-      unsigned int fInvCount = 0;
-
-    public:
-      ParticleCut(const units::si::HEPEnergyType vCut);
-
-      EProcessReturn DoSecondaries(corsika::setup::StackView&);
-
-      EProcessReturn DoContinuous(Particle& vParticle, Track const& vTrajectory);
-
-      corsika::units::si::LengthType MaxStepLength(
-          corsika::setup::Stack::ParticleType const&, corsika::setup::Trajectory const&) {
-        return units::si::meter * std::numeric_limits<double>::infinity();
-      }
-
-      void ShowResults();
-
-      units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
-      units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
-      units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-      unsigned int GetNumberEmParticles() const { return fEmCount; }
-      unsigned int GetNumberInvParticles() const { return fInvCount; }
-
-    protected:
-      template <typename TParticle>
-      bool checkCutParticle(const TParticle& p);
-
-      template <typename TParticle>
-      bool ParticleIsBelowEnergyCut(TParticle const&) const;
-
-      bool ParticleIsEmParticle(particles::Code) const;
-      bool ParticleIsInvisible(particles::Code) const;
-    };
-  } // namespace particle_cut
-} // namespace corsika::process

Processes/ParticleCut/testParticleCut.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/particle_cut/ParticleCut.h>
-
-#include <corsika/environment/Environment.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::particle_cut;
-using namespace corsika::units;
-using namespace corsika::units::si;
-
-TEST_CASE("ParticleCut", "[processes]") {
-  feenableexcept(FE_INVALID);
-  using EnvType = environment::Environment<setup::IEnvironmentModel>;
-  EnvType env;
-  const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  // setup empty particle stack
-  setup::Stack stack;
-  stack.Clear();
-  // two energies
-  const HEPEnergyType Eabove = 1_TeV;
-  const HEPEnergyType Ebelow = 10_GeV;
-  // list of arbitrary particles
-  std::vector<particles::Code> particleList = {
-      particles::Code::PiPlus,   particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus,   particles::Code::K0Long,  particles::Code::K0Short,
-      particles::Code::Electron, particles::Code::MuPlus,  particles::Code::NuE,
-      particles::Code::Neutron};
-
-  SECTION("cut on particle type") {
-
-    ParticleCut cut(20_GeV);
-
-    // 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, Eabove,
-            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
-    corsika::stack::SecondaryView 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
-    for (auto proType : particleList)
-      projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                         corsika::stack::MomentumVector, geometry::Point,
-                                         units::si::TimeType>{
-          proType, Eabove, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
-          geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-
-    cut.DoSecondaries(view);
-
-    REQUIRE(view.GetSize() == 8);
-  }
-
-  SECTION("cut low energy") {
-    ParticleCut cut(20_GeV);
-
-    // 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, Eabove,
-            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
-    corsika::stack::SecondaryView 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 below the threshold
-    // only cut is by species
-    for (auto proType : particleList)
-      projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                         corsika::stack::MomentumVector, geometry::Point,
-                                         units::si::TimeType>{
-          proType, Ebelow, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
-          geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-
-    cut.DoSecondaries(view);
-
-    REQUIRE(view.GetSize() == 0);
-  }
-}

Processes/Pythia/CMakeLists.txt

-set (
-  MODEL_SOURCES
-  Decay.cc
-  Random.cc
-  Interaction.cc
-  )
-
-set (
-  MODEL_HEADERS
-  Decay.h
-  Random.h
-  Interaction.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/pythia
-  )
-
-add_library (ProcessPythia8 STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessPythia8 ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-
-set_target_properties (
-  ProcessPythia8
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessPythia8
-  CORSIKAprocesssequence
-  CORSIKAparticles
-  CORSIKAutilities
-  CORSIKAunits
-  CORSIKAthirdparty
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAsetup
-  CORSIKArandom
-  C8::ext::pythia8
-  )
-
-target_include_directories (
-  ProcessPythia8
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessPythia8
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  )
-
-
-# --------------------
-# code unit testing
-
-CORSIKA_ADD_TEST(testPythia8
-  SOURCES
-  testPythia8.cc
-  ${MODEL_HEADERS}
-)
-target_link_libraries (
-  testPythia8
-  ProcessPythia8
-  CORSIKAtesting
-  )
-

Processes/Pythia/Decay.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/Random.h>
-
-#include <corsika/geometry/FourVector.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/utl/COMBoost.h>
-
-using std::cout;
-using std::endl;
-using std::tuple;
-using std::vector;
-
-using namespace corsika;
-using namespace corsika::setup;
-using Projectile = corsika::setup::StackView::ParticleType;
-using Particle = corsika::setup::Stack::ParticleType;
-using Track = Trajectory;
-
-namespace corsika::process::pythia {
-
-  Decay::Decay() {
-
-    // set random number generator in pythia
-    Pythia8::RndmEngine* rndm = new corsika::process::pythia::Random();
-    fPythia.setRndmEnginePtr(rndm);
-
-    /*
-       issue xyz: definition of particles and decay channels use the same mechanism in
-       corsika and pythia we should force pythia to use the file in corsika.
-     */
-    // bool ParticleData::reInit(string startFile, bool xmlFormat = true)
-    // read in particle data from Corsika 8
-    // fPythia.particleData.reInit("/home/felix/ngcorsika/corsika-build/include/corsika/particles/ParticleData.xml");
-    // fPythia.particleData.checkTable();
-
-    fPythia.readString("Next:numberShowInfo = 0");
-    fPythia.readString("Next:numberShowProcess = 0");
-    fPythia.readString("Next:numberShowEvent = 0");
-
-    fPythia.readString("Print:quiet = on");
-    fPythia.readString("Check:particleData = 0");
-
-    /*
-       switching off event check in pythia is needed to allow decays that are off-shell
-       according to the mass definition in pythia.
-       the consistency of particle masses between event generators is an unsolved issues
-    */
-    cout << "Pythia::Init: switching off event checking in pythia.." << endl;
-    fPythia.readString("Check:event = 1");
-
-    fPythia.readString("ProcessLevel:all = off");
-    fPythia.readString("ProcessLevel:resonanceDecays = off");
-
-    // making sure
-    SetStable(particles::Code::Pi0);
-
-    //    fPythia.particleData.readString("59:m0 = 101.00");
-
-    if (!fPythia.init())
-      throw std::runtime_error("Pythia::Decay: Initialization failed!");
-  }
-
-  Decay::Decay(std::set<particles::Code> vHandled)
-      : handleAllDecays_(false)
-      , handledDecays_(vHandled) {}
-
-  Decay::~Decay() { cout << "Pythia::Decay n=" << fCount << endl; }
-
-  bool Decay::CanHandleDecay(const particles::Code vParticleCode) {
-    using namespace corsika::particles;
-    // if known to pythia and not proton, electron or neutrino it can decay
-    if (vParticleCode == Code::Proton || vParticleCode == Code::AntiProton ||
-        vParticleCode == Code::NuE || vParticleCode == Code::NuMu ||
-        vParticleCode == Code::NuTau || vParticleCode == Code::NuEBar ||
-        vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
-        vParticleCode == Code::Electron || vParticleCode == Code::Positron)
-      return false;
-    else if (CanDecay(vParticleCode)) // non-zero for particles known to sibyll
-      return true;
-    else
-      return false;
-  }
-
-  void Decay::SetHandleDecay(const particles::Code vParticleCode) {
-    handleAllDecays_ = false;
-    cout << "Pythia::Decay: set to handle decay of " << vParticleCode << endl;
-    if (Decay::CanHandleDecay(vParticleCode))
-      handledDecays_.insert(vParticleCode);
-    else
-      throw std::runtime_error("this decay can not be handled by pythia!");
-  }
-
-  void Decay::SetHandleDecay(const vector<particles::Code> vParticleList) {
-    handleAllDecays_ = false;
-    for (auto p : vParticleList) SetHandleDecay(p);
-  }
-
-  bool Decay::IsDecayHandled(const corsika::particles::Code vParticleCode) {
-    if (handleAllDecays_ && CanHandleDecay(vParticleCode))
-      return true;
-    else
-      return handledDecays_.find(vParticleCode) != Decay::handledDecays_.end();
-  }
-
-  bool Decay::IsStable(const particles::Code vCode) {
-    return fPythia.particleData.canDecay(static_cast<int>(particles::GetPDG(vCode)));
-  }
-
-  void Decay::PrintDecayConfig(const particles::Code vCode) {
-    cout << "Decay: Pythia decay configuration:" << endl;
-    cout << vCode << " is ";
-    if (IsStable(vCode))
-      cout << "stable" << endl;
-    else
-      cout << "unstable" << endl;
-  }
-
-  void Decay::PrintDecayConfig() {
-    cout << "Pythia::Decay: decay configuration:" << endl;
-    if (handleAllDecays_)
-      cout << " all particles known to Pythia are handled by Pythia::Decay!" << endl;
-    else
-      for (auto& pCode : handledDecays_)
-        cout << "Decay of " << pCode << " is handled by Pythia!" << endl;
-  }
-
-  void Decay::SetStable(const vector<particles::Code> particleList) {
-    for (auto p : particleList) Decay::SetStable(p);
-  }
-
-  bool Decay::CanDecay(const particles::Code pCode) {
-    std::cout << "Pythia::Decay: checking if particle: " << pCode
-              << " can decay in PYTHIA? ";
-    const bool ans =
-        fPythia.particleData.canDecay(static_cast<int>(particles::GetPDG(pCode)));
-    std::cout << ans << std::endl;
-    return ans;
-  }
-
-  void Decay::SetUnstable(const particles::Code pCode) {
-    cout << "Pythia::Decay: setting " << pCode << " unstable.." << endl;
-    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
-  }
-
-  void Decay::SetStable(const particles::Code pCode) {
-    cout << "Pythia::Decay: setting " << pCode << " stable.." << endl;
-    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), false);
-  }
-
-  template <>
-  units::si::TimeType Decay::GetLifetime(Particle const& vP) {
-    using namespace units::si;
-
-    const auto pid = vP.GetPID();
-    if (CanDecay(pid)) {
-      HEPEnergyType E = vP.GetEnergy();
-      HEPMassType m = vP.GetMass();
-
-      const double gamma = E / m;
-
-      const TimeType t0 = particles::GetLifetime(pid);
-      auto const lifetime = gamma * t0;
-      cout << "Pythia::Decay: code: " << vP.GetPID() << endl;
-      cout << "Pythia::Decay: MinStep: t0: " << t0 << endl;
-      cout << "Pythia::Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
-      cout << "Pythia::Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV
-           << " GeV" << endl;
-      cout << "Pythia::Decay: MinStep: gamma: " << gamma << endl;
-      cout << "Pythia::Decay: MinStep: tau: " << lifetime << endl;
-
-      return lifetime;
-    } else
-      return std::numeric_limits<double>::infinity() * 1_s;
-  }
-
-  template <>
-  void Decay::DoDecay(Projectile& vP) {
-    using geometry::Point;
-    using namespace units;
-    using namespace units::si;
-
-    auto const& decayPoint = vP.GetPosition();
-    auto const t0 = vP.GetTime();
-
-    auto const& labMomentum = vP.GetMomentum();
-    geometry::CoordinateSystem const& labCS = labMomentum.GetCoordinateSystem();
-
-    // define target kinematics in lab frame
-    // define boost to and from CoM frame
-    // CoM frame definition in Pythia projectile: +z
-    utl::COMBoost const boost(labMomentum, vP.GetMass());
-    auto const& rotatedCS = boost.GetRotatedCS();
-
-    fCount++;
-
-    // pythia stack
-    Pythia8::Event& event = fPythia.event;
-    event.reset();
-
-    auto const particleId = vP.GetPID();
-
-    // set particle unstable
-    Decay::SetUnstable(particleId);
-
-    // input particle PDG
-    auto const pdgCode = static_cast<int>(particles::GetPDG(particleId));
-
-    double constexpr px = 0;
-    double constexpr py = 0;
-    double constexpr pz = 0;
-    double const en = vP.GetMass() / 1_GeV;
-    double const m = en;
-
-    // add particle to pythia stack
-    event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
-
-    if (!fPythia.next())
-      throw std::runtime_error("Pythia::Decay: decay failed!");
-    else
-      cout << "Pythia::Decay: particles after decay: " << event.size() << endl;
-
-    // list final state
-    event.list();
-
-    // loop over final state
-    for (int i = 0; i < event.size(); ++i)
-      if (event[i].isFinal()) {
-        auto const pyId =
-            particles::ConvertFromPDG(static_cast<particles::PDGCode>(event[i].id()));
-        HEPEnergyType const Erest = event[i].e() * 1_GeV;
-        MomentumVector const pRest(
-            rotatedCS,
-            {event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
-        geometry::FourVector const fourMomRest{Erest, pRest};
-        auto const fourMomLab = boost.fromCoM(fourMomRest);
-
-        cout << "particle: id=" << pyId << " momentum="
-             << fourMomLab.GetSpaceLikeComponents().GetComponents(labCS) / 1_GeV
-             << " energy=" << fourMomLab.GetTimeLikeComponent() << endl;
-
-        vP.AddSecondary(
-            tuple<particles::Code, units::si::HEPEnergyType,
-                  corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-                pyId, fourMomLab.GetTimeLikeComponent(),
-                fourMomLab.GetSpaceLikeComponents(), decayPoint, t0});
-      }
-
-    // set particle stable
-    Decay::SetStable(particleId);
-  }
-
-} // namespace corsika::process::pythia

Processes/Pythia/Decay.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 <Pythia8/Pythia.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/DecayProcess.h>
-
-namespace corsika::process {
-
-  namespace pythia {
-
-    typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
-
-    class Decay : public corsika::process::DecayProcess<Decay> {
-      int fCount = 0;
-      bool handleAllDecays_ = true;
-
-    public:
-      Decay();
-      Decay(std::set<particles::Code>);
-      ~Decay();
-
-      // is Pythia::Decay set to handle the decay of this particle?
-      bool IsDecayHandled(const corsika::particles::Code);
-
-      // is decay possible in principle?
-      bool CanHandleDecay(const corsika::particles::Code);
-
-      // set Pythia::Decay to handle the decay of this particle!
-      void SetHandleDecay(const corsika::particles::Code);
-      // set Pythia::Decay to handle the decay of this list of particles!
-      void SetHandleDecay(const std::vector<particles::Code>);
-      // set Pythia::Decay to handle all particle decays
-      void SetHandleAllDecays();
-
-      // print internal configuration for this particle
-      void PrintDecayConfig(const corsika::particles::Code);
-      // print configuration of decays in corsika
-      void PrintDecayConfig();
-
-      bool CanDecay(const corsika::particles::Code);
-
-      /**
-       In this function PYTHIA is asked for the lifetime of the input particle.
-       Unknown particles should return an infinite lifetime so that another decay process
-       can act on the particle later on.
-     */
-
-      template <typename TParticle>
-      corsika::units::si::TimeType GetLifetime(TParticle const&);
-
-      /**
-       In this function PYTHIA is called to execute the decay of the input particle.
-     */
-
-      template <typename TProjectile>
-      void DoDecay(TProjectile&);
-
-    private:
-      void SetUnstable(const corsika::particles::Code);
-      void SetStable(const corsika::particles::Code);
-      void SetStable(const std::vector<particles::Code>);
-      bool IsStable(const corsika::particles::Code);
-
-      Pythia8::Pythia fPythia;
-      std::set<particles::Code> handledDecays_;
-    };
-
-  } // namespace pythia
-} // namespace corsika::process

Processes/Pythia/Interaction.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/pythia/Interaction.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/FourVector.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/utl/COMBoost.h>
-
-#include <tuple>
-
-using std::cout;
-using std::endl;
-using std::tuple;
-
-using namespace corsika;
-using namespace corsika::setup;
-using Projectile = corsika::setup::StackView::ParticleType;
-using Particle = corsika::setup::Stack::ParticleType;
-
-namespace corsika::process::pythia {
-
-  typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
-
-  Interaction::~Interaction() {}
-
-  Interaction::Interaction() {
-    cout << "Pythia::Interaction n=" << fCount << endl;
-
-    using random::RNGManager;
-
-    // initialize Pythia
-    if (!fInitialized) {
-
-      fPythia.readString("Print:quiet = on");
-      // TODO: proper process initialization for MinBias needed
-      fPythia.readString("HardQCD:all = on");
-      fPythia.readString("ProcessLevel:resonanceDecays = off");
-
-      fPythia.init();
-
-      // any decays in pythia? if yes need to define which particles
-      if (fInternalDecays) {
-        // define which particles are passed to corsika, i.e. which particles make it into
-        // history even very shortlived particles like charm or pi0 are of interest here
-        const std::vector<particles::Code> HadronsWeWantTrackedByCorsika = {
-            particles::Code::PiPlus,     particles::Code::PiMinus,
-            particles::Code::Pi0,        particles::Code::KMinus,
-            particles::Code::KPlus,      particles::Code::K0Long,
-            particles::Code::K0Short,    particles::Code::SigmaPlus,
-            particles::Code::SigmaMinus, particles::Code::Lambda0,
-            particles::Code::Xi0,        particles::Code::XiMinus,
-            particles::Code::OmegaMinus, particles::Code::DPlus,
-            particles::Code::DMinus,     particles::Code::D0,
-            particles::Code::D0Bar};
-
-        Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
-      }
-
-      // basic initialization of cross section routines
-      fSigma.init(&fPythia.info, fPythia.settings, &fPythia.particleData, &fPythia.rndm);
-
-      fInitialized = true;
-    }
-  }
-
-  void Interaction::SetParticleListStable(
-      std::vector<particles::Code> const& particleList) {
-    for (auto p : particleList) Interaction::SetStable(p);
-  }
-
-  void Interaction::SetUnstable(const particles::Code pCode) {
-    cout << "Pythia::Interaction: setting " << pCode << " unstable.." << endl;
-    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
-  }
-
-  void Interaction::SetStable(const particles::Code pCode) {
-    cout << "Pythia::Interaction: setting " << pCode << " stable.." << endl;
-    fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), false);
-  }
-
-  void Interaction::ConfigureLabFrameCollision(
-      const particles::Code BeamId, const particles::Code TargetId,
-      const units::si::HEPEnergyType BeamEnergy) {
-    using namespace units::si;
-    // Pythia configuration of the current event
-    // very clumsy. I am sure this can be done better..
-
-    // set beam
-    // beam id for pythia
-    auto const pdgBeam = static_cast<int>(particles::GetPDG(BeamId));
-    std::stringstream stBeam;
-    stBeam << "Beams:idA = " << pdgBeam;
-    fPythia.readString(stBeam.str());
-    // set target
-    auto pdgTarget = static_cast<int>(particles::GetPDG(TargetId));
-    // replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
-    if (TargetId == particles::Code::Hydrogen)
-      pdgTarget = static_cast<int>(particles::GetPDG(particles::Code::Proton));
-    std::stringstream stTarget;
-    stTarget << "Beams:idB = " << pdgTarget;
-    fPythia.readString(stTarget.str());
-    // set frame to lab. frame
-    fPythia.readString("Beams:frameType = 2");
-    // set beam energy
-    const double Elab = BeamEnergy / 1_GeV;
-    std::stringstream stEnergy;
-    stEnergy << "Beams:eA = " << Elab;
-    fPythia.readString(stEnergy.str());
-    // target at rest
-    fPythia.readString("Beams:eB = 0.");
-    // initialize this config
-    fPythia.init();
-  }
-
-  bool Interaction::CanInteract(const corsika::particles::Code pCode) {
-    return pCode == corsika::particles::Code::Proton ||
-           pCode == corsika::particles::Code::Neutron ||
-           pCode == corsika::particles::Code::AntiProton ||
-           pCode == corsika::particles::Code::AntiNeutron ||
-           pCode == corsika::particles::Code::PiMinus ||
-           pCode == corsika::particles::Code::PiPlus;
-  }
-
-  tuple<units::si::CrossSectionType, units::si::CrossSectionType>
-  Interaction::GetCrossSection(const particles::Code BeamId,
-                               const particles::Code TargetId,
-                               const units::si::HEPEnergyType CoMenergy) {
-    using namespace units::si;
-
-    // interaction possible in pythia?
-    if (TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen) {
-      if (CanInteract(BeamId) && ValidCoMEnergy(CoMenergy)) {
-        // input particle PDG
-        auto const pdgCodeBeam = static_cast<int>(particles::GetPDG(BeamId));
-        auto const pdgCodeTarget = static_cast<int>(particles::GetPDG(TargetId));
-        const double ecm = CoMenergy / 1_GeV;
-
-        // calculate cross section
-        fSigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
-        if (fSigma.hasSigmaTot()) {
-          const double sigEla = fSigma.sigmaEl();
-          const double sigProd = fSigma.sigmaTot() - sigEla;
-
-          return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
-
-        } else
-          throw std::runtime_error("pythia cross section init failed");
-
-      } else {
-        return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
-                               std::numeric_limits<double>::infinity() * 1_mb);
-      }
-    } else {
-      throw std::runtime_error("invalid target for pythia");
-    }
-  }
-
-  template <>
-  units::si::GrammageType Interaction::GetInteractionLength(Particle& p) {
-
-    using namespace units;
-    using namespace units::si;
-    using namespace geometry;
-
-    // coordinate system, get global frame of reference
-    CoordinateSystem& rootCS =
-        RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-    const particles::Code corsikaBeamId = p.GetPID();
-
-    // beam particles for pythia : 1, 2, 3 for p, pi, k
-    // read from cross section code table
-    const bool kInteraction = CanInteract(corsikaBeamId);
-
-    // FOR NOW: assume target is at rest
-    process::pythia::MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
-
-    // total momentum and energy
-    HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
-    process::pythia::MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
-    pTotLab += p.GetMomentum();
-    pTotLab += pTarget;
-    auto const pTotLabNorm = pTotLab.norm();
-    // calculate cm. energy
-    const HEPEnergyType ECoM = sqrt(
-        (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
-
-    cout << "Interaction: LambdaInt: \n"
-         << " input energy: " << p.GetEnergy() / 1_GeV << endl
-         << " beam can interact:" << kInteraction << endl
-         << " beam pid:" << p.GetPID() << endl;
-
-    // TODO: move limits into variables
-    if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM)) {
-
-      // get target from environment
-      /*
-        the target should be defined by the Environment,
-        ideally as full particle object so that the four momenta
-        and the boosts can be defined..
-      */
-      const auto* currentNode = p.GetNode();
-      const auto mediumComposition =
-          currentNode->GetModelProperties().GetNuclearComposition();
-      // determine average interaction length
-
-      auto const weightedProdCrossSection =
-          mediumComposition.WeightedSum([=](auto vTargetID) {
-            return std::get<0>(this->GetCrossSection(corsikaBeamId, vTargetID, ECoM));
-          });
-
-      cout << "Interaction: IntLength: weighted CrossSection (mb): "
-           << weightedProdCrossSection / 1_mb << endl
-           << "Interaction: IntLength: average mass number: "
-           << mediumComposition.GetAverageMassNumber() << endl;
-
-      // calculate interaction length in medium
-      GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
-                                      units::constants::u / weightedProdCrossSection;
-      cout << "Interaction: "
-           << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
-           << endl;
-
-      return int_length;
-    }
-
-    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-  }
-
-  /**
-     In this function PYTHIA is called to produce one event. The
-     event is copied (and boosted) into the shower lab frame.
-   */
-
-  template <>
-  process::EProcessReturn Interaction::DoInteraction(Projectile& vP) {
-
-    using namespace units;
-    using namespace utl;
-    using namespace units::si;
-    using namespace geometry;
-
-    const auto corsikaBeamId = vP.GetPID();
-    cout << "Pythia::Interaction: "
-         << "DoInteraction: " << corsikaBeamId << " interaction? "
-         << process::pythia::Interaction::CanInteract(corsikaBeamId) << endl;
-
-    if (particles::IsNucleus(corsikaBeamId)) {
-      // nuclei handled by different process, this should not happen
-      throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
-    }
-
-    if (process::pythia::Interaction::CanInteract(corsikaBeamId)) {
-
-      const CoordinateSystem& rootCS =
-          RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-      // position and time of interaction, not used in Sibyll
-      Point pOrig = vP.GetPosition();
-      TimeType tOrig = vP.GetTime();
-
-      // define target
-      // FOR NOW: target is always at rest
-      const auto eTargetLab = 0_GeV + constants::nucleonMass;
-      const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
-      const FourVector PtargLab(eTargetLab, pTargetLab);
-
-      // define projectile
-      HEPEnergyType const eProjectileLab = vP.GetEnergy();
-      auto const pProjectileLab = vP.GetMomentum();
-
-      cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
-           << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
-           << endl;
-      cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
-           << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV << endl;
-
-      const FourVector PprojLab(eProjectileLab, pProjectileLab);
-
-      // define target kinematics in lab frame
-      // define boost to and from CoM frame
-      // CoM frame definition in Pythia projectile: +z
-      COMBoost const boost(PprojLab, constants::nucleonMass);
-
-      // just for show:
-      // boost projecticle
-      auto const PprojCoM = boost.toCoM(PprojLab);
-
-      // boost target
-      auto const PtargCoM = boost.toCoM(PtargLab);
-
-      cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
-           << endl
-           << "Interaction: pbeam CoM: "
-           << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
-      cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
-           << endl
-           << "Interaction: ptarget CoM: "
-           << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
-
-      cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
-      cout << "Interaction: time: " << tOrig << endl;
-
-      HEPEnergyType Etot = eProjectileLab + eTargetLab;
-      MomentumVector Ptot = vP.GetMomentum();
-      // invariant mass, i.e. cm. energy
-      HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
-
-      // sample target mass number
-      const auto* currentNode = vP.GetNode();
-      const auto& mediumComposition =
-          currentNode->GetModelProperties().GetNuclearComposition();
-      // get cross sections for target materials
-      /*
-        Here we read the cross section from the interaction model again,
-        should be passed from GetInteractionLength if possible
-       */
-      //#warning reading interaction cross section again, should not be necessary
-      auto const& compVec = mediumComposition.GetComponents();
-      std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
-
-      for (size_t i = 0; i < compVec.size(); ++i) {
-        auto const targetId = compVec[i];
-        const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
-        [[maybe_unused]] const auto& dummy_sigEla = sigEla;
-        cross_section_of_components[i] = sigProd;
-      }
-
-      const auto corsikaTargetId =
-          mediumComposition.SampleTarget(cross_section_of_components, fRNG);
-      cout << "Interaction: target selected: " << corsikaTargetId << endl;
-
-      if (corsikaTargetId != particles::Code::Hydrogen &&
-          corsikaTargetId != particles::Code::Neutron &&
-          corsikaTargetId != particles::Code::Proton)
-        throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
-
-      cout << "Interaction: "
-           << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
-           << " Ecm(GeV): " << Ecm / 1_GeV << endl;
-
-      if (eProjectileLab < 8.5_GeV || !ValidCoMEnergy(Ecm)) {
-        cout << "Interaction: "
-             << " DoInteraction: should have dropped particle.. "
-             << "THIS IS AN ERROR" << endl;
-        throw std::runtime_error("energy too low for PYTHIA");
-
-      } else {
-        fCount++;
-
-        ConfigureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
-
-        // create event in pytia
-        if (!fPythia.next()) throw std::runtime_error("Pythia::DoInteraction: failed!");
-
-        // link to pythia stack
-        Pythia8::Event& event = fPythia.event;
-        // print final state
-        event.list();
-
-        MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
-        HEPEnergyType Elab_final = 0_GeV;
-        for (int i = 0; i < event.size(); ++i) {
-          Pythia8::Particle& p8p = event[i];
-          // skip particles that have decayed in pythia
-          if (!p8p.isFinal()) continue;
-
-          auto const pyId =
-              particles::ConvertFromPDG(static_cast<particles::PDGCode>(p8p.id()));
-
-          const MomentumVector pyPlab(
-              rootCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
-          HEPEnergyType const pyEn = p8p.e() * 1_GeV;
-
-          // add to corsika stack
-          auto pnew = vP.AddSecondary(
-              tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
-                    geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig,
-                                                          tOrig});
-
-          Plab_final += pnew.GetMomentum();
-          Elab_final += pnew.GetEnergy();
-        }
-        cout << "conservation (all GeV): "
-             << "Elab_final=" << Elab_final / 1_GeV
-             << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
-      }
-    }
-    return process::EProcessReturn::eOk;
-  }
-
-} // namespace corsika::process::pythia

Processes/Pythia/Interaction.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 <Pythia8/Pythia.h>
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/process/InteractionProcess.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <tuple>
-
-namespace corsika::process::pythia {
-
-  class Interaction : public corsika::process::InteractionProcess<Interaction> {
-
-    int fCount = 0;
-    bool fInitialized = false;
-
-  public:
-    Interaction();
-    ~Interaction();
-
-    void SetParticleListStable(std::vector<particles::Code> const&);
-    void SetUnstable(const corsika::particles::Code);
-    void SetStable(const corsika::particles::Code);
-
-    bool WasInitialized() { return fInitialized; }
-    bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
-      using namespace corsika::units::si;
-      return (10_GeV < ecm) && (ecm < 1_PeV);
-    }
-
-    bool CanInteract(const corsika::particles::Code);
-    void ConfigureLabFrameCollision(const corsika::particles::Code,
-                                    const corsika::particles::Code,
-                                    const corsika::units::si::HEPEnergyType);
-    std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
-    GetCrossSection(const corsika::particles::Code BeamId,
-                    const corsika::particles::Code TargetId,
-                    const corsika::units::si::HEPEnergyType CoMenergy);
-
-    template <typename TParticle>
-    corsika::units::si::GrammageType GetInteractionLength(TParticle&);
-
-    /**
-       In this function PYTHIA is called to produce one event. The
-       event is copied (and boosted) into the shower lab frame.
-     */
-
-    template <typename TProjectile>
-    corsika::process::EProcessReturn DoInteraction(TProjectile&);
-
-  private:
-    corsika::random::RNG& fRNG =
-        corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
-    Pythia8::Pythia fPythia;
-    Pythia8::SigmaTotal fSigma;
-    const bool fInternalDecays = true;
-  };
-
-} // namespace corsika::process::pythia

Processes/Pythia/Random.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/pythia/Random.h>
-
-namespace corsika::process::pythia {
-
-  double Random::flat() { return fDist(fRNG); }
-
-} // namespace corsika::process::pythia

Processes/Pythia/Random.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 <Pythia8/Pythia.h>
-#include <corsika/random/RNGManager.h>
-
-namespace corsika::process {
-
-  namespace pythia {
-
-    class Random : public Pythia8::RndmEngine {
-      double flat();
-
-    private:
-      std::uniform_real_distribution<double> fDist;
-      corsika::random::RNG& fRNG =
-          corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
-    };
-
-  } // namespace pythia
-} // namespace corsika::process

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/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.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") {
-
-  // setup environment, geometry
-  environment::Environment<environment::IMediumModel> env;
-
-  geometry::CoordinateSystem const& cs = env.GetCoordinateSystem();
-
-  auto theMedium =
-      environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
-          geometry::Point{cs, 0_m, 0_m, 0_m},
-          1_km * std::numeric_limits<double>::infinity());
-
-  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
-  theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_kg / (1_m * 1_m * 1_m),
-      environment::NuclearComposition(
-          std::vector<particles::Code>{particles::Code::Hydrogen},
-          std::vector<float>{1.}));
-
-  auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
-
-  SECTION("pythia decay") {
-    feenableexcept(FE_INVALID);
-    setup::Stack stack;
-    const HEPEnergyType E0 = 10_GeV;
-    HEPMomentumType P0 =
-        sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
-    auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, P0});
-    geometry::Point pos(cs, 0_m, 0_m, 0_m);
-    auto particle = stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            particles::Code::PiPlus, E0, plab, pos, 0_ns});
-
-    random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-
-    corsika::stack::SecondaryView view(particle);
-    auto projectile = view.GetProjectile();
-
-    process::pythia::Decay model;
-
-    [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
-    model.DoDecay(projectile);
-    CHECK(stack.GetSize() == 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") {
-
-    setup::Stack stack;
-    const HEPEnergyType E0 = 100_GeV;
-    HEPMomentumType P0 =
-        sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
-    auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
-    geometry::Point pos(cs, 0_m, 0_m, 0_m);
-    auto particle = stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            particles::Code::PiPlus, E0, plab, pos, 0_ns});
-    particle.SetNode(nodePtr);
-    corsika::stack::SecondaryView view(particle);
-    auto projectile = view.GetProjectile();
-
-    process::pythia::Interaction model;
-
-    [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
-    [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
-  }
-}