Test

Framework/Geometry/oCoordinateSystem.h

+/*
+ * (c) Copyright 2018 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.
+ */
+
+#ifndef _include_COORDINATESYSTEM_H_
+#define _include_COORDINATESYSTEM_H_
+
+#include <corsika/geometry/QuantityVector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/sgn.h>
+#include <Eigen/Dense>
+#include <stdexcept>
+
+typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
+typedef Eigen::Translation<double, 3> EigenTranslation;
+
+namespace corsika::geometry {
+
+  class RootCoordinateSystem;
+  template <typename T>
+  class Vector;
+
+  using corsika::units::si::length_d;
+
+  class CoordinateSystem {
+    CoordinateSystem const* reference = nullptr;
+    EigenTransform transf;
+
+    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
+        : reference(&reference)
+        , transf(transf) {}
+
+    CoordinateSystem()
+        : // for creating the root CS
+        transf(EigenTransform::Identity()) {}
+
+  protected:
+    static auto CreateCS() { return CoordinateSystem(); }
+    friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
+                                                    /// create ONE unique root CS
+
+  public:
+    static EigenTransform GetTransformation(CoordinateSystem const& c1,
+                                            CoordinateSystem const& c2);
+
+    auto& operator=(const CoordinateSystem& pCS) {
+      reference = pCS.reference;
+      transf = pCS.transf;
+      return *this;
+    }
+
+    auto translate(QuantityVector<length_d> vector) const {
+      EigenTransform const translation{EigenTranslation(vector.eVector)};
+
+      return CoordinateSystem(*this, translation);
+    }
+
+    template <typename TDim>
+    auto RotateToZ(Vector<TDim> vVec) const {
+      auto const a = vVec.normalized().GetComponents(*this).eVector;
+      auto const a1 = a(0), a2 = a(1);
+
+      auto const s = utl::sgn(a(2));
+      auto const c = 1 / (1 + s * a(2));
+
+      Eigen::Matrix3d A, B;
+
+      if (s > 0) {
+        A << 1, 0, a1,                      // comment to prevent clang-format
+            0, 1, a2,                       // .
+            -a1, -a2, 1;                    // .
+        B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
+            -a1 * a2 * c, -a2 * a2 * c, 0,  // .
+            0, 0, -(a1 * a1 + a2 * a2) * c; // .
+
+      } else {
+        A << 1, 0, a1,                      // .
+            0, -1, a2,                      // .
+            a1, -a2, -1;                    // .
+        B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
+            -a1 * a2 * c, +a2 * a2 * c, 0,  // .
+            0, 0, (a1 * a1 + a2 * a2) * c;  // .
+      }
+
+      return CoordinateSystem(*this, EigenTransform(A + B));
+    }
+
+    template <typename TDim>
+    auto rotate(QuantityVector<TDim> axis, double angle) const {
+      if (axis.eVector.isZero()) {
+        throw std::runtime_error("null-vector given as axis parameter");
+      }
+
+      EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
+
+      return CoordinateSystem(*this, rotation);
+    }
+
+    template <typename TDim>
+    auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
+                            QuantityVector<TDim> axis, double angle) {
+      if (axis.eVector.isZero()) {
+        throw std::runtime_error("null-vector given as axis parameter");
+      }
+
+      EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
+                                  EigenTranslation(translation.eVector)};
+
+      return CoordinateSystem(*this, transf);
+    }
+
+    auto const* GetReference() const { return reference; }
+
+    auto const& GetTransform() const { return transf; }
+  };
+
+} // namespace corsika::geometry
+
+#endif

Framework/Utilities/otestCOMBoost.cc

+/*
+ * (c) Copyright 2018 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 <catch2/catch.hpp>
+
+#include <corsika/geometry/FourVector.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+
+#include <iostream>
+
+using namespace corsika::geometry;
+using namespace corsika::utl;
+using namespace corsika::units::si;
+using corsika::units::constants::c;
+using corsika::units::constants::cSquared;
+
+double constexpr absMargin = 1e-6;
+
+CoordinateSystem const& rootCS =
+    RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+// helper function for energy-momentum
+// relativistic energy
+auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
+  return sqrt(m * m + p.squaredNorm());
+};
+
+// helper function for mandelstam-s
+auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
+  return E * E - p.squaredNorm();
+};
+
+TEST_CASE("boosts") {
+  // define target kinematics in lab frame
+  HEPMassType const targetMass = 1_GeV;
+  Vector<hepmomentum_d> pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};
+  HEPEnergyType const eTargetLab = energy(targetMass, pTargetLab);
+
+  /*
+    General tests check the interface and basic operation
+   */
+
+  SECTION("General tests") {
+
+    // define projectile kinematics in lab frame
+    HEPMassType const projectileMass = 1._GeV;
+    Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
+    HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+    const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+    // define boost to com frame
+    COMBoost boost(PprojLab, targetMass);
+
+    // boost projecticle
+    auto const PprojCoM = boost.toCoM(PprojLab);
+
+    // boost target
+    auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+    // sum of momenta in CoM, should be 0
+    auto const sumPCoM =
+        PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+    CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+
+    // mandelstam-s should be invariant under transformation
+    CHECK(s(eProjectileLab + eTargetLab,
+            pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
+              1_GeV / 1_GeV ==
+          Approx(s(PprojCoM.GetTimeLikeComponent() + PtargCoM.GetTimeLikeComponent(),
+                   PprojCoM.GetSpaceLikeComponents().GetComponents() +
+                       PtargCoM.GetSpaceLikeComponents().GetComponents()) /
+                 1_GeV / 1_GeV));
+
+    // boost back...
+    auto const PprojBack = boost.fromCoM(PprojCoM);
+
+    // ...should yield original values before the boosts
+    CHECK(PprojBack.GetTimeLikeComponent() / PprojLab.GetTimeLikeComponent() ==
+          Approx(1));
+    CHECK(
+        (PprojBack.GetSpaceLikeComponents() - PprojLab.GetSpaceLikeComponents()).norm() /
+            PprojLab.GetSpaceLikeComponents().norm() ==
+        Approx(0).margin(absMargin));
+  }
+
+  /*
+    special case: projectile along -z
+   */
+
+  SECTION("Test boost along z-axis") {
+
+    // define projectile kinematics in lab frame
+    HEPMassType const projectileMass = 1_GeV;
+    Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -1_PeV}};
+    HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+    const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+    // define boost to com frame
+    COMBoost boost(PprojLab, targetMass);
+
+    // boost projecticle
+    auto const PprojCoM = boost.toCoM(PprojLab);
+
+    // boost target
+    auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+    // sum of momenta in CoM, should be 0
+    auto const sumPCoM =
+        PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+    CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+  }
+
+  /*
+    special case: projectile with arbitrary direction
+   */
+
+  SECTION("Test boost along tilted axis") {
+
+    const HEPMomentumType P0 = 1_PeV;
+    double theta = 33.;
+    double phi = -10.;
+    auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+      return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+                             -ptot * cos(theta));
+    };
+    auto const [px, py, pz] =
+        momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+
+    // define projectile kinematics in lab frame
+    HEPMassType const projectileMass = 1_GeV;
+    Vector<hepmomentum_d> pProjectileLab(rootCS, {px, py, pz});
+    HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+    const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+    // define boost to com frame
+    COMBoost boost(PprojLab, targetMass);
+
+    // boost projecticle
+    auto const PprojCoM = boost.toCoM(PprojLab);
+
+    // boost target
+    auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+    // sum of momenta in CoM, should be 0
+    auto const sumPCoM =
+        PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+    CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+  }
+
+  /*
+    test the ultra-high energy behaviour: E=ZeV
+   */
+
+  SECTION("High energy") {
+    // define projectile kinematics in lab frame
+    HEPMassType const projectileMass = 1_GeV;
+    HEPMomentumType P0 = 1_ZeV;
+    Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -P0}};
+    HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+    const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+    // define boost to com frame
+    COMBoost boost(PprojLab, targetMass);
+
+    // boost projecticle
+    auto const PprojCoM = boost.toCoM(PprojLab);
+
+    // boost target
+    auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+
+    // sum of momenta in CoM, should be 0
+    auto const sumPCoM =
+        PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
+    CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
+  }
+
+  SECTION("rest frame") {
+    HEPMassType const projectileMass = 1_GeV;
+    HEPMomentumType const P0 = 1_TeV;
+    Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
+    HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+    const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+    COMBoost boostRest(pProjectileLab, projectileMass);
+    auto const& csPrime = boostRest.GetRotatedCS();
+    FourVector const rest4Mom = boostRest.toCoM(PprojLab);
+
+    CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
+    CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV ==
+          Approx(0).margin(absMargin));
+
+    FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
+    FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
+    auto const aLab = boostRest.fromCoM(a);
+    auto const bLab = boostRest.fromCoM(b);
+
+    CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
+    CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
+          Approx((5_GeV).magnitude()));
+    CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
+          Approx((3_GeV).magnitude()));
+  }
+}

Processes/Sibyll/oInteraction.cc

+/*
+ * (c) Copyright 2018 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/sibyll/Interaction.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/SibStack.h>
+#include <corsika/process/sibyll/sibyll2.3d.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.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 SetupParticle = setup::Stack::StackIterator;
+using SetupProjectile = setup::StackView::StackIterator;
+using Track = Trajectory;
+
+namespace corsika::process::sibyll {
+
+  Interaction::Interaction() {}
+
+  Interaction::~Interaction() {
+    cout << "Sibyll::Interaction n=" << count_ << " Nnuc=" << nucCount_ << endl;
+  }
+
+  void Interaction::Init() {
+
+    using random::RNGManager;
+
+    // initialize Sibyll
+    if (!initialized_) {
+      sibyll_ini_();
+      initialized_ = true;
+    }
+  }
+
+  void Interaction::SetAllStable() {
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
+  }
+
+  tuple<units::si::CrossSectionType, units::si::CrossSectionType>
+  Interaction::GetCrossSection(const particles::Code BeamId,
+                               const particles::Code TargetId,
+                               const units::si::HEPEnergyType CoMenergy) const {
+    using namespace units::si;
+    double sigProd, sigEla, dummy, dum1, dum3, dum4;
+    double dumdif[3];
+    const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
+    if (!IsValidCoMEnergy(CoMenergy)) {
+      throw std::runtime_error(
+          "Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
+    }
+    const double dEcm = CoMenergy / 1_GeV;
+    if (particles::IsNucleus(TargetId)) {
+      const int iTarget = particles::GetNucleusA(TargetId);
+      if (iTarget > maxTargetMassNumber_ || iTarget == 0)
+        throw std::runtime_error(
+            "Sibyll target outside range. Only nuclei with A<18 are allowed.");
+      sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
+    } else if (TargetId == particles::Proton::GetCode()) {
+      sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+    } else {
+      // no interaction in sibyll possible, return infinite cross section? or throw?
+      sigProd = std::numeric_limits<double>::infinity();
+      sigEla = std::numeric_limits<double>::infinity();
+    }
+    return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
+  }
+
+  template <>
+  units::si::GrammageType Interaction::GetInteractionLength(
+      SetupParticle const& vP) const {
+
+    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 = vP.GetPID();
+
+    // beam particles for sibyll : 1, 2, 3 for p, pi, k
+    // read from cross section code table
+    const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
+
+    // FOR NOW: assume target is at rest
+    MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
+
+    // total momentum and energy
+    HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
+    MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
+    pTotLab += vP.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: " << vP.GetEnergy() / 1_GeV << endl
+         << " beam can interact:" << kInteraction << endl
+         << " beam pid:" << vP.GetPID() << endl;
+
+    // TODO: move limits into variables
+    // FR: removed && Elab >= 8.5_GeV
+    if (kInteraction && IsValidCoMEnergy(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..
+      */
+
+      auto const* currentNode = vP.GetNode();
+      const auto& mediumComposition =
+          currentNode->GetModelProperties().GetNuclearComposition();
+
+      si::CrossSectionType weightedProdCrossSection = mediumComposition.WeightedSum(
+          [=](particles::Code targetID) -> si::CrossSectionType {
+            return std::get<0>(this->GetCrossSection(corsikaBeamId, targetID, ECoM));
+          });
+
+      cout << "Interaction: "
+           << "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
+           << 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 SIBYLL is called to produce one event. The
+     event is copied (and boosted) into the shower lab frame.
+   */
+
+  template <>
+  process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
+
+    using namespace units;
+    using namespace utl;
+    using namespace units::si;
+    using namespace geometry;
+
+    const auto corsikaBeamId = vP.GetPID();
+    cout << "ProcessSibyll: "
+         << "DoInteraction: " << corsikaBeamId << " interaction? "
+         << process::sibyll::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 SIBYLL!");
+    }
+
+    if (process::sibyll::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 Sibyll is always a single nucleon
+      // 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 Sibyll 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
+      auto const* currentNode = vP.GetNode();
+      auto const& 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 targetCode =
+          mediumComposition.SampleTarget(cross_section_of_components, RNG_);
+      cout << "Interaction: target selected: " << targetCode << endl;
+      /*
+        FOR NOW: allow nuclei with A<18 or protons only.
+        when medium composition becomes more complex, approximations will have to be
+        allowed air in atmosphere also contains some Argon.
+      */
+      int targetSibCode = -1;
+      if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
+      if (targetCode == particles::Proton::GetCode()) targetSibCode = 1;
+      cout << "Interaction: sibyll code: " << targetSibCode << endl;
+      if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
+        throw std::runtime_error(
+            "Sibyll target outside range. Only nuclei with A<18 or protons are "
+            "allowed.");
+
+      // beam id for sibyll
+      const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
+
+      cout << "Interaction: "
+           << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
+           << " Ecm(GeV): " << Ecm / 1_GeV << endl;
+      if (Ecm > GetMaxEnergyCoM())
+        throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
+      // FR: removed eProjectileLab < 8.5_GeV ||
+      if (Ecm < GetMinEnergyCoM()) {
+        cout << "Interaction: "
+             << " DoInteraction: should have dropped particle.. "
+             << "THIS IS AN ERROR" << endl;
+        throw std::runtime_error("energy too low for SIBYLL");
+      } else {
+        count_++;
+        // Sibyll does not know about units..
+        const double sqs = Ecm / 1_GeV;
+        // running sibyll, filling stack
+        sibyll_(kBeam, targetSibCode, sqs);
+
+        // print final state
+        int print_unit = 6;
+        sib_list_(print_unit);
+        nucCount_ += get_nwounded() - 1;
+
+        // add particles from sibyll to stack
+        // link to sibyll stack
+        SibStack ss;
+
+        MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+        HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
+        for (auto& psib : ss) {
+
+          // abort on particles that have decayed in Sibyll. Should not happen!
+          if (psib.HasDecayed())
+            throw std::runtime_error("found particle that decayed in SIBYLL!");
+
+          // transform energy to lab. frame
+          auto const pCoM = psib.GetMomentum();
+          HEPEnergyType const eCoM = psib.GetEnergy();
+          auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
+
+          // add to corsika stack
+          auto pnew = vP.AddSecondary(
+              tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+                    geometry::Point, units::si::TimeType>{
+                  process::sibyll::ConvertFromSibyll(psib.GetPID()),
+                  Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+                  tOrig});
+
+          Plab_final += pnew.GetMomentum();
+          Elab_final += pnew.GetEnergy();
+          Ecm_final += psib.GetEnergy();
+        }
+        cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV << endl
+             << "Elab_final=" << Elab_final / 1_GeV
+             << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
+      }
+    }
+    return process::EProcessReturn::eOk;
+  }
+
+} // namespace corsika::process::sibyll

Processes/Sibyll/otestSibyll.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/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::sibyll;
+using namespace corsika::units;
+using namespace corsika::units::si;
+
+TEST_CASE("Sibyll", "[processes]") {
+
+  SECTION("Sibyll -> Corsika") {
+    REQUIRE(particles::Electron::GetCode() ==
+            process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
+  }
+
+  SECTION("Corsika -> Sibyll") {
+    REQUIRE(process::sibyll::ConvertToSibyll(particles::Electron::GetCode()) ==
+            process::sibyll::SibyllCode::Electron);
+    REQUIRE(process::sibyll::ConvertToSibyllRaw(particles::Proton::GetCode()) == 13);
+    REQUIRE(process::sibyll::ConvertToSibyll(particles::XiStarC0::GetCode()) ==
+            process::sibyll::SibyllCode::XiStarC0);
+  }
+
+  SECTION("canInteractInSibyll") {
+
+    REQUIRE(process::sibyll::CanInteract(particles::Proton::GetCode()));
+    REQUIRE(process::sibyll::CanInteract(particles::Code::XiCPlus));
+
+    REQUIRE_FALSE(process::sibyll::CanInteract(particles::Electron::GetCode()));
+    REQUIRE_FALSE(process::sibyll::CanInteract(particles::SigmaC0::GetCode()));
+
+    REQUIRE_FALSE(process::sibyll::CanInteract(particles::Nucleus::GetCode()));
+    REQUIRE_FALSE(process::sibyll::CanInteract(particles::Helium::GetCode()));
+  }
+
+  SECTION("cross-section type") {
+
+    REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::Electron) == 0);
+    REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::K0Long) == 3);
+    REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::SigmaPlus) == 1);
+    REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::PiMinus) == 2);
+  }
+
+  SECTION("sibyll mass") {
+
+    REQUIRE_FALSE(process::sibyll::GetSibyllMass(particles::Code::Electron) == 0_GeV);
+  }
+}
+
+#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>
+#include <corsika/process/sibyll/sibyll2.3d.h>
+
+using namespace corsika::units::si;
+using namespace corsika::units;
+
+TEST_CASE("SibyllInterface", "[processes]") {
+
+  // setup environment, geometry
+  environment::Environment<environment::IMediumModel> env;
+  auto& universe = *(env.GetUniverse());
+
+  auto theMedium =
+      environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+          geometry::Point{env.GetCoordinateSystem(), 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::Oxygen}, std::vector<float>{1.}));
+
+  auto const* nodePtr = theMedium.get();
+  universe.AddChild(std::move(theMedium));
+
+  const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+
+  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+
+  SECTION("InteractionInterface") {
+
+    setup::Stack stack;
+    const HEPEnergyType E0 = 100_GeV;
+    HEPMomentumType P0 =
+        sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::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::Proton, E0, plab, pos, 0_ns});
+    particle.SetNode(nodePtr);
+    corsika::stack::SecondaryView view(particle);
+    auto projectile = view.GetProjectile();
+
+    Interaction model;
+
+    model.Init();
+    [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
+    [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+  }
+
+  SECTION("NuclearInteractionInterface") {
+
+    setup::Stack stack;
+    const HEPEnergyType E0 = 400_GeV;
+    HEPMomentumType P0 =
+        sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::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, unsigned short, unsigned short>{
+            particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
+    particle.SetNode(nodePtr);
+    corsika::stack::SecondaryView view(particle);
+    auto projectile = view.GetProjectile();
+
+    Interaction hmodel;
+    NuclearInteraction model(hmodel, env);
+
+    model.Init();
+    [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
+    [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+  }
+
+  SECTION("DecayInterface") {
+
+    setup::Stack stack;
+    const HEPEnergyType E0 = 10_GeV;
+    HEPMomentumType P0 =
+        sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::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::Lambda0, E0, plab, pos, 0_ns});
+    corsika::stack::SecondaryView view(particle);
+    auto projectile = view.GetProjectile();
+
+    Decay model;
+
+    model.PrintDecayConfig();
+
+    model.Init();
+
+    [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
+
+    /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
+
+    // run checks
+    // lambda decays into proton and pi- or neutron and pi+
+    REQUIRE(stack.GetSize() == 3);
+  }
+
+  SECTION("DecayConfiguration") {
+
+    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));
+  }
+}