diff --git a/Framework/Geometry/oCoordinateSystem.h b/Framework/Geometry/oCoordinateSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..614c76da9d2f32925fa9d8c18828196fbbe8f19f
--- /dev/null
+++ b/Framework/Geometry/oCoordinateSystem.h
@@ -0,0 +1,125 @@
+/*
+ * (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
diff --git a/Framework/Utilities/otestCOMBoost.cc b/Framework/Utilities/otestCOMBoost.cc
new file mode 100644
index 0000000000000000000000000000000000000000..572839296b07634b219f518eac5192e086cda5e4
--- /dev/null
+++ b/Framework/Utilities/otestCOMBoost.cc
@@ -0,0 +1,213 @@
+/*
+ * (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()));
+ }
+}
diff --git a/Processes/Sibyll/oInteraction.cc b/Processes/Sibyll/oInteraction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..916ea2ec0829f6de41b4055428bac17b78f0fb62
--- /dev/null
+++ b/Processes/Sibyll/oInteraction.cc
@@ -0,0 +1,339 @@
+/*
+ * (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
diff --git a/Processes/Sibyll/otestSibyll.cc b/Processes/Sibyll/otestSibyll.cc
new file mode 100644
index 0000000000000000000000000000000000000000..779529e765e42aa97e10d4ec64ea7c992979f38d
--- /dev/null
+++ b/Processes/Sibyll/otestSibyll.cc
@@ -0,0 +1,216 @@
+/*
+ * (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));
+ }
+}