diff --git a/Framework/Utilities/otestCOMBoost.cc b/Framework/Utilities/otestCOMBoost.cc
deleted file mode 100644
index 572839296b07634b219f518eac5192e086cda5e4..0000000000000000000000000000000000000000
--- a/Framework/Utilities/otestCOMBoost.cc
+++ /dev/null
@@ -1,213 +0,0 @@
-/*
- * (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/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index b135f9c7a466a144111735f59b052d12a0ad05ca..c77ce9b6c0528abfe6008545611ba68a78fe0963 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -88,6 +88,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
(plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
+
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return std::numeric_limits<double>::infinity() * 1_m;
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
diff --git a/Processes/ObservationPlane/ObservationPlane_interpolation.cc b/Processes/ObservationPlane/ObservationPlane_interpolation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8460c5f5ba34947adc2e2693de897e69dbbebf7d
--- /dev/null
+++ b/Processes/ObservationPlane/ObservationPlane_interpolation.cc
@@ -0,0 +1,68 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/observation_plane/ObservationPlane.h>
+
+#include <fstream>
+
+using namespace corsika::process::observation_plane;
+using namespace corsika::units::si;
+
+ObservationPlane::ObservationPlane(
+ geometry::Plane const& obsPlane,
+ geometry::Vector<units::si::dimensionless_d> const& x_axis,
+ std::string const& filename, bool deleteOnHit)
+ : plane_(obsPlane)
+ , outputStream_(filename)
+ , deleteOnHit_(deleteOnHit)
+ , xAxis_(x_axis.normalized())
+ , yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
+ outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;
+}
+
+corsika::process::EProcessReturn ObservationPlane::DoContinuous(
+ setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
+
+ if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
+
+ if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ return process::EProcessReturn::eOk;
+ }
+
+ auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
+
+ outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
+ << particle.GetEnergy() * (1 / 1_eV) << ' '
+ << displacement.dot(xAxis_) / 1_m << ' ' << displacement.dot(yAxis_) / 1_m
+ << ' ' << std::endl;
+
+ if (deleteOnHit_) {
+ return process::EProcessReturn::eParticleAbsorbed;
+ } else {
+ return process::EProcessReturn::eOk;
+ }
+}
+
+LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+ setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+ trajectory.GetV0().dot(plane_.GetNormal());
+
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+
+ auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+ return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+}
diff --git a/Processes/Sibyll/oInteraction.cc b/Processes/Sibyll/oInteraction.cc
deleted file mode 100644
index 916ea2ec0829f6de41b4055428bac17b78f0fb62..0000000000000000000000000000000000000000
--- a/Processes/Sibyll/oInteraction.cc
+++ /dev/null
@@ -1,339 +0,0 @@
-/*
- * (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
deleted file mode 100644
index 779529e765e42aa97e10d4ec64ea7c992979f38d..0000000000000000000000000000000000000000
--- a/Processes/Sibyll/otestSibyll.cc
+++ /dev/null
@@ -1,216 +0,0 @@
-/*
- * (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));
- }
-}
diff --git a/Processes/TrackingLine/TrackingLine_interpolation.h b/Processes/TrackingLine/TrackingLine_interpolation.h
index ea0c39d21b027fe6978e8a4283e1d3ba17be1ba9..b712b6b06a8be56871143a342441065b6a51c44a 100644
--- a/Processes/TrackingLine/TrackingLine_interpolation.h
+++ b/Processes/TrackingLine/TrackingLine_interpolation.h
@@ -54,20 +54,6 @@ namespace corsika::process {
geometry::Vector<SpeedType::dimension_type> velocity =
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
- std::complex<double>* solutions = solve_quartic(1, 0, 1, -20);
- std::vector<double> tmp;
- for(int i = 0; i < 4; i++) {
- if(solutions[i].imag() == 0 && solutions[i].real() >= 0) {
- tmp.push_back(solutions[i].real());
- }
- }
- double s = *std::min_element(tmp.begin(),tmp.end());
- std::cout << "s = " << s << std::endl;
- std::cout << "x1 = " << (solutions[0].real()>=0. ? " " : "") << solutions[0].real(); if(solutions[0].imag()!=0.0) std::cout << " + i * " << solutions[0].imag(); std::cout << std::endl;
- std::cout << "x2 = " << (solutions[1].real()>=0. ? " " : "") << solutions[1].real(); if(solutions[1].imag()!=0.0) std::cout << " - i * " << -solutions[1].imag(); std::cout << std::endl;
- std::cout << "x3 = " << (solutions[2].real()>=0. ? " " : "") << solutions[2].real(); if(solutions[2].imag()!=0.0) std::cout << " + i * " << solutions[2].imag(); std::cout << std::endl;
- std::cout << "x4 = " << (solutions[3].real()>=0. ? " " : "") << solutions[3].real(); if(solutions[3].imag()!=0.0) std::cout << " - i * " << -solutions[3].imag(); std::cout << std::endl;
-
auto const currentPosition = p.GetPosition();
std::cout << "TrackingLine pid: " << p.GetPID()
<< " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
@@ -97,8 +83,11 @@ namespace corsika::process {
(corsika::units::constants::cSquared * abs(chargeNumber) *
magneticfield.GetNorm() * 1_eV);
//steplength should consider more things than just gyroradius
- double maxAngle = 0.1;
+ double maxAngle = 1e-5;
LengthType const Steplength = 2 * sin(maxAngle * M_PI / 180) * gyroradius;
+
+ std::cout << "Test: " << Steplength << " " << gyroradius << std::endl;
+
// First Movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field at position
@@ -114,11 +103,6 @@ namespace corsika::process {
velocity = direction * velocity.GetNorm();
std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
<< " GeV " << std::endl;
-
- auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
- geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
- double test =((directionBefore.cross(magneticfield)).dot(position-currentPosition) * k + 1) / (1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
- std::cout << "Test: " << test << k << std::endl;
} else {
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV