diff --git a/corsika/detail/framework/utility/CorsikaData.inl b/corsika/detail/framework/utility/CorsikaData.inl
index e2ea1e99bc3742783020312f725419fbd0d63846..fc638b81a7f2e1a7c66908835d06182d99e69610 100644
--- a/corsika/detail/framework/utility/CorsikaData.inl
+++ b/corsika/detail/framework/utility/CorsikaData.inl
@@ -7,14 +7,15 @@
*/
#pragma once
+#include <boost/filesystem/path.hpp>
+
#include <cstdlib>
#include <stdexcept>
#include <string>
-inline std::string corsika::corsika_data(std::string const& key) {
+inline boost::filesystem::path corsika::corsika_data(boost::filesystem::path const& key) {
if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
- auto const path = std::string(p) + "/" + key;
- return path;
+ return boost::filesystem::path(p) / key;
} else {
throw std::runtime_error("CORSIKA_DATA not set");
}
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
index 828b9fcfeba8ae193ac622249a74ef317eaacc13..217fe30b72c3fa27d1ad3794a4a6b1ec8852ffd5 100644
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ b/corsika/detail/modules/qgsjetII/Interaction.inl
@@ -239,10 +239,14 @@ namespace corsika::qgsjetII {
if (is_nucleus(targetCode)) { // nucleus
targetMassNumber = get_nucleus_A(targetCode);
if (targetMassNumber > maxMassNumber_)
- throw std::runtime_error("QgsjetII target mass outside range.");
+ throw std::runtime_error(
+ "QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
} else {
- if (targetCode != Proton::code)
- throw std::runtime_error("QgsjetII Taget not possible.");
+ if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
+ throw std::runtime_error(
+ "QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
+ // here
}
CORSIKA_LOG_DEBUG("Interaction: target qgsjetII code/A: {}", targetMassNumber);
@@ -252,7 +256,9 @@ namespace corsika::qgsjetII {
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
projectileMassNumber = projectile.getNuclearA();
if (projectileMassNumber > maxMassNumber_)
- throw std::runtime_error("QgsjetII projectile mass outside range.");
+ throw std::runtime_error(
+ "QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1);
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index a1abaa6abf994b7f07a165b4f42599c3b4aaa89f..ee29642e00bbd408c355a5c2c645abef312564ec 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -15,15 +15,112 @@
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <boost/filesystem.hpp>
+#include <boost/multi_array.hpp>
+
#include <algorithm>
+#include <cassert>
#include <functional>
#include <iostream>
+#include <fstream>
+#include <sstream>
#include <urqmd.hpp>
namespace corsika::urqmd {
- inline UrQMD::UrQMD() { ::urqmd::iniurqmdc8_(); }
+ inline UrQMD::UrQMD(boost::filesystem::path const& xs_file) {
+ readXSFile(xs_file);
+ ::urqmd::iniurqmdc8_();
+ }
+
+ inline CrossSectionType UrQMD::getTabulatedCrossSection(Code projectileCode,
+ Code targetCode,
+ HEPEnergyType labEnergy) const {
+ // translated to C++ from CORSIKA 7 subroutine cxtot_u
+
+ auto const kinEnergy = labEnergy - get_mass(projectileCode);
+
+ assert(kinEnergy >= HEPEnergyType::zero());
+
+ double const logKinEnergy = std::log10(kinEnergy * (1 / 1_GeV));
+ double const ye = std::max(10 * logKinEnergy + 10.5, 1.);
+ int const je = std::min(int(ye), int(xs_interp_support_table_.shape()[2] - 2));
+ std::array<double, 3> w;
+ w[2 - 1] = ye - je;
+ w[3 - 1] = w[2 - 1] * (w[2 - 1] - 1.) * .5;
+ w[1 - 1] = 1 - w[2 - 1] + w[3 - 1];
+ w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
+
+ int projectileIndex;
+ switch (projectileCode) {
+ case Code::Proton:
+ projectileIndex = 0;
+ break;
+ case Code::AntiProton:
+ projectileIndex = 1;
+ break;
+ case Code::Neutron:
+ projectileIndex = 2;
+ break;
+ case Code::AntiNeutron:
+ projectileIndex = 3;
+ break;
+ case Code::PiPlus:
+ projectileIndex = 4;
+ break;
+ case Code::PiMinus:
+ projectileIndex = 5;
+ break;
+ case Code::KPlus:
+ projectileIndex = 6;
+ break;
+ case Code::KMinus:
+ projectileIndex = 7;
+ break;
+ case Code::K0Short:
+ case Code::K0Long:
+ /* since K0Short and K0Long are treated the same, we can also add K0 and K0Bar
+ * to the list. This is a deviation from CORSIKA 7. */
+ case Code::K0:
+ case Code::K0Bar:
+ projectileIndex = 8;
+ break;
+ default: { // LCOV_EXCL_START since this can never happen due to canInteract
+ CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
+ return CrossSectionType::zero(); // LCOV_EXCL_STOP
+ }
+ }
+
+ int targetIndex;
+ switch (targetCode) {
+ case Code::Nitrogen:
+ targetIndex = 0;
+ break;
+ case Code::Oxygen:
+ targetIndex = 1;
+ break;
+ case Code::Argon:
+ targetIndex = 2;
+ break;
+ default:
+ std::stringstream ss;
+ ss << "UrQMD cross-section not tabluated for target " << targetCode;
+ throw std::runtime_error(ss.str().data());
+ }
+
+ auto result = CrossSectionType::zero();
+ for (int i = 0; i < 3; ++i) {
+ result +=
+ xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i];
+ }
+
+ CORSIKA_LOG_TRACE(
+ "UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
+ get_name(projectileCode), get_name(targetCode), labEnergy / 1_GeV, result);
+
+ return result;
+ }
inline CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
HEPEnergyType vLabEnergy,
@@ -63,20 +160,19 @@ namespace corsika::urqmd {
template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile
- inline CrossSectionType UrQMD::getCrossSection(TParticle const& vProjectile,
- Code vTargetCode) const {
- // TODO: return 0 for non-hadrons?
-
- auto const projectileCode = vProjectile.getPID();
- auto const projectileEnergyLab = vProjectile.getEnergy();
-
- if (projectileCode == Code::K0Long) {
- return 0.5 * (getCrossSection(Code::K0, vTargetCode, projectileEnergyLab) +
- getCrossSection(Code::K0Bar, vTargetCode, projectileEnergyLab));
+ inline CrossSectionType UrQMD::getCrossSection(TParticle const& projectile,
+ Code targetCode) const {
+ auto const projectileCode = projectile.getPID();
+
+ if (projectileCode == Code::Nucleus) {
+ /*
+ * unfortunately unavoidable at the moment until we have tools to get the actual
+ * inealstic cross-section from UrQMD
+ */
+ return CrossSectionType::zero();
}
- int const Ap = (projectileCode == Code::Nucleus) ? vProjectile.getNuclearA() : 1;
- return getCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
+ return getTabulatedCrossSection(projectileCode, targetCode, projectile.getEnergy());
}
inline bool UrQMD::canInteract(Code vCode) const {
@@ -85,10 +181,9 @@ namespace corsika::urqmd {
// only the usual long-lived species as input.
// TODO: Charmed mesons should be added to the list, too
- static Code const validProjectileCodes[] = {
- Code::Nucleus, Code::Proton, Code::AntiProton, Code::Neutron,
- Code::AntiNeutron, Code::PiPlus, Code::PiMinus, Code::KPlus,
- Code::KMinus, Code::K0, Code::K0Bar, Code::K0Long};
+ static std::array constexpr validProjectileCodes{
+ Code::Proton, Code::AntiProton, Code::Neutron, Code::AntiNeutron, Code::PiPlus,
+ Code::PiMinus, Code::KPlus, Code::KMinus, Code::K0Short, Code::K0Long};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
@@ -301,4 +396,41 @@ namespace corsika::urqmd {
return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
}
+ inline void UrQMD::readXSFile(boost::filesystem::path const& filename) {
+ boost::filesystem::ifstream file(filename, std::ios::in);
+
+ if (!file.is_open()) {
+ throw std::runtime_error(
+ filename.native() +
+ " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
+ }
+
+ std::string line;
+
+ std::getline(file, line);
+ std::stringstream ss(line);
+
+ char dummy;
+ int nTargets, nProjectiles, nSupports;
+ ss >> dummy >> nTargets >> nProjectiles >> nSupports;
+
+ decltype(xs_interp_support_table_)::extent_gen extents;
+ xs_interp_support_table_.resize(extents[nProjectiles][nTargets][nSupports]);
+
+ for (int i = 0; i < nTargets; ++i) {
+ for (int j = 0; j < nProjectiles; ++j) {
+ for (int k = 0; k < nSupports; ++k) {
+ std::getline(file, line);
+ std::stringstream s(line);
+ double energy, sigma;
+ s >> energy >> sigma;
+ xs_interp_support_table_[j][i][k] = sigma * 1_mb;
+ }
+
+ std::getline(file, line);
+ std::getline(file, line);
+ }
+ }
+ }
+
} // namespace corsika::urqmd
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index befd3025035105a184f33c1c2a872c750364665e..43a123cf18fa62fe62baded114b319525b856d55 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -19,14 +19,36 @@
namespace corsika {
/**
+ @defgroup Utilities
+
+ Collection of classes and methods to perform recurring tasks.
+ **/
+
+ /**
+ @class
+ @ingroup Utilities
+
This utility class handles Lorentz boost between different
referenence frames, using FourVector.
+
+ The class is initialized with projectile and optionally target
+ energy/momentum data. During initialization, a rotation matrix is
+ calculated to represent the projectile movement (and thus the
+ boost) along the z-axis. Also the inverse of this rotation is
+ calculated. The Lorentz boost matrix and its inverse are
+ determined as 2x2 matrices considering the energy and
+ pz-momentum.
+
+ Different constructors are offered with different specialization
+ for the cases of collisions (projectile-target) or just decays
+ (projectile only).
*/
class COMBoost {
public:
- //! construct a COMBoost given four-vector of projectile and mass of target
+ //! construct a COMBoost given four-vector of projectile and mass of target (target at
+ //! rest)
COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
HEPEnergyType const massTarget);
@@ -41,9 +63,11 @@ namespace corsika {
template <typename FourVector>
FourVector fromCoM(FourVector const& p) const;
+ //! returns the rotated coordinate system
CoordinateSystemPtr getRotatedCS() const;
protected:
+ //! internal method
void setBoost(double coshEta, double sinhEta);
private:
diff --git a/corsika/framework/utility/CorsikaData.hpp b/corsika/framework/utility/CorsikaData.hpp
index d117db8da59e3399e77ff7304335c7e9eaa21711..2ab2dec7dabcf7a81bcbe976d689bd92f72df1b3 100644
--- a/corsika/framework/utility/CorsikaData.hpp
+++ b/corsika/framework/utility/CorsikaData.hpp
@@ -8,13 +8,18 @@
#pragma once
-#include <string>
+#include <boost/filesystem/path.hpp>
namespace corsika {
/**
+ * @file CorsikaData.hpp
+ * @ingroup Utilities
+ * @{
* returns the full path of the file \p filename within the CORSIKA_DATA directory
*/
- std::string corsika_data(std::string const& filename);
+ boost::filesystem::path corsika_data(boost::filesystem::path const& filename);
+
+ //! @}
} // namespace corsika
#include <corsika/detail/framework/utility/CorsikaData.inl>
diff --git a/corsika/framework/utility/QuarticSolver.hpp b/corsika/framework/utility/QuarticSolver.hpp
index 0bcdfcedf1a9a9055306acacfd756012726cf8b0..ad9484939d10691e38c7e74d9e2632b3d2e60463 100644
--- a/corsika/framework/utility/QuarticSolver.hpp
+++ b/corsika/framework/utility/QuarticSolver.hpp
@@ -11,6 +11,9 @@
#include <complex>
/**
+ * @file QuarticSolver.hpp
+ * @ingroup Utilities
+ * @{
* \todo convert to class
*/
@@ -56,6 +59,8 @@ namespace corsika::quartic_solver {
// afterwards.
DComplex* solve_quartic(double a, double b, double c, double d);
+ //! @}
+
} // namespace corsika::quartic_solver
#include <corsika/detail/framework/utility/QuarticSolver.inl>
diff --git a/corsika/framework/utility/SaveBoostHistogram.hpp b/corsika/framework/utility/SaveBoostHistogram.hpp
index 14e0379966c447cd0c2f8cfa7ccffcf3c29a0478..30cbdda2a6c9e95ce63dbfbaf7c7a120d4bbe727 100644
--- a/corsika/framework/utility/SaveBoostHistogram.hpp
+++ b/corsika/framework/utility/SaveBoostHistogram.hpp
@@ -13,6 +13,8 @@
namespace corsika {
/**
+ * @ingroup Utilities
+ *
* This functions saves a boost::histogram into a numpy file. Only rather basic axis
* types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
* counts (i.e. a double or int) are supported, nothing fancy like profiles.
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 8f2fd8edeed26f84be2b6049fb760c26b4729ac4..3f6d3e9936f1fd13fcb831359f8faee6768f2948 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -12,21 +12,28 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
#include <corsika/setup/SetupStack.hpp>
+#include <boost/filesystem/path.hpp>
+#include <boost/multi_array.hpp>
+
#include <array>
#include <utility>
+#include <string>
namespace corsika::urqmd {
class UrQMD : public InteractionProcess<UrQMD> {
public:
- UrQMD();
+ UrQMD(boost::filesystem::path const& path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"));
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
+ CrossSectionType getTabulatedCrossSection(Code, Code, HEPEnergyType) const;
+
template <typename TParticle>
CrossSectionType getCrossSection(TParticle const&, Code) const;
@@ -39,11 +46,12 @@ namespace corsika::urqmd {
private:
static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int);
+ void readXSFile(boost::filesystem::path const&);
// data members
default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("urqmd");
-
std::uniform_int_distribution<int> booleanDist_{0, 1};
+ boost::multi_array<CrossSectionType, 3> xs_interp_support_table_;
};
/**
diff --git a/documentation/index.rst b/documentation/index.rst
index c1be80a37ee4a481cd3c67638772b11309c5a5d3..79440c503bd0b500d8235d5c1edd77c2bf459822 100644
--- a/documentation/index.rst
+++ b/documentation/index.rst
@@ -13,6 +13,7 @@ Welcome to the CORSIKA 8 air shower simulation framework.
units
environment
stack
+ utilities
api
diff --git a/documentation/utilities.rst b/documentation/utilities.rst
new file mode 100644
index 0000000000000000000000000000000000000000..37c289fde6a68f36548ff2fcdab0ce0dbda2114b
--- /dev/null
+++ b/documentation/utilities.rst
@@ -0,0 +1,9 @@
+Utilities
+==========
+
+.. doxygengroup:: Utilities
+ :project: CORSIKA8
+ :members:
+
+
+
diff --git a/modules/urqmd/urqmdInterface.F b/modules/urqmd/urqmdInterface.F
index ec2bd092a8386393368a6789df73b3efa3253160..446221dfc2b2f311b11a4c0080e42ba02fab7296 100644
--- a/modules/urqmd/urqmdInterface.F
+++ b/modules/urqmd/urqmdInterface.F
@@ -432,16 +432,3 @@ c~ xsegymin=0.25d0
end
-c
-c M. Reininghaus, 2020-04-08
-c
- integer function ReadSigmaLn(ia, ib, ic)
- implicit none
-
- include 'comres.f'
-
- integer :: ia, ib, ic
-
- ReadSigmaLn = SigmaLn(ia, ib, ic)
-
- end function ReadSigmaLn
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 477dedc6ca15a123b26f5524cbeb888fb2700cd7..4a604fd8bbf32f440db2eb0e073aed684b412fe4 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -175,7 +175,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
SECTION("InteractionInterface Nuclei") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Nucleus, 20, 10, 10100_GeV,
+ Code::Nucleus, 60, 30, 20100_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
@@ -183,10 +183,15 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto const projectileMomentum = projectile.getMomentum();
corsika::qgsjetII::Interaction model;
- model.doInteraction(view);
+ model.doInteraction(view); // this also should produce some fragments
+ CHECK(view.getSize() == Approx(188).margin(2)); // this is not physics validation
+ int countFragments = 0;
+ for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); }
+ CHECK(countFragments == Approx(2).margin(1)); // this is not physics validation
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
- CHECK(length / (1_g / square(1_cm)) == Approx(20.13).margin(0.1));
+ CHECK(length / (1_g / square(1_cm)) ==
+ Approx(12).margin(2)); // this is not physics validation
}
SECTION("Heavy nuclei") {
@@ -210,17 +215,50 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
}
SECTION("Allowed Particles") {
-
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Electron, 0, 0, 1100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
- auto particle = stackPtr->first();
- auto projectile = secViewPtr->getProjectile();
- auto const projectileMomentum = projectile.getMomentum();
-
- corsika::qgsjetII::Interaction model;
-
- GrammageType const length = model.getInteractionLength(particle);
- CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
+ { // electron
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ auto particle = stackPtr->first();
+ corsika::qgsjetII::Interaction model;
+ GrammageType const length = model.getInteractionLength(particle);
+ CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
+ }
+ { // pi0 is internally converted into pi+/pi-
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(18).margin(2)); // this is not physics validation
+ }
+ { // rho0 is internally converted into pi-/pi+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(7).margin(2)); // this is not physics validation
+ }
+ { // Lambda is internally converted into neutron
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Lambda0, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ }
+ { // AntiLambda is internally converted into anti neutron
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Lambda0Bar, 0, 0, 100_GeV,
+ (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ }
}
}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 28ed39a159c9ceb8bb0960d8adac735fc67a95a2..0f9f4a86243ec0f7f5b0c16e989b9459a5e0a293 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -77,24 +77,41 @@ TEST_CASE("UrQMD") {
auto const& cs = *csPtr;
{ [[maybe_unused]] auto const& env_dummy = env; }
- Code validProjectileCodes[] = {Code::PiPlus, Code::PiMinus, Code::Proton,
- Code::Neutron, Code::KPlus, Code::KMinus,
- Code::K0, Code::K0Bar, Code::K0Long};
+ Code validProjectileCodes[] = {Code::PiPlus, Code::PiMinus, Code::Proton,
+ Code::AntiProton, Code::AntiNeutron, Code::Neutron,
+ Code::KPlus, Code::KMinus, Code::K0,
+ Code::K0Bar, Code::K0Long};
for (auto code : validProjectileCodes) {
- auto [stack, view] = setup::testing::setup_stack(
- code, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ auto [stack, view] = setup::testing::setup_stack(code, 0, 0, 100_GeV, nodePtr, cs);
CHECK(stack->getEntries() == 1);
CHECK(view->getEntries() == 0);
// simple check whether the cross-section is non-vanishing
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Proton) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Nitrogen) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Oxygen) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Argon) / 1_mb > 0);
+ // only nuclei with available tabluated data so far
+ CHECK(urqmd.getInteractionLength(stack->getNextParticle()) > 1_g / square(1_cm));
}
}
+ SECTION("targets options") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
+ auto const& cs = *csPtr;
+ { [[maybe_unused]] auto const& env_dummy = env; }
+ auto [stack, view] =
+ setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs);
+ CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) ==
+ Approx(105).margin(5));
+ }
+
+ SECTION("invalid targets options") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Omega);
+ auto const& cs = *csPtr;
+ { [[maybe_unused]] auto const& env_dummy = env; }
+ auto [stack, view] =
+ setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs);
+ CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
+ }
+
SECTION("nucleus projectile") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
@@ -146,6 +163,45 @@ TEST_CASE("UrQMD") {
Approx(0).margin(1e-2));
}
+ SECTION("\"special\" projectile and target") {
+ {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target
+ }
+
+ {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK(stackPtr->getEntries() == 1);
+ CHECK(secViewPtr->getEntries() == 0);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
+
+ urqmd.doInteraction(*secViewPtr);
+
+ CHECK(sumCharge(*secViewPtr) ==
+ get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+ Approx(0).margin(1e-2));
+ }
+ }
+
SECTION("K0Long projectile") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
@@ -171,4 +227,22 @@ TEST_CASE("UrQMD") {
CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
Approx(0).margin(1e-2));
}
+
+ SECTION("K0Short projectile") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::K0Short, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK(stackPtr->getEntries() == 1);
+ CHECK(secViewPtr->getEntries() == 0);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
+
+ CHECK_THROWS(urqmd.doInteraction(*secViewPtr));
+ }
}