diff --git a/CMakeLists.txt b/CMakeLists.txt
index d29d6566625b54239eda5b76578b97498cc97c5d..7ba7c0f9242ac697378101681271151ac56a4306 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -93,6 +93,7 @@ endif(Pythia_FOUND)
# get UrQMD
#
#find_package( UrQMD REQUIRED )
+add_subdirectory(dependencies/urqmd)
if(UrQMD_FOUND)
include_directories(${UrQMD_INCLUDE_DIR})
endif(UrQMD_FOUND)
diff --git a/corsika/detail/framework/geometry/Plane.inl b/corsika/detail/framework/geometry/Plane.inl
index d1a41c3eafebcdf6a8f62c15a6734030a39b1163..c63749215efbd087b2ec70670fd95676714209c0 100644
--- a/corsika/detail/framework/geometry/Plane.inl
+++ b/corsika/detail/framework/geometry/Plane.inl
@@ -17,22 +17,22 @@
namespace corsika {
- bool Plane::IsAbove(Point const& vP) const
+ inline bool Plane::IsAbove(Point const& vP) const
{
return fNormal.dot(vP - fCenter) > corsika::units::si::LengthType::zero();
}
- units::si::LengthType Plane::DistanceTo(corsika::Point const& vP) const
+ inline units::si::LengthType Plane::DistanceTo(corsika::Point const& vP) const
{
return (fNormal * (vP - fCenter).dot(fNormal)).norm();
}
- Point const& Plane::GetCenter() const
+ inline Point const& Plane::GetCenter() const
{
return fCenter;
}
- Plane::DimLessVec const& Plane::GetNormal() const
+ inline Plane::DimLessVec const& Plane::GetNormal() const
{
return fNormal;
}
diff --git a/corsika/detail/modules/tracking_line/TrackingLine.inl b/corsika/detail/modules/tracking_line/TrackingLine.inl
new file mode 100644
index 0000000000000000000000000000000000000000..fa9d86db8fb574491cc99433a6bb5215258355c6
--- /dev/null
+++ b/corsika/detail/modules/tracking_line/TrackingLine.inl
@@ -0,0 +1,65 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/modules/tracking_line/TrackingLine.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/media/Environment.hpp>
+
+#include <limits>
+#include <stdexcept>
+#include <utility>
+
+namespace corsika::tracking_line {
+
+ std::optional<std::pair<units::si::TimeType, units::si::TimeType>> TimeOfIntersection(
+ corsika::Line const& line, corsika::Sphere const& sphere) {
+ auto const delta = line.GetR0() - sphere.GetCenter();
+ auto const v = line.GetV0();
+ auto const vSqNorm =
+ v.squaredNorm(); // todo: get rid of this by having V0 normalized always
+ auto const R = sphere.GetRadius();
+
+ auto const vDotDelta = v.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
+
+ if (discriminant.magnitude() > 0) {
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+ return std::make_pair((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
+ } else {
+ return {};
+ }
+ }
+
+ units::si::TimeType TimeOfIntersection(Line const& vLine, Plane const& vPlane) {
+
+ using namespace units::si;
+
+ auto const delta = vPlane.GetCenter() - vLine.GetR0();
+ auto const v = vLine.GetV0();
+ auto const n = vPlane.GetNormal();
+ auto const c = n.dot(v);
+
+ if (c.magnitude() == 0) {
+ return std::numeric_limits<units::si::TimeType::value_type>::infinity() * 1_s;
+ } else {
+ return n.dot(delta) / c;
+ }
+ }
+
+} // namespace corsika::tracking_line
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
new file mode 100644
index 0000000000000000000000000000000000000000..eaacb3fb8e9b53fbd017100bdc013a4e643cf457
--- /dev/null
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -0,0 +1,327 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/modules/urqmd/UrQMD.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <algorithm>
+#include <functional>
+#include <iostream>
+#include <random>
+
+namespace corsika::urqmd {
+
+ UrQMD::UrQMD() { iniurqmdc8_(); }
+
+ using SetupStack = corsika::setup::Stack;
+ using SetupParticle = corsika::setup::Stack::StackIterator;
+ using SetupProjectile = corsika::setup::StackView::StackIterator;
+
+ corsika::units::si::CrossSectionType UrQMD::GetCrossSection(
+ corsika::Code vProjectileCode, corsika::Code vTargetCode,
+ corsika::units::si::HEPEnergyType vLabEnergy, int vAProjectile = 1) {
+
+ using namespace units::si;
+
+ // the following is a translation of ptsigtot() into C++
+ if (vProjectileCode != corsika::Code::Nucleus &&
+ !IsNucleus(vTargetCode)) { // both particles are "special"
+ auto const mProj = corsika::GetMass(vProjectileCode);
+ auto const mTar = corsika::GetMass(vTargetCode);
+ double sqrtS =
+ sqrt(units::si::detail::static_pow<2>(mProj) +
+ units::si::detail::static_pow<2>(mTar) + 2 * vLabEnergy * mTar) *
+ (1 / 1_GeV);
+
+ // we must set some UrQMD globals first...
+ auto const [ityp, iso3] = ConvertToUrQMD(vProjectileCode);
+ inputs_.spityp[0] = ityp;
+ inputs_.spiso3[0] = iso3;
+
+ auto const [itypTar, iso3Tar] = ConvertToUrQMD(vTargetCode);
+ inputs_.spityp[1] = itypTar;
+ inputs_.spiso3[1] = iso3Tar;
+
+ int one = 1;
+ int two = 2;
+ return sigtot_(one, two, sqrtS) * 1_mb;
+ } else {
+ int const Ap = vAProjectile;
+ int const At = IsNucleus(vTargetCode) ? corsika::GetNucleusA(vTargetCode) : 1;
+
+ double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1];
+ return 10_mb * M_PI * units::si::detail::static_pow<2>(maxImpact);
+ // is a constant cross-section really reasonable?
+ }
+ }
+
+ template <typename TParticle> // need template here, as this is called both with
+ // SetupParticle as well as SetupProjectile
+ corsika::units::si::CrossSectionType UrQMD::GetCrossSection(
+ TParticle const& vProjectile, corsika::Code vTargetCode) const {
+ // TODO: return 0 for non-hadrons?
+
+ auto const projectileCode = vProjectile.GetPID();
+ auto const projectileEnergyLab = vProjectile.GetEnergy();
+
+ if (projectileCode == corsika::Code::K0Long) {
+ return 0.5 *
+ (GetCrossSection(corsika::Code::K0, vTargetCode, projectileEnergyLab) +
+ GetCrossSection(corsika::Code::K0Bar, vTargetCode, projectileEnergyLab));
+ }
+
+ int const Ap =
+ (projectileCode == corsika::Code::Nucleus) ? vProjectile.GetNuclearA() : 1;
+ return GetCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
+ }
+
+ bool UrQMD::CanInteract(corsika::Code vCode) const {
+ // According to the manual, UrQMD can use all mesons, baryons and nucleons
+ // which are modeled also as input particles. I think it is safer to accept
+ // only the usual long-lived species as input.
+ // TODO: Charmed mesons should be added to the list, too
+
+ static corsika::Code const validProjectileCodes[] = {
+ corsika::Code::Nucleus, corsika::Code::Proton, corsika::Code::AntiProton,
+ corsika::Code::Neutron, corsika::Code::AntiNeutron, corsika::Code::PiPlus,
+ corsika::Code::PiMinus, corsika::Code::KPlus, corsika::Code::KMinus,
+ corsika::Code::K0, corsika::Code::K0Bar, corsika::Code::K0Long};
+
+ return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
+ vCode) != std::cend(validProjectileCodes);
+ }
+
+ corsika::units::si::GrammageType UrQMD::GetInteractionLength(
+ SetupParticle& vParticle) const {
+
+ using namespace units::si;
+
+ if (!CanInteract(vParticle.GetPID())) {
+ // we could do the canInteract check in GetCrossSection, too but if
+ // we do it here we have the advantage of avoiding the loop
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+
+ auto const& mediumComposition =
+ vParticle.GetNode()->GetModelProperties().GetNuclearComposition();
+ using namespace std::placeholders;
+
+ corsika::units::si::CrossSectionType const weightedProdCrossSection =
+ mediumComposition.WeightedSum(
+ std::bind(&UrQMD::GetCrossSection<decltype(vParticle)>, this, vParticle, _1));
+
+ return mediumComposition.GetAverageMassNumber() * units::constants::u /
+ weightedProdCrossSection;
+ }
+
+ corsika::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjectile) {
+
+ using namespace units::si;
+
+ auto projectileCode = vProjectile.GetPID();
+ auto const projectileEnergyLab = vProjectile.GetEnergy();
+ auto const& projectileMomentumLab = vProjectile.GetMomentum();
+ auto const& projectilePosition = vProjectile.GetPosition();
+ auto const projectileTime = vProjectile.GetTime();
+
+ // sample target particle
+ auto const& mediumComposition =
+ vProjectile.GetNode()->GetModelProperties().GetNuclearComposition();
+ auto const componentCrossSections = std::invoke([&]() {
+ auto const& components = mediumComposition.GetComponents();
+ std::vector<corsika::units::si::CrossSectionType> crossSections;
+ crossSections.reserve(components.size());
+
+ for (auto const c : components) {
+ crossSections.push_back(GetCrossSection(vProjectile, c));
+ }
+
+ return crossSections;
+ });
+
+ auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, fRNG);
+ auto const targetA = corsika::GetNucleusA(targetCode);
+ auto const targetZ = corsika::GetNucleusZ(targetCode);
+
+ inputs_.nevents = 1;
+ sys_.eos = 0; // could be configurable in principle
+ inputs_.outsteps = 1;
+ sys_.nsteps = 1;
+
+ // initialization regarding projectile
+ if (corsika::Code::Nucleus == projectileCode) {
+ // is this everything?
+ inputs_.prspflg = 0;
+
+ sys_.Ap = vProjectile.GetNuclearA();
+ sys_.Zp = vProjectile.GetNuclearZ();
+ rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV) /
+ vProjectile.GetNuclearA();
+
+ rsys_.bdist = nucrad_(targetA) + nucrad_(sys_.Ap) + 2 * options_.CTParam[30 - 1];
+
+ int const id = 1;
+ cascinit_(sys_.Zp, sys_.Ap, id);
+ } else {
+ inputs_.prspflg = 1;
+ sys_.Ap = 1; // even for non-baryons this has to be set, see vanilla UrQMD.f
+ rsys_.bdist = nucrad_(targetA) + nucrad_(1) + 2 * options_.CTParam[30 - 1];
+ rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV);
+
+ if (projectileCode == corsika::Code::K0Long) {
+ projectileCode = fBooleanDist(fRNG) ? corsika::Code::K0 : corsika::Code::K0Bar;
+ } else if (projectileCode == corsika::Code::K0Short) {
+ throw std::runtime_error("K0Short should not interact");
+ }
+
+ auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
+ // todo: conversion of K_long/short into strong eigenstates;
+ inputs_.spityp[0] = ityp;
+ inputs_.spiso3[0] = iso3;
+ }
+
+ // initilazation regarding target
+ if (corsika::IsNucleus(targetCode)) {
+ sys_.Zt = targetZ;
+ sys_.At = targetA;
+ inputs_.trspflg = 0; // nucleus as target
+ int const id = 2;
+ cascinit_(sys_.Zt, sys_.At, id);
+ } else {
+ inputs_.trspflg = 1; // special particle as target
+ auto const [ityp, iso3] = ConvertToUrQMD(targetCode);
+ inputs_.spityp[1] = ityp;
+ inputs_.spiso3[1] = iso3;
+ }
+
+ int iflb = 0; // flag for retrying interaction in case of empty event, 0 means retry
+ urqmd_(iflb);
+
+ // now retrieve secondaries from UrQMD
+ auto const& originalCS = projectileMomentumLab.GetCoordinateSystem();
+ corsika::CoordinateSystem const zAxisFrame =
+ originalCS.RotateToZ(projectileMomentumLab);
+
+ for (int i = 0; i < sys_.npart; ++i) {
+ auto code = ConvertFromUrQMD(isys_.ityp[i], isys_.iso3[i]);
+ if (code == corsika::Code::K0 || code == corsika::Code::K0Bar) {
+ code = fBooleanDist(fRNG) ? corsika::Code::K0Short : corsika::Code::K0Long;
+ }
+
+ // "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
+ auto momentum = corsika::Vector(
+ zAxisFrame, corsika::QuantityVector<dimensionless_d>{coor_.px[i], coor_.py[i],
+ coor_.pz[i]} *
+ 1_GeV);
+
+ auto const energy = sqrt(momentum.squaredNorm() + square(corsika::GetMass(code)));
+
+ momentum.rebase(originalCS); // transform back into standard lab frame
+ std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl;
+
+ vProjectile.AddSecondary(
+ std::tuple<corsika::Code, corsika::units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point, corsika::units::si::TimeType>{
+ code, energy, momentum, projectilePosition, projectileTime});
+ }
+
+ std::cout << "UrQMD generated " << sys_.npart << " secondaries!" << std::endl;
+
+ return corsika::EProcessReturn::eOk;
+ }
+
+ /**
+ * the random number generator function of UrQMD
+ */
+ double ranf_(int&) {
+ static corsika::RNG& rng =
+ corsika::RNGManager::GetInstance().GetRandomStream("UrQMD");
+ static std::uniform_real_distribution<double> dist;
+
+ return dist(rng);
+ }
+
+ corsika::Code ConvertFromUrQMD(int vItyp, int vIso3) {
+ int const pdgInt =
+ pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
+ if (pdgInt == 0) { // pdgid_ returns 0 on error
+ throw std::runtime_error("UrQMD pdgid() returned 0");
+ }
+ auto const pdg = static_cast<corsika::PDGCode>(pdgInt);
+ return corsika::ConvertFromPDG(pdg);
+ }
+
+ std::pair<int, int> ConvertToUrQMD(corsika::Code code) {
+ static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
+ // data mostly from github.com/afedynitch/ParticleDataTool
+ {22, {100, 0}}, // gamma
+ {111, {101, 0}}, // pi0
+ {211, {101, 2}}, // pi+
+ {-211, {101, -2}}, // pi-
+ {321, {106, 1}}, // K+
+ {-321, {-106, -1}}, // K-
+ {311, {106, -1}}, // K0
+ {-311, {-106, 1}}, // K0bar
+ {2212, {1, 1}}, // p
+ {2112, {1, -1}}, // n
+ {-2212, {-1, -1}}, // pbar
+ {-2112, {-1, 1}}, // nbar
+ {221, {102, 0}}, // eta
+ {213, {104, 2}}, // rho+
+ {-213, {104, -2}}, // rho-
+ {113, {104, 0}}, // rho0
+ {323, {108, 2}}, // K*+
+ {-323, {108, -2}}, // K*-
+ {313, {108, 0}}, // K*0
+ {-313, {-108, 0}}, // K*0-bar
+ {223, {103, 0}}, // omega
+ {333, {109, 0}}, // phi
+ {3222, {40, 2}}, // Sigma+
+ {3212, {40, 0}}, // Sigma0
+ {3112, {40, -2}}, // Sigma-
+ {3322, {49, 0}}, // Xi0
+ {3312, {49, -1}}, // Xi-
+ {3122, {27, 0}}, // Lambda0
+ {2224, {17, 4}}, // Delta++
+ {2214, {17, 2}}, // Delta+
+ {2114, {17, 0}}, // Delta0
+ {1114, {17, -2}}, // Delta-
+ {3224, {41, 2}}, // Sigma*+
+ {3214, {41, 0}}, // Sigma*0
+ {3114, {41, -2}}, // Sigma*-
+ {3324, {50, 0}}, // Xi*0
+ {3314, {50, -1}}, // Xi*-
+ {3334, {55, 0}}, // Omega-
+ {411, {133, 2}}, // D+
+ {-411, {133, -2}}, // D-
+ {421, {133, 0}}, // D0
+ {-421, {-133, 0}}, // D0-bar
+ {441, {107, 0}}, // etaC
+ {431, {138, 1}}, // Ds+
+ {-431, {138, -1}}, // Ds-
+ {433, {139, 1}}, // Ds*+
+ {-433, {139, -1}}, // Ds*-
+ {413, {134, 1}}, // D*+
+ {-413, {134, -1}}, // D*-
+ {10421, {134, 0}}, // D*0
+ {-10421, {-134, 0}}, // D*0-bar
+ {443, {135, 0}}, // jpsi
+ };
+
+ return mapPDGToUrQMD.at(static_cast<int>(GetPDG(code)));
+ }
+
+} // namespace corsika::urqmd
\ No newline at end of file
diff --git a/corsika/modules/switch_process/SwitchProcess.hpp b/corsika/modules/switch_process/SwitchProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3d083535342f86cf13ac1701e4f55a994d9b7767
--- /dev/null
+++ b/corsika/modules/switch_process/SwitchProcess.hpp
@@ -0,0 +1,106 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/sequence/InteractionProcess.hpp>
+#include <corsika/framework/sequence/ProcessSequence.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika::switch_process {
+
+ /**
+ * This process provides an energy-based switch between two interaction processes P1 and
+ * P1. For energies below the threshold, P1 is invoked, otherwise P2. Both can be either
+ * single interaction processes or multiple ones combined in a ProcessSequence. A
+ * SwitchProcess itself will always be regarded as a distinct case when assembled into a
+ * (greater) ProcessSequence.
+ */
+
+ template <class TLowEProcess, class THighEProcess>
+ class SwitchProcess : public BaseProcess<SwitchProcess<TLowEProcess, THighEProcess>> {
+ TLowEProcess& fLowEProcess;
+ THighEProcess& fHighEProcess;
+ units::si::HEPEnergyType const fThresholdEnergy;
+
+ public:
+ SwitchProcess(TLowEProcess& vLowEProcess, THighEProcess& vHighEProcess,
+ units::si::HEPEnergyType vThresholdEnergy)
+ : fLowEProcess(vLowEProcess)
+ , fHighEProcess(vHighEProcess)
+ , fThresholdEnergy(vThresholdEnergy) {}
+
+ void Init() {
+ fLowEProcess.Init();
+ fHighEProcess.Init();
+ }
+
+ template <typename TParticle>
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
+ return 1 / GetInteractionLength(p);
+ }
+
+ template <typename TParticle>
+ units::si::GrammageType GetInteractionLength(TParticle& vParticle) {
+ if (vParticle.GetEnergy() < fThresholdEnergy) {
+ if constexpr (is_process_sequence_v<TLowEProcess>) {
+ return fLowEProcess.GetTotalInteractionLength(vParticle);
+ } else {
+ return fLowEProcess.GetInteractionLength(vParticle);
+ }
+ } else {
+ if constexpr (is_process_sequence_v<THighEProcess>) {
+ return fHighEProcess.GetTotalInteractionLength(vParticle);
+ } else {
+ return fHighEProcess.GetInteractionLength(vParticle);
+ }
+ }
+ }
+
+ // required to conform to ProcessSequence interface. We cannot just
+ // implement DoInteraction() because we want to call SelectInteraction
+ // in case a member process is a ProcessSequence.
+ template <typename TParticle, typename TSecondaries>
+ EProcessReturn SelectInteraction(
+ TParticle& vP, TSecondaries& vS,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
+ corsika::units::si::InverseGrammageType& lambda_inv_count) {
+ if (vP.GetEnergy() < fThresholdEnergy) {
+ if constexpr (is_process_sequence_v<TLowEProcess>) {
+ return fLowEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
+ } else {
+ lambda_inv_count += fLowEProcess.GetInverseInteractionLength(vP);
+ // check if we should execute THIS process and then EXIT
+ if (lambda_select < lambda_inv_count) {
+ fLowEProcess.DoInteraction(vS);
+ return EProcessReturn::eInteracted;
+ } else {
+ return EProcessReturn::eOk;
+ }
+ }
+ } else {
+ if constexpr (is_process_sequence_v<THighEProcess>) {
+ return fHighEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
+ } else {
+ lambda_inv_count += fHighEProcess.GetInverseInteractionLength(vP);
+ // check if we should execute THIS process and then EXIT
+ if (lambda_select < lambda_inv_count) {
+ fHighEProcess.DoInteraction(vS);
+ return EProcessReturn::eInteracted;
+ } else {
+ return EProcessReturn::eOk;
+ }
+ }
+ }
+ }
+ };
+} // namespace corsika::switch_process
+
diff --git a/corsika/modules/tracking_line/TrackingLine.hpp b/corsika/modules/tracking_line/TrackingLine.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f93f931f198eb97e3acf0d01f98ebd2d70dfd2ea
--- /dev/null
+++ b/corsika/modules/tracking_line/TrackingLine.hpp
@@ -0,0 +1,138 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <optional>
+#include <type_traits>
+#include <utility>
+
+namespace corsika {
+ template <typename IEnvironmentModel>
+ class Environment;
+ template <typename IEnvironmentModel>
+ class VolumeTreeNode;
+} // namespace corsika
+
+namespace corsika::tracking_line {
+
+ std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+ TimeOfIntersection(corsika::Line const&, corsika::Sphere const&);
+
+ corsika::units::si::TimeType TimeOfIntersection(corsika::Line const&,
+ corsika::Plane const&);
+
+ class TrackingLine {
+
+ public:
+ TrackingLine(){};
+
+ template <typename Particle> // was Stack previously, and argument was
+ // Stack::StackIterator
+ auto GetTrack(Particle const& p) {
+ using namespace corsika::units::si;
+ corsika::Vector<SpeedType::dimension_type> const velocity =
+ p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
+
+ auto const currentPosition = p.GetPosition();
+ std::cout << "TrackingLine pid: " << p.GetPID()
+ << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine pos: "
+ << currentPosition.GetCoordinates()
+ // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
+ << std::endl;
+ std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
+ << " GeV " << std::endl;
+ std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+
+ // to do: include effect of magnetic field
+ corsika::Line line(currentPosition, velocity);
+
+ auto const* currentLogicalVolumeNode = p.GetNode();
+ //~ auto const* currentNumericalVolumeNode =
+ //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+ auto const numericallyInside =
+ currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
+
+ std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
+
+ auto const& children = currentLogicalVolumeNode->GetChildNodes();
+ auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
+
+ std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
+
+ // for entering from outside
+ auto addIfIntersects = [&](auto const& vtn) {
+ auto const& volume = vtn.GetVolume();
+ auto const& sphere = dynamic_cast<corsika::Sphere const&>(
+ volume); // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+
+ if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+ auto const [t1, t2] = *opt;
+ std::cout << "intersection times: " << t1 / 1_s << "; "
+ << t2 / 1_s
+ // << " " << vtn.GetModelProperties().GetName()
+ << std::endl;
+ if (t1.magnitude() > 0)
+ intersections.emplace_back(t1, &vtn);
+ else if (t2.magnitude() > 0)
+ std::cout << "inside other volume" << std::endl;
+ }
+ };
+
+ for (auto const& child : children) { addIfIntersects(*child); }
+ for (auto const* ex : excluded) { addIfIntersects(*ex); }
+
+ {
+ auto const& sphere =
+ dynamic_cast<corsika::Sphere const&>(currentLogicalVolumeNode->GetVolume());
+ // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+ [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+ [[maybe_unused]] auto dummy_t1 = t1;
+ intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
+ }
+
+ auto const minIter = std::min_element(
+ intersections.cbegin(), intersections.cend(),
+ [](auto const& a, auto const& b) { return a.first < b.first; });
+
+ TimeType min;
+
+ if (minIter == intersections.cend()) {
+ min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+ // to handle the numerics properly
+ throw std::runtime_error("no intersection with anything!");
+ } else {
+ min = minIter->first;
+ }
+
+ std::cout << " t-intersect: "
+ << min
+ // << " " << minIter->second->GetModelProperties().GetName()
+ << std::endl;
+
+ return std::make_tuple(corsika::Trajectory<corsika::Line>(line, min),
+ velocity.norm() * min, minIter->second);
+ }
+ };
+
+} // namespace corsika::tracking_line
+
+#include <corsika/detail/modules/tracking_line/TrackingLine.inl>
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..71abe44adfd72eafd26b28fbf60d683c9fe312fd
--- /dev/null
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -0,0 +1,141 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/sequence/InteractionProcess.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <array>
+#include <utility>
+
+namespace corsika::urqmd {
+
+ class UrQMD : public corsika::InteractionProcess<UrQMD> {
+ public:
+ UrQMD();
+ void Init() {}
+ corsika::units::si::GrammageType GetInteractionLength(
+ corsika::setup::Stack::StackIterator&) const;
+
+ template <typename TParticle>
+ corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,
+ corsika::Code) const;
+
+ corsika::EProcessReturn DoInteraction(
+ corsika::setup::StackView::StackIterator&);
+
+ bool CanInteract(corsika::Code) const;
+
+ private:
+ static corsika::units::si::CrossSectionType GetCrossSection(
+ corsika::Code, corsika::Code, corsika::units::si::HEPEnergyType, int);
+ corsika::RNG& fRNG =
+ corsika::RNGManager::GetInstance().GetRandomStream("UrQMD");
+
+ std::uniform_int_distribution<int> fBooleanDist{0, 1};
+ };
+
+ namespace constants {
+ // from coms.f
+ int constexpr nmax = 500;
+
+ // from options.f
+ int constexpr numcto = 400;
+ int constexpr numctp = 400;
+
+ // from inputs.f
+ int constexpr aamax = 300;
+
+ } // namespace constants
+
+ template <typename T>
+ using nmaxArray = std::array<T, constants::nmax>;
+ using nmaxIntArray = nmaxArray<int>;
+ using nmaxDoubleArray = nmaxArray<double>;
+
+ extern "C" {
+ void iniurqmdc8_();
+ double ranf_(int&);
+ void cascinit_(int const&, int const&, int const&);
+ double nucrad_(int const&);
+ void urqmd_(int&);
+ int pdgid_(int&, int&);
+ double sigtot_(int&, int&, double&);
+
+ // defined in coms.f
+ extern struct {
+ int npart, nbar, nmes, ctag, nsteps, uid_cnt, ranseed, event;
+ int Ap; // projectile mass number (in case of nucleus)
+ int At; // target mass number (in case of nucleus)
+ int Zp; // projectile charge number (in case of nucleus)
+ int Zt; // target charge number (in case of nucleus)
+ int eos, dectag, NHardRes, NSoftRes, NDecRes, NElColl, NBlColl;
+ } sys_;
+
+ extern struct {
+ double time, acttime, bdist, bimp, bmin;
+ double ebeam; // lab-frame energy of projectile
+ double ecm;
+ } rsys_;
+
+ // defined in coms.f
+ extern struct {
+ nmaxIntArray spin, ncoll, charge, ityp, lstcoll, iso3, origin, strid, uid;
+ } isys_;
+
+ // defined in coor.f
+ extern struct {
+ nmaxDoubleArray r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime;
+ } coor_;
+
+ // defined in inputs.f
+ extern struct {
+ int nevents;
+ std::array<int, 2> spityp; // particle codes of: [0]: projectile, [1]: target
+ int prspflg; // projectile special flag
+ int trspflg; // target special flag, set to 1 unless target is nucleus > H
+ std::array<int, 2> spiso3; // particle codes of: [0]: projectile, [1]: target
+ int outsteps, bflag, srtflag, efuncflag, nsrt, npb, firstev;
+ } inputs_;
+
+ // defined in inputs.f
+ extern struct {
+ double srtmin, srtmax, pbeam, betann, betatar, betapro, pbmin, pbmax;
+ } input2_;
+
+ // defined in options.f
+ extern struct {
+ std::array<double, constants::numcto> CTOption;
+ std::array<double, constants::numctp> CTParam;
+ } options_;
+
+ extern struct {
+ int fixedseed, bf13, bf14, bf15, bf16, bf17, bf18, bf19, bf20;
+ } loptions_;
+
+ // defined in urqmdInterface.F
+ extern struct { std::array<double, 3> xs, bim; } cxs_u2_;
+ }
+
+ /**
+ * convert CORSIKA code to UrQMD code tuple
+ *
+ * In the current implementation a detour via the PDG code is made.
+ */
+ std::pair<int, int> ConvertToUrQMD(corsika::Code);
+ corsika::Code ConvertFromUrQMD(int vItyp, int vIso3);
+
+} // namespace corsika::urqmd
+
+#include <corsika/detail/modules/urqmd/UrQMD.inl>
diff --git a/dependencies/UrQMD/CMakeLists.txt b/dependencies/urqmd/CMakeLists.txt
similarity index 78%
rename from dependencies/UrQMD/CMakeLists.txt
rename to dependencies/urqmd/CMakeLists.txt
index 6c538f7ff17f4e4d0fcb1a37f2b71c6ab6cf04c7..45b7c244a4ede4722c9a1df0dc9feac755eec9bf 100644
--- a/dependencies/UrQMD/CMakeLists.txt
+++ b/dependencies/urqmd/CMakeLists.txt
@@ -1,6 +1,5 @@
-cmake_minimum_required(VERSION 3.1)
-
-project(libUrQMD)
+#cmake_minimum_required(VERSION 3.1)
+#project(libUrQMD)
set (
MODEL_SOURCES
@@ -49,3 +48,6 @@ add_library (UrQMD_static STATIC ${MODEL_SOURCES})
set_property(TARGET UrQMD_static PROPERTY POSITION_INDEPENDENT_CODE 1)
+# add UrQMD to CORSIKA8 build
+add_dependencies (CORSIKA8 UrQMD_static)
+target_link_libraries (CORSIKA8 INTERFACE UrQMD_static)
diff --git a/dependencies/UrQMD/Copyright b/dependencies/urqmd/Copyright
similarity index 100%
rename from dependencies/UrQMD/Copyright
rename to dependencies/urqmd/Copyright
diff --git a/dependencies/UrQMD/README b/dependencies/urqmd/README
similarity index 100%
rename from dependencies/UrQMD/README
rename to dependencies/urqmd/README
diff --git a/dependencies/UrQMD/UrQMD.cc b/dependencies/urqmd/UrQMD.cc
similarity index 100%
rename from dependencies/UrQMD/UrQMD.cc
rename to dependencies/urqmd/UrQMD.cc
diff --git a/dependencies/UrQMD/UrQMD.h b/dependencies/urqmd/UrQMD.h
similarity index 100%
rename from dependencies/UrQMD/UrQMD.h
rename to dependencies/urqmd/UrQMD.h
diff --git a/dependencies/UrQMD/addpart.f b/dependencies/urqmd/addpart.f
similarity index 100%
rename from dependencies/UrQMD/addpart.f
rename to dependencies/urqmd/addpart.f
diff --git a/dependencies/UrQMD/angdis.f b/dependencies/urqmd/angdis.f
similarity index 100%
rename from dependencies/UrQMD/angdis.f
rename to dependencies/urqmd/angdis.f
diff --git a/dependencies/UrQMD/anndec.f b/dependencies/urqmd/anndec.f
similarity index 100%
rename from dependencies/UrQMD/anndec.f
rename to dependencies/urqmd/anndec.f
diff --git a/dependencies/UrQMD/blockres.f b/dependencies/urqmd/blockres.f
similarity index 100%
rename from dependencies/UrQMD/blockres.f
rename to dependencies/urqmd/blockres.f
diff --git a/dependencies/UrQMD/boxinc.f b/dependencies/urqmd/boxinc.f
similarity index 100%
rename from dependencies/UrQMD/boxinc.f
rename to dependencies/urqmd/boxinc.f
diff --git a/dependencies/UrQMD/boxprg.f b/dependencies/urqmd/boxprg.f
similarity index 100%
rename from dependencies/UrQMD/boxprg.f
rename to dependencies/urqmd/boxprg.f
diff --git a/dependencies/UrQMD/cascinit.f b/dependencies/urqmd/cascinit.f
similarity index 100%
rename from dependencies/UrQMD/cascinit.f
rename to dependencies/urqmd/cascinit.f
diff --git a/dependencies/UrQMD/colltab.f b/dependencies/urqmd/colltab.f
similarity index 100%
rename from dependencies/UrQMD/colltab.f
rename to dependencies/urqmd/colltab.f
diff --git a/dependencies/UrQMD/coload.f b/dependencies/urqmd/coload.f
similarity index 100%
rename from dependencies/UrQMD/coload.f
rename to dependencies/urqmd/coload.f
diff --git a/dependencies/UrQMD/comnorm.f b/dependencies/urqmd/comnorm.f
similarity index 100%
rename from dependencies/UrQMD/comnorm.f
rename to dependencies/urqmd/comnorm.f
diff --git a/dependencies/UrQMD/comres.f b/dependencies/urqmd/comres.f
similarity index 100%
rename from dependencies/UrQMD/comres.f
rename to dependencies/urqmd/comres.f
diff --git a/dependencies/UrQMD/coms.f b/dependencies/urqmd/coms.f
similarity index 100%
rename from dependencies/UrQMD/coms.f
rename to dependencies/urqmd/coms.f
diff --git a/dependencies/UrQMD/comstr.f b/dependencies/urqmd/comstr.f
similarity index 100%
rename from dependencies/UrQMD/comstr.f
rename to dependencies/urqmd/comstr.f
diff --git a/dependencies/UrQMD/comwid.f b/dependencies/urqmd/comwid.f
similarity index 100%
rename from dependencies/UrQMD/comwid.f
rename to dependencies/urqmd/comwid.f
diff --git a/dependencies/UrQMD/dectim.f b/dependencies/urqmd/dectim.f
similarity index 100%
rename from dependencies/UrQMD/dectim.f
rename to dependencies/urqmd/dectim.f
diff --git a/dependencies/UrQMD/delpart.f b/dependencies/urqmd/delpart.f
similarity index 100%
rename from dependencies/UrQMD/delpart.f
rename to dependencies/urqmd/delpart.f
diff --git a/dependencies/UrQMD/detbal.f b/dependencies/urqmd/detbal.f
similarity index 100%
rename from dependencies/UrQMD/detbal.f
rename to dependencies/urqmd/detbal.f
diff --git a/dependencies/UrQMD/dwidth.f b/dependencies/urqmd/dwidth.f
similarity index 100%
rename from dependencies/UrQMD/dwidth.f
rename to dependencies/urqmd/dwidth.f
diff --git a/dependencies/UrQMD/error.f b/dependencies/urqmd/error.f
similarity index 100%
rename from dependencies/UrQMD/error.f
rename to dependencies/urqmd/error.f
diff --git a/dependencies/UrQMD/freezeout.f b/dependencies/urqmd/freezeout.f
similarity index 100%
rename from dependencies/UrQMD/freezeout.f
rename to dependencies/urqmd/freezeout.f
diff --git a/dependencies/UrQMD/getmass.f b/dependencies/urqmd/getmass.f
similarity index 100%
rename from dependencies/UrQMD/getmass.f
rename to dependencies/urqmd/getmass.f
diff --git a/dependencies/UrQMD/getspin.f b/dependencies/urqmd/getspin.f
similarity index 100%
rename from dependencies/UrQMD/getspin.f
rename to dependencies/urqmd/getspin.f
diff --git a/dependencies/UrQMD/init.f b/dependencies/urqmd/init.f
similarity index 100%
rename from dependencies/UrQMD/init.f
rename to dependencies/urqmd/init.f
diff --git a/dependencies/UrQMD/inputs.f b/dependencies/urqmd/inputs.f
similarity index 100%
rename from dependencies/UrQMD/inputs.f
rename to dependencies/urqmd/inputs.f
diff --git a/dependencies/UrQMD/iso.f b/dependencies/urqmd/iso.f
similarity index 100%
rename from dependencies/UrQMD/iso.f
rename to dependencies/urqmd/iso.f
diff --git a/dependencies/UrQMD/ityp2pdg.f b/dependencies/urqmd/ityp2pdg.f
similarity index 100%
rename from dependencies/UrQMD/ityp2pdg.f
rename to dependencies/urqmd/ityp2pdg.f
diff --git a/dependencies/UrQMD/jdecay2.f b/dependencies/urqmd/jdecay2.f
similarity index 100%
rename from dependencies/UrQMD/jdecay2.f
rename to dependencies/urqmd/jdecay2.f
diff --git a/dependencies/UrQMD/make22.f b/dependencies/urqmd/make22.f
similarity index 100%
rename from dependencies/UrQMD/make22.f
rename to dependencies/urqmd/make22.f
diff --git a/dependencies/UrQMD/newpart.f b/dependencies/urqmd/newpart.f
similarity index 100%
rename from dependencies/UrQMD/newpart.f
rename to dependencies/urqmd/newpart.f
diff --git a/dependencies/UrQMD/numrec.f b/dependencies/urqmd/numrec.f
similarity index 100%
rename from dependencies/UrQMD/numrec.f
rename to dependencies/urqmd/numrec.f
diff --git a/dependencies/UrQMD/options.f b/dependencies/urqmd/options.f
similarity index 100%
rename from dependencies/UrQMD/options.f
rename to dependencies/urqmd/options.f
diff --git a/dependencies/UrQMD/outcom.f b/dependencies/urqmd/outcom.f
similarity index 100%
rename from dependencies/UrQMD/outcom.f
rename to dependencies/urqmd/outcom.f
diff --git a/dependencies/UrQMD/output.f b/dependencies/urqmd/output.f
similarity index 100%
rename from dependencies/UrQMD/output.f
rename to dependencies/urqmd/output.f
diff --git a/dependencies/UrQMD/paulibl.f b/dependencies/urqmd/paulibl.f
similarity index 100%
rename from dependencies/UrQMD/paulibl.f
rename to dependencies/urqmd/paulibl.f
diff --git a/dependencies/UrQMD/proppot.f b/dependencies/urqmd/proppot.f
similarity index 100%
rename from dependencies/UrQMD/proppot.f
rename to dependencies/urqmd/proppot.f
diff --git a/dependencies/UrQMD/saveinfo.f b/dependencies/urqmd/saveinfo.f
similarity index 100%
rename from dependencies/UrQMD/saveinfo.f
rename to dependencies/urqmd/saveinfo.f
diff --git a/dependencies/UrQMD/scatter.f b/dependencies/urqmd/scatter.f
similarity index 100%
rename from dependencies/UrQMD/scatter.f
rename to dependencies/urqmd/scatter.f
diff --git a/dependencies/UrQMD/siglookup.f b/dependencies/urqmd/siglookup.f
similarity index 100%
rename from dependencies/UrQMD/siglookup.f
rename to dependencies/urqmd/siglookup.f
diff --git a/dependencies/UrQMD/string.f b/dependencies/urqmd/string.f
similarity index 100%
rename from dependencies/UrQMD/string.f
rename to dependencies/urqmd/string.f
diff --git a/dependencies/UrQMD/tabinit.f b/dependencies/urqmd/tabinit.f
similarity index 100%
rename from dependencies/UrQMD/tabinit.f
rename to dependencies/urqmd/tabinit.f
diff --git a/dependencies/UrQMD/testUrQMD.cc b/dependencies/urqmd/testUrQMD.cc
similarity index 100%
rename from dependencies/UrQMD/testUrQMD.cc
rename to dependencies/urqmd/testUrQMD.cc
diff --git a/dependencies/UrQMD/urqmd.f b/dependencies/urqmd/urqmd.f
similarity index 100%
rename from dependencies/UrQMD/urqmd.f
rename to dependencies/urqmd/urqmd.f
diff --git a/dependencies/UrQMD/urqmdInterface.F b/dependencies/urqmd/urqmdInterface.F
similarity index 100%
rename from dependencies/UrQMD/urqmdInterface.F
rename to dependencies/urqmd/urqmdInterface.F
diff --git a/dependencies/UrQMD/urqmd_xs.cc b/dependencies/urqmd/urqmd_xs.cc
similarity index 100%
rename from dependencies/UrQMD/urqmd_xs.cc
rename to dependencies/urqmd/urqmd_xs.cc
diff --git a/dependencies/UrQMD/whichres.f b/dependencies/urqmd/whichres.f
similarity index 100%
rename from dependencies/UrQMD/whichres.f
rename to dependencies/urqmd/whichres.f
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 41a0b6ca65c3c6fb11def46aa76b7448b263aa30..0dcfeec74692af74b4418e34489b5b2a2aed769f 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -9,8 +9,8 @@ set (test_modules_sources
#testQGSJetII.cpp
testStackInspector.cpp
testSwitchProcess.cpp
- #testTrackingLine.cpp
- #testUrQMD.cpp
+ testTrackingLine.cpp
+ testUrQMD.cpp
)
add_executable (testModules ${test_modules_sources})
diff --git a/tests/modules/testSwitchProcess.cpp b/tests/modules/testSwitchProcess.cpp
index 510c731ff633331b8996326d3c940c68025b4796..0502bd3b9831993f6b2f3887665ff256837d1a3e 100644
--- a/tests/modules/testSwitchProcess.cpp
+++ b/tests/modules/testSwitchProcess.cpp
@@ -8,10 +8,10 @@
* the license.
*/
-#include <corsika/process/switch_process/SwitchProcess.h>
-#include <corsika/stack/SecondaryView.h>
-#include <corsika/stack/Stack.h>
-#include <corsika/units/PhysicalUnits.h>
+#include <corsika/modules/switch_process/SwitchProcess.hpp>
+#include <corsika/framework/stack/SecondaryView.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <catch2/catch.hpp>
@@ -19,7 +19,6 @@
#include <random>
using namespace corsika;
-using namespace corsika::process;
using namespace corsika::units::si;
class TestStackData {
@@ -61,11 +60,11 @@ private:
*/
template <typename StackIteratorInterface>
class TestParticleInterface
- : public corsika::stack::ParticleBase<StackIteratorInterface> {
+ : public corsika::ParticleBase<StackIteratorInterface> {
public:
- using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
- using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+ using corsika::ParticleBase<StackIteratorInterface>::GetStackData;
+ using corsika::ParticleBase<StackIteratorInterface>::GetIndex;
/*
The SetParticleData methods are called for creating new entries
@@ -85,13 +84,13 @@ public:
HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
};
-using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
+using SimpleStack = corsika::Stack<TestStackData, TestParticleInterface>;
// see issue 161
#if defined(__clang__)
-using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
+using StackTestView = corsika::SecondaryView<TestStackData, TestParticleInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
-using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
+using StackTestView = corsika::MakeView<SimpleStack>::type;
#endif
auto constexpr kgMSq = 1_kg / (1_m * 1_m);
@@ -106,7 +105,7 @@ struct DummyProcess : InteractionProcess<DummyProcess<N>> {
}
template <typename TSecondaries>
- corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
+ corsika::EProcessReturn DoInteraction(TSecondaries& vSec) {
// to figure out which process was selected in the end, we produce N
// secondaries for DummyProcess<N>
diff --git a/tests/modules/testTrackingLine.cpp b/tests/modules/testTrackingLine.cpp
index 638e43bc7f83f8f45d1446b75dbfdb1d5092dee2..1d6c8b8a20d89dc4969f77a288e8515e4141b1c8 100644
--- a/tests/modules/testTrackingLine.cpp
+++ b/tests/modules/testTrackingLine.cpp
@@ -8,32 +8,31 @@
* the license.
*/
-#include <corsika/process/tracking_line/TrackingLine.h>
-#include <testTrackingLineStack.h> // test-build, and include file is obtained from CMAKE_CURRENT_SOURCE_DIR
+#include <corsika/modules/tracking_line/TrackingLine.hpp>
-#include <corsika/environment/Environment.h>
-#include <corsika/particles/ParticleProperties.h>
+#include <testTrackingLineStack.hpp> // test-build, and include file is obtained from CMAKE_CURRENT_SOURCE_DIR
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Vector.h>
+#include <corsika/media/Environment.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/setup/SetupTrajectory.hpp>
using corsika::setup::Trajectory;
#include <catch2/catch.hpp>
using namespace corsika;
-using namespace corsika::process;
using namespace corsika::units;
-using namespace corsika::geometry;
#include <iostream>
using namespace std;
using namespace corsika::units::si;
TEST_CASE("TrackingLine") {
- environment::Environment<environment::Empty> env; // dummy environment
+ corsika::Environment<corsika::Empty> env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
tracking_line::TrackingLine tracking;
@@ -66,16 +65,16 @@ TEST_CASE("TrackingLine") {
auto const radius = 20_m;
- auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
+ auto theMedium = corsika::Environment<corsika::Empty>::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
auto const* theMediumPtr = theMedium.get();
universe.AddChild(std::move(theMedium));
TestTrackingLineStack stack;
stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::MuPlus,
+ std::tuple<corsika::Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point, units::si::TimeType>{
+ corsika::Code::MuPlus,
1_GeV,
{cs, {0_GeV, 0_GeV, 1_GeV}},
{cs, {0_m, 0_m, 0_km}},
diff --git a/tests/modules/testTrackingLineStack.hpp b/tests/modules/testTrackingLineStack.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..411178f5f670f08b3ad5020dbc60509a921dd01b
--- /dev/null
+++ b/tests/modules/testTrackingLineStack.hpp
@@ -0,0 +1,33 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/setup/SetupStack.hpp>
+
+using TestEnvironmentType = corsika::Environment<corsika::Empty>;
+
+template <typename T>
+using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>;
+
+// combine particle data stack with geometry information for tracking
+template <typename StackIter>
+using StackWithGeometryInterface =
+ corsika::CombinedParticleInterface<corsika::setup::detail::ParticleDataStack::PIType,
+ SetupGeometryDataInterface, StackIter>;
+using TestTrackingLineStack =
+ corsika::CombinedStack<typename corsika::setup::detail::ParticleDataStack::StackImpl,
+ GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
+
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 3339fb9d32cc84939a3d9c4b0bbb6dbc031c4012..a76025124840dca8f18186b906e292e9b13c67ef 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -8,25 +8,25 @@
* the license.
*/
-#include <corsika/process/urqmd/UrQMD.h>
-#include <corsika/random/RNGManager.h>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/modules/urqmd/UrQMD.hpp>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/units/PhysicalConstants.h>
-#include <corsika/units/PhysicalUnits.h>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/NuclearComposition.hpp>
#include <tuple>
#include <utility>
@@ -34,43 +34,43 @@
#include <catch2/catch.hpp>
using namespace corsika;
-using namespace corsika::process::UrQMD;
+using namespace corsika::urqmd;
using namespace corsika::units::si;
template <typename TStackView>
auto sumCharge(TStackView const& view) {
int totalCharge = 0;
- for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); }
+ for (auto const& p : view) { totalCharge += corsika::GetChargeNumber(p.GetPID()); }
return totalCharge;
}
template <typename TStackView>
-auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
- geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
+auto sumMomentum(TStackView const& view, corsika::CoordinateSystem const& vCS) {
+ corsika::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
for (auto const& p : view) { sum += p.GetMomentum(); }
return sum;
}
-auto setupEnvironment(particles::Code vTargetCode) {
+auto setupEnvironment(corsika::Code vTargetCode) {
// setup environment, geometry
- auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
+ auto env = std::make_unique<corsika::Environment<corsika::IMediumModel>>();
auto& universe = *(env->GetUniverse());
- const geometry::CoordinateSystem& cs = env->GetCoordinateSystem();
+ const corsika::CoordinateSystem& cs = env->GetCoordinateSystem();
auto theMedium =
- environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
- geometry::Point{cs, 0_m, 0_m, 0_m},
+ corsika::Environment<corsika::IMediumModel>::CreateNode<corsika::Sphere>(
+ corsika::Point{cs, 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+ using MyHomogeneousModel = corsika::HomogeneousMedium<corsika::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(std::vector<particles::Code>{vTargetCode},
- std::vector<float>{1.}));
+ corsika::NuclearComposition(std::vector<corsika::Code>{vTargetCode},
+ std::vector<float>{1.}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
@@ -80,102 +80,98 @@ auto setupEnvironment(particles::Code vTargetCode) {
template <typename TNodeType>
auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr,
- geometry::CoordinateSystem const& cs) {
+ corsika::CoordinateSystem const& cs) {
auto stack = std::make_unique<setup::Stack>();
auto constexpr mN = corsika::units::constants::nucleonMass;
- geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+ corsika::Point const origin(cs, {0_m, 0_m, 0_m});
+ corsika::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
sqrt(units::si::detail::static_pow<2>(mN * vA) + pLab.squaredNorm());
- 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, origin, 0_ns, vA, vZ});
+ auto particle = stack->AddParticle(
+ std::tuple<corsika::Code, units::si::HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, units::si::TimeType, unsigned short, unsigned short>{
+ corsika::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ});
particle.SetNode(vNodePtr);
return std::make_tuple(
std::move(stack),
- std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
+ std::make_unique<decltype(corsika::SecondaryView(particle))>(particle));
}
template <typename TNodeType>
-auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum,
- TNodeType* vNodePtr, geometry::CoordinateSystem const& cs) {
+auto setupStack(corsika::Code vProjectileType, HEPEnergyType vMomentum,
+ TNodeType* vNodePtr, corsika::CoordinateSystem const& cs) {
auto stack = std::make_unique<setup::Stack>();
- geometry::Point const origin(cs, {0_m, 0_m, 0_m});
- corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+ corsika::Point const origin(cs, {0_m, 0_m, 0_m});
+ corsika::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
- sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) +
+ sqrt(units::si::detail::static_pow<2>(corsika::GetMass(vProjectileType)) +
pLab.squaredNorm());
auto particle = stack->AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- vProjectileType, E0, pLab, origin, 0_ns});
+ std::tuple<corsika::Code, units::si::HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, units::si::TimeType>{vProjectileType, E0, pLab, origin,
+ 0_ns});
particle.SetNode(vNodePtr);
return std::make_tuple(
std::move(stack),
- std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
+ std::make_unique<decltype(corsika::SecondaryView(particle))>(particle));
}
TEST_CASE("UrQMD") {
SECTION("conversion") {
- REQUIRE_THROWS(process::UrQMD::ConvertFromUrQMD(106, 0));
- REQUIRE(process::UrQMD::ConvertFromUrQMD(101, 0) == particles::Code::Pi0);
- REQUIRE(process::UrQMD::ConvertToUrQMD(particles::Code::PiPlus) ==
+ REQUIRE_THROWS(corsika::urqmd::ConvertFromUrQMD(106, 0));
+ REQUIRE(corsika::urqmd::ConvertFromUrQMD(101, 0) == corsika::Code::Pi0);
+ REQUIRE(corsika::urqmd::ConvertToUrQMD(corsika::Code::PiPlus) ==
std::make_pair<int, int>(101, 2));
}
feenableexcept(FE_INVALID);
- corsika::random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
+ corsika::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
UrQMD urqmd;
SECTION("cross sections") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Unknown);
+ auto [env, csPtr, nodePtr] = setupEnvironment(corsika::Code::Unknown);
auto const& cs = *csPtr;
- particles::Code validProjectileCodes[] = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Proton,
- particles::Code::Neutron, particles::Code::KPlus, particles::Code::KMinus,
- particles::Code::K0, particles::Code::K0Bar, particles::Code::K0Long};
+ corsika::Code validProjectileCodes[] = {
+ corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::Proton,
+ corsika::Code::Neutron, corsika::Code::KPlus, corsika::Code::KMinus,
+ corsika::Code::K0, corsika::Code::K0Bar, corsika::Code::K0Long};
for (auto code : validProjectileCodes) {
auto [stack, view] = setupStack(code, 100_GeV, nodePtr, cs);
REQUIRE(stack->GetSize() == 1);
// simple check whether the cross-section is non-vanishing
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Proton) /
- 1_mb >
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Proton) / 1_mb >
0);
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Nitrogen) /
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Nitrogen) /
1_mb >
0);
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Oxygen) /
- 1_mb >
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Oxygen) / 1_mb >
0);
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Argon) /
- 1_mb >
+ REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Argon) / 1_mb >
0);
}
}
SECTION("nucleon projectile") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] = setupEnvironment(corsika::Code::Oxygen);
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
auto const projectileMomentum = projectile.GetMomentum();
- [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+ [[maybe_unused]] corsika::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
- Z + particles::GetChargeNumber(particles::Code::Oxygen));
+ Z + corsika::GetChargeNumber(corsika::Code::Oxygen));
auto const secMomSum =
sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
@@ -184,19 +180,19 @@ TEST_CASE("UrQMD") {
}
SECTION("\"special\" projectile") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] = setupEnvironment(corsika::Code::Oxygen);
auto [stackPtr, secViewPtr] =
- setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
+ setupStack(corsika::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
auto const projectileMomentum = projectile.GetMomentum();
- [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+ [[maybe_unused]] corsika::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
- particles::GetChargeNumber(particles::Code::PiPlus) +
- particles::GetChargeNumber(particles::Code::Oxygen));
+ corsika::GetChargeNumber(corsika::Code::PiPlus) +
+ corsika::GetChargeNumber(corsika::Code::Oxygen));
auto const secMomSum =
sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
@@ -205,19 +201,19 @@ TEST_CASE("UrQMD") {
}
SECTION("K0Long projectile") {
- auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] = setupEnvironment(corsika::Code::Oxygen);
auto [stackPtr, secViewPtr] =
- setupStack(particles::Code::K0Long, 400_GeV, nodePtr, *csPtr);
+ setupStack(corsika::Code::K0Long, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
auto const projectileMomentum = projectile.GetMomentum();
- [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+ [[maybe_unused]] corsika::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
- particles::GetChargeNumber(particles::Code::K0Long) +
- particles::GetChargeNumber(particles::Code::Oxygen));
+ corsika::GetChargeNumber(corsika::Code::K0Long) +
+ corsika::GetChargeNumber(corsika::Code::Oxygen));
auto const secMomSum =
sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());