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 e85a6abc302824e80f1188e57d7fda4cd56216a8..ee29642e00bbd408c355a5c2c645abef312564ec 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -86,9 +86,9 @@ namespace corsika::urqmd {
case Code::K0Bar:
projectileIndex = 8;
break;
- default: {
+ 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();
+ return CrossSectionType::zero(); // LCOV_EXCL_STOP
}
}
@@ -400,7 +400,9 @@ namespace corsika::urqmd {
boost::filesystem::ifstream file(filename, std::ios::in);
if (!file.is_open()) {
- throw std::runtime_error(filename.native() + " could not be opened.");
+ throw std::runtime_error(
+ filename.native() +
+ " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
}
std::string line;
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 9f24d5036a34f96ab10a3422c96b735b6338f0df..2ab2dec7dabcf7a81bcbe976d689bd92f72df1b3 100644
--- a/corsika/framework/utility/CorsikaData.hpp
+++ b/corsika/framework/utility/CorsikaData.hpp
@@ -12,9 +12,14 @@
namespace corsika {
/**
+ * @file CorsikaData.hpp
+ * @ingroup Utilities
+ * @{
* returns the full path of the file \p filename within the CORSIKA_DATA directory
*/
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/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/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 92a544932db254827fb2cd70518214d00b4d7db6..da788b797ffc8785b2fdf0275075d33d28b1abc8 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -84,15 +84,34 @@ TEST_CASE("UrQMD") {
for (auto code : validProjectileCodes) {
auto [stack, view] = setup::testing::setup_stack(code, 0, 0, 100_GeV, nodePtr, cs);
- REQUIRE(stack->getEntries() == 1);
- REQUIRE(view->getEntries() == 0);
+ CHECK(stack->getEntries() == 1);
+ CHECK(view->getEntries() == 0);
// simple check whether the cross-section is non-vanishing
// only nuclei with available tabluated data so far
- REQUIRE(urqmd.getInteractionLength(stack->getNextParticle()) > 1_g / square(1_cm));
+ 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
@@ -144,6 +163,32 @@ 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(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
@@ -169,4 +214,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));
+ }
}