IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 9f89b5c0 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan: Committed by ralfulrich
Browse files

QGSJetII

parent dadf106e
No related branches found
No related tags found
1 merge request!280Refactory 2020
...@@ -67,7 +67,7 @@ namespace corsika::qgsjetII { ...@@ -67,7 +67,7 @@ namespace corsika::qgsjetII {
const int iBeam = corsika::qgsjetII::GetQgsjetIIXSCode(beamId); const int iBeam = corsika::qgsjetII::GetQgsjetIIXSCode(beamId);
int iTarget = 1; int iTarget = 1;
if (corsika::IsNucleus(targetId)) { if (corsika::is_nucleus(targetId)) {
iTarget = targetA; iTarget = targetA;
if (iTarget > maxMassNumber_ || iTarget <= 0) { if (iTarget > maxMassNumber_ || iTarget <= 0) {
std::ostringstream txt; std::ostringstream txt;
...@@ -76,7 +76,7 @@ namespace corsika::qgsjetII { ...@@ -76,7 +76,7 @@ namespace corsika::qgsjetII {
} }
} }
int iProjectile = 1; int iProjectile = 1;
if (corsika::IsNucleus(beamId)) { if (corsika::is_nucleus(beamId)) {
iProjectile = Abeam; iProjectile = Abeam;
if (iProjectile > maxMassNumber_ || iProjectile <= 0) if (iProjectile > maxMassNumber_ || iProjectile <= 0)
throw std::runtime_error("QgsjetII target outside range. "); throw std::runtime_error("QgsjetII target outside range. ");
...@@ -118,7 +118,7 @@ namespace corsika::qgsjetII { ...@@ -118,7 +118,7 @@ namespace corsika::qgsjetII {
if (kInteraction) { if (kInteraction) {
int Abeam = 0; int Abeam = 0;
if (corsika::IsNucleus(vP.GetPID())) Abeam = vP.GetNuclearA(); if (corsika::is_nucleus(vP.GetPID())) Abeam = vP.GetNuclearA();
// get target from environment // get target from environment
/* /*
...@@ -134,7 +134,7 @@ namespace corsika::qgsjetII { ...@@ -134,7 +134,7 @@ namespace corsika::qgsjetII {
CrossSectionType weightedProdCrossSection = CrossSectionType weightedProdCrossSection =
mediumComposition.WeightedSum([=](corsika::Code targetID) -> CrossSectionType { mediumComposition.WeightedSum([=](corsika::Code targetID) -> CrossSectionType {
int targetA = 0; 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); return GetCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA);
}); });
...@@ -189,7 +189,7 @@ namespace corsika::qgsjetII { ...@@ -189,7 +189,7 @@ namespace corsika::qgsjetII {
auto const projectileMomentumLab = vP.GetMomentum(); auto const projectileMomentumLab = vP.GetMomentum();
int beamA = 0; 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 std::cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << std::endl
<< "Interaction: pbeam lab: " << "Interaction: pbeam lab: "
...@@ -217,7 +217,7 @@ namespace corsika::qgsjetII { ...@@ -217,7 +217,7 @@ namespace corsika::qgsjetII {
for (size_t i = 0; i < compVec.size(); ++i) { for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i]; auto const targetId = compVec[i];
int targetA = 0; int targetA = 0;
if (corsika::IsNucleus(targetId)) targetA = corsika::GetNucleusA(targetId); if (corsika::is_nucleus(targetId)) targetA = corsika::nucleus_A(targetId);
const auto sigProd = const auto sigProd =
GetCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA); GetCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA);
cross_section_of_components[i] = sigProd; cross_section_of_components[i] = sigProd;
...@@ -228,15 +228,14 @@ namespace corsika::qgsjetII { ...@@ -228,15 +228,14 @@ namespace corsika::qgsjetII {
std::cout << "Interaction: target selected: " << targetCode << std::endl; std::cout << "Interaction: target selected: " << targetCode << std::endl;
int targetQgsCode = -1; int targetQgsCode = -1;
if (corsika::IsNucleus(targetCode)) if (corsika::is_nucleus(targetCode)) targetQgsCode = corsika::nucleus_A(targetCode);
targetQgsCode = corsika::GetNucleusA(targetCode); if (targetCode == corsika::Code::Proton) targetQgsCode = 1;
if (targetCode == corsika::Proton::GetCode()) targetQgsCode = 1;
std::cout << "Interaction: target qgsjetII code/A: " << targetQgsCode << std::endl; std::cout << "Interaction: target qgsjetII code/A: " << targetQgsCode << std::endl;
if (targetQgsCode > maxMassNumber_ || targetQgsCode < 1) if (targetQgsCode > maxMassNumber_ || targetQgsCode < 1)
throw std::runtime_error("QgsjetII target outside range."); throw std::runtime_error("QgsjetII target outside range.");
int projQgsCode = 1; 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 << " " std::cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " "
<< corsikaBeamId << std::endl; << corsikaBeamId << std::endl;
if (projQgsCode > maxMassNumber_ || projQgsCode < 1) if (projQgsCode > maxMassNumber_ || projQgsCode < 1)
...@@ -290,14 +289,13 @@ namespace corsika::qgsjetII { ...@@ -290,14 +289,13 @@ namespace corsika::qgsjetII {
idFragm = corsika::Code::Proton; idFragm = corsika::Code::Proton;
auto momentum = corsika::Vector( auto momentum = corsika::Vector(
zAxisFrame, zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
corsika::QuantityVector<hepmomentum_d>{ 0.0_GeV, 0.0_GeV,
0.0_GeV, 0.0_GeV, sqrt((projectileEnergyLab + corsika::Proton::mass()) *
sqrt((projectileEnergyLab + corsika::Proton::GetMass()) * (projectileEnergyLab - corsika::Proton::mass()))});
(projectileEnergyLab - corsika::Proton::GetMass()))});
auto const energy = 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 momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << "secondary fragment> id=" << idFragm std::cout << "secondary fragment> id=" << idFragm
<< " p=" << momentum.GetComponents() << std::endl; << " p=" << momentum.GetComponents() << std::endl;
......
...@@ -35,8 +35,8 @@ namespace corsika::qgsjetII { ...@@ -35,8 +35,8 @@ namespace corsika::qgsjetII {
bool WasInitialized() { return initialized_; } bool WasInitialized() { return initialized_; }
int GetMaxTargetMassNumber() const { return maxMassNumber_; } int GetMaxTargetMassNumber() const { return maxMassNumber_; }
bool IsValidTarget(corsika::Code TargetId) const { bool IsValidTarget(corsika::Code TargetId) const {
return (corsika::GetNucleusA(TargetId) < maxMassNumber_) && return corsika::is_nucleus(TargetId) &&
corsika::IsNucleus(TargetId); (corsika::nucleus_A(TargetId) < maxMassNumber_);
} }
CrossSectionType GetCrossSection(const corsika::Code, const corsika::Code, CrossSectionType GetCrossSection(const corsika::Code, const corsika::Code,
......
...@@ -22,25 +22,25 @@ using namespace corsika::qgsjetII; ...@@ -22,25 +22,25 @@ using namespace corsika::qgsjetII;
TEST_CASE("QgsjetII", "[processes]") { TEST_CASE("QgsjetII", "[processes]") {
SECTION("QgsjetII -> Corsika") { SECTION("QgsjetII -> Corsika") {
REQUIRE(corsika::PiPlus::GetCode() == corsika::qgsjetII::ConvertFromQgsjetII( REQUIRE(corsika::Code::PiPlus == corsika::qgsjetII::ConvertFromQgsjetII(
corsika::qgsjetII::QgsjetIICode::PiPlus)); corsika::qgsjetII::QgsjetIICode::PiPlus));
} }
SECTION("Corsika -> QgsjetII") { SECTION("Corsika -> QgsjetII") {
REQUIRE(corsika::qgsjetII::ConvertToQgsjetII(corsika::PiMinus::GetCode()) == REQUIRE(corsika::qgsjetII::ConvertToQgsjetII(corsika::Code::PiMinus) ==
corsika::qgsjetII::QgsjetIICode::PiMinus); corsika::qgsjetII::QgsjetIICode::PiMinus);
REQUIRE(corsika::qgsjetII::ConvertToQgsjetIIRaw(corsika::Proton::GetCode()) == 2); REQUIRE(corsika::qgsjetII::ConvertToQgsjetIIRaw(corsika::Code::Proton) == 2);
} }
SECTION("canInteractInQgsjetII") { 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::Code::KPlus));
REQUIRE(corsika::qgsjetII::CanInteract(corsika::Nucleus::GetCode())); REQUIRE(corsika::qgsjetII::CanInteract(corsika::Code::Nucleus));
// REQUIRE(corsika::qgsjetII::CanInteract(corsika::Helium::GetCode())); // REQUIRE(corsika::qgsjetII::CanInteract(corsika::Helium::GetCode()));
REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::EtaC::GetCode())); REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::Code::EtaC));
REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::SigmaC0::GetCode())); REQUIRE_FALSE(corsika::qgsjetII::CanInteract(corsika::Code::SigmaC0));
} }
SECTION("cross-section type") { SECTION("cross-section type") {
...@@ -95,7 +95,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { ...@@ -95,7 +95,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
setup::Stack stack; setup::Stack stack;
const HEPEnergyType E0 = 100_GeV; const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 = 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}); auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
corsika::Point pos(cs, 0_m, 0_m, 0_m); corsika::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle( auto particle = stack.AddParticle(
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment