From 9f89b5c03d97346b1180cc4b822827222a3150a3 Mon Sep 17 00:00:00 2001 From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu> Date: Tue, 17 Nov 2020 13:19:41 +0100 Subject: [PATCH] QGSJetII --- .../detail/modules/qgsjetII/Interaction.inl | 30 +++++++++---------- corsika/modules/qgsjetII/Interaction.hpp | 4 +-- tests/modules/testQGSJetII.cpp | 18 +++++------ 3 files changed, 25 insertions(+), 27 deletions(-) diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl index 3900c95d6..0f5a55591 100644 --- a/corsika/detail/modules/qgsjetII/Interaction.inl +++ b/corsika/detail/modules/qgsjetII/Interaction.inl @@ -67,7 +67,7 @@ namespace corsika::qgsjetII { const int iBeam = corsika::qgsjetII::GetQgsjetIIXSCode(beamId); int iTarget = 1; - if (corsika::IsNucleus(targetId)) { + if (corsika::is_nucleus(targetId)) { iTarget = targetA; if (iTarget > maxMassNumber_ || iTarget <= 0) { std::ostringstream txt; @@ -76,7 +76,7 @@ namespace corsika::qgsjetII { } } int iProjectile = 1; - if (corsika::IsNucleus(beamId)) { + if (corsika::is_nucleus(beamId)) { iProjectile = Abeam; if (iProjectile > maxMassNumber_ || iProjectile <= 0) throw std::runtime_error("QgsjetII target outside range. "); @@ -118,7 +118,7 @@ namespace corsika::qgsjetII { if (kInteraction) { int Abeam = 0; - if (corsika::IsNucleus(vP.GetPID())) Abeam = vP.GetNuclearA(); + if (corsika::is_nucleus(vP.GetPID())) Abeam = vP.GetNuclearA(); // get target from environment /* @@ -134,7 +134,7 @@ namespace corsika::qgsjetII { CrossSectionType weightedProdCrossSection = mediumComposition.WeightedSum([=](corsika::Code targetID) -> CrossSectionType { int targetA = 0; - if (corsika::IsNucleus(targetID)) targetA = corsika::GetNucleusA(targetID); + if (corsika::is_nucleus(targetID)) targetA = corsika::nucleus_A(targetID); return GetCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA); }); @@ -189,7 +189,7 @@ namespace corsika::qgsjetII { auto const projectileMomentumLab = vP.GetMomentum(); int beamA = 0; - if (corsika::IsNucleus(corsikaBeamId)) beamA = vP.GetNuclearA(); + if (corsika::is_nucleus(corsikaBeamId)) beamA = vP.GetNuclearA(); std::cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << std::endl << "Interaction: pbeam lab: " @@ -217,7 +217,7 @@ namespace corsika::qgsjetII { for (size_t i = 0; i < compVec.size(); ++i) { auto const targetId = compVec[i]; int targetA = 0; - if (corsika::IsNucleus(targetId)) targetA = corsika::GetNucleusA(targetId); + if (corsika::is_nucleus(targetId)) targetA = corsika::nucleus_A(targetId); const auto sigProd = GetCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA); cross_section_of_components[i] = sigProd; @@ -228,15 +228,14 @@ namespace corsika::qgsjetII { std::cout << "Interaction: target selected: " << targetCode << std::endl; int targetQgsCode = -1; - if (corsika::IsNucleus(targetCode)) - targetQgsCode = corsika::GetNucleusA(targetCode); - if (targetCode == corsika::Proton::GetCode()) targetQgsCode = 1; + if (corsika::is_nucleus(targetCode)) targetQgsCode = corsika::nucleus_A(targetCode); + if (targetCode == corsika::Code::Proton) targetQgsCode = 1; std::cout << "Interaction: target qgsjetII code/A: " << targetQgsCode << std::endl; if (targetQgsCode > maxMassNumber_ || targetQgsCode < 1) throw std::runtime_error("QgsjetII target outside range."); int projQgsCode = 1; - if (corsika::IsNucleus(corsikaBeamId)) projQgsCode = vP.GetNuclearA(); + if (corsika::is_nucleus(corsikaBeamId)) projQgsCode = vP.GetNuclearA(); std::cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " " << corsikaBeamId << std::endl; if (projQgsCode > maxMassNumber_ || projQgsCode < 1) @@ -290,14 +289,13 @@ namespace corsika::qgsjetII { idFragm = corsika::Code::Proton; auto momentum = corsika::Vector( - zAxisFrame, - corsika::QuantityVector<hepmomentum_d>{ - 0.0_GeV, 0.0_GeV, - sqrt((projectileEnergyLab + corsika::Proton::GetMass()) * - (projectileEnergyLab - corsika::Proton::GetMass()))}); + zAxisFrame, corsika::QuantityVector<hepmomentum_d>{ + 0.0_GeV, 0.0_GeV, + sqrt((projectileEnergyLab + corsika::Proton::mass()) * + (projectileEnergyLab - corsika::Proton::mass()))}); auto const energy = - sqrt(momentum.squaredNorm() + square(corsika::GetMass(idFragm))); + sqrt(momentum.squaredNorm() + square(corsika::mass(idFragm))); momentum.rebase(originalCS); // transform back into standard lab frame std::cout << "secondary fragment> id=" << idFragm << " p=" << momentum.GetComponents() << std::endl; diff --git a/corsika/modules/qgsjetII/Interaction.hpp b/corsika/modules/qgsjetII/Interaction.hpp index 6d94312e4..f55bae8b2 100644 --- a/corsika/modules/qgsjetII/Interaction.hpp +++ b/corsika/modules/qgsjetII/Interaction.hpp @@ -35,8 +35,8 @@ namespace corsika::qgsjetII { bool WasInitialized() { return initialized_; } int GetMaxTargetMassNumber() const { return maxMassNumber_; } bool IsValidTarget(corsika::Code TargetId) const { - return (corsika::GetNucleusA(TargetId) < maxMassNumber_) && - corsika::IsNucleus(TargetId); + return corsika::is_nucleus(TargetId) && + (corsika::nucleus_A(TargetId) < maxMassNumber_); } CrossSectionType GetCrossSection(const corsika::Code, const corsika::Code, diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp index 5ecfe1d63..e145b73d1 100644 --- a/tests/modules/testQGSJetII.cpp +++ b/tests/modules/testQGSJetII.cpp @@ -22,25 +22,25 @@ using namespace corsika::qgsjetII; TEST_CASE("QgsjetII", "[processes]") { SECTION("QgsjetII -> Corsika") { - REQUIRE(corsika::PiPlus::GetCode() == corsika::qgsjetII::ConvertFromQgsjetII( - corsika::qgsjetII::QgsjetIICode::PiPlus)); + REQUIRE(corsika::Code::PiPlus == corsika::qgsjetII::ConvertFromQgsjetII( + corsika::qgsjetII::QgsjetIICode::PiPlus)); } SECTION("Corsika -> QgsjetII") { - REQUIRE(corsika::qgsjetII::ConvertToQgsjetII(corsika::PiMinus::GetCode()) == + REQUIRE(corsika::qgsjetII::ConvertToQgsjetII(corsika::Code::PiMinus) == corsika::qgsjetII::QgsjetIICode::PiMinus); - REQUIRE(corsika::qgsjetII::ConvertToQgsjetIIRaw(corsika::Proton::GetCode()) == 2); + REQUIRE(corsika::qgsjetII::ConvertToQgsjetIIRaw(corsika::Code::Proton) == 2); } SECTION("canInteractInQgsjetII") { - REQUIRE(corsika::qgsjetII::CanInteract(corsika::Proton::GetCode())); + REQUIRE(corsika::qgsjetII::CanInteract(corsika::Code::Proton)); REQUIRE(corsika::qgsjetII::CanInteract(corsika::Code::KPlus)); - REQUIRE(corsika::qgsjetII::CanInteract(corsika::Nucleus::GetCode())); + REQUIRE(corsika::qgsjetII::CanInteract(corsika::Code::Nucleus)); // REQUIRE(corsika::qgsjetII::CanInteract(corsika::Helium::GetCode())); - REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::EtaC::GetCode())); - REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::SigmaC0::GetCode())); + REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::Code::EtaC)); + REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::Code::SigmaC0)); } SECTION("cross-section type") { @@ -95,7 +95,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { setup::Stack stack; const HEPEnergyType E0 = 100_GeV; HEPMomentumType P0 = - sqrt(E0 * E0 - corsika::Proton::GetMass() * corsika::Proton::GetMass()); + sqrt(E0 * E0 - corsika::Proton::mass() * corsika::Proton::mass()); auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0}); corsika::Point pos(cs, 0_m, 0_m, 0_m); auto particle = stack.AddParticle( -- GitLab