IAP GITLAB

Skip to content
Snippets Groups Projects
Commit a932b6c6 authored by Felix Riehn's avatar Felix Riehn Committed by ralfulrich
Browse files

fixed nuclear projectile cross section, cleanup test

parent 488e5403
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
corsika/detail/modules/epos/Interaction.inl
+ 37
− 9
View file @ a932b6c6
......@@ -62,6 +62,22 @@ namespace corsika::epos {
}
}
inline bool Interaction::isValidTarget(Code const TargetId) const {
if (is_nucleus(TargetId))
if (TargetId == Code::Nucleus) {
// nuclearExtension for projectiles only
CORSIKA_LOGGER_WARN(logger_,
"Invalid target!"
" Code::Nucleus only allowed for "
"projectiles! "
"This should not happen!");
return false;
} else {
return (get_nucleus_Z(TargetId) < maxTargetMassNumber_ ? true : false);
}
return false;
}
inline void Interaction::initialize_eposlhc_c7() const {
CORSIKA_LOGGER_DEBUG(logger_, "initializing...");
......@@ -81,8 +97,8 @@ namespace corsika::epos {
::epos::cseed_.seedj = 1;
::epos::cseed_.seedc = 1;
::epos::enrgy_.egymin = 6.;
::epos::enrgy_.egymax = 2.e6;
::epos::enrgy_.egymin = minEnergyCoM_ / 1_GeV; // 6.;
::epos::enrgy_.egymax = maxEnergyCoM_ / 1_GeV; // 2.e6;
::epos::lhcparameters_();
......@@ -211,9 +227,13 @@ namespace corsika::epos {
int const iTargetZ) const {
CORSIKA_LOGGER_TRACE(logger_,
"configure_particles: setting "
"Beam={} "
"Target={}",
idBeam, idTarget);
"Beam={}, "
"BeamA={}, "
"BeamZ={}, "
"Target={}"
"TargetA={}, "
"TargetZ={} ",
idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
if (is_nucleus(idBeam)) {
::epos::hadr25_.idprojin = convertToEposRaw(Code::Proton);
......@@ -272,10 +292,11 @@ namespace corsika::epos {
const corsika::HEPEnergyType EnergyCOM) const {
CORSIKA_LOGGER_DEBUG(logger_,
"getCrossSection: input:"
" beam={},"
" target={},"
" beamId={}, beamA={}, beamZ={}"
" target={}, targetA={}, targetZ={}"
" Ecm={:4.3f} GeV,",
BeamId, TargetId, EnergyCOM / 1_GeV);
BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM / 1_GeV);
const int iBeam = corsika::epos::getEposXSCode(
BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
......@@ -284,7 +305,11 @@ namespace corsika::epos {
"getCrossSection: interaction of beam hadron not defined in "
"Epos!");
// reset beam particle // (1: proton-like, 2: pion-like, 3:kaon-like)
CORSIKA_LOGGER_TRACE(logger_,
"projectile cross section type={} "
"(0: cannot interact, 1:baryon, 2:pion, 3:kaon, 4:nucleus)",
iBeam);
// reset beam particle // (1: proton-like, 2: pion-like, 3:kaon-like, 4:nucleus)
if (iBeam == 1)
initialize_event_CoM(Code::Proton, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM);
......@@ -294,6 +319,9 @@ namespace corsika::epos {
else if (iBeam == 3)
initialize_event_CoM(Code::KPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM);
else if (iBeam == 4)
initialize_event_CoM(Code::Nucleus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM);
else
throw std::runtime_error(
"getCrossSection: interaction of beam hadron not defined in "
......
corsika/modules/epos/Interaction.hpp
+ 5
− 7
View file @ a932b6c6
......@@ -48,12 +48,10 @@ namespace corsika::epos {
bool isValidCoMEnergy(HEPEnergyType const ecm) const {
return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
}
//! eposlhc only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
//! eposlhc only accepts nuclei with X<=A<=Y as targets, or protons aka Hydrogen or
//! neutrons (p,n == nucleon)
bool isValidTarget(Code const TargetId) const {
return false;
}
bool isValidTarget(Code const) const;
void initialize_eposlhc_c7() const;
void initialize_event_CoM(Code const, int const, int const, Code const, int const,
int const, HEPEnergyType const) const;
......@@ -66,8 +64,8 @@ namespace corsika::epos {
private:
default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("epos");
std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_epos_Interaction");
HEPEnergyType const minEnergyCoM_ = -10. * 1e9 * electronvolt;
HEPEnergyType const maxEnergyCoM_ = -1.e6 * 1e9 * electronvolt;
HEPEnergyType const minEnergyCoM_ = 6 * 1e9 * electronvolt;
HEPEnergyType const maxEnergyCoM_ = 2.e6 * 1e9 * electronvolt;
int const maxTargetMassNumber_ = 20;
int const minNuclearTargetA_ = 4;
};
......
corsika/modules/epos/ParticleConversion.hpp
+ 1
− 0
View file @ a932b6c6
......@@ -28,6 +28,7 @@ namespace corsika::epos {
Baryon = 1,
Pion = 2,
Kaon = 3,
Nucleus = 4
};
using EposXSClassIntType = std::underlying_type<EposXSClass>::type;
......
src/modules/epos/code_generator.py
+ 1
− 1
View file @ a932b6c6
......@@ -57,7 +57,7 @@ def set_default_epos_definition(particle_db):
xsType = "CannotInteract"
hadronType = "UndefinedType"
if (pData['isNucleus']):
xsType = "Baryon"
xsType = "Nucleus"
hadronType = "NucleusType"
pData['epos_xsType'] = xsType
......
tests/modules/testEpos.cpp
+ 24
− 18
View file @ a932b6c6
......@@ -53,14 +53,14 @@ TEST_CASE("Epos", "[processes]") {
SECTION("canInteractInEpos") {
CHECK(corsika::epos::canInteract(Code::Proton));
CHECK_FALSE(corsika::epos::canInteract(Code::Electron));
CHECK_FALSE(corsika::epos::canInteract(Code::Nucleus));
CHECK_FALSE(corsika::epos::canInteract(Code::Helium));
CHECK(corsika::epos::canInteract(Code::Nucleus));
CHECK(corsika::epos::canInteract(Code::Helium));
}
SECTION("cross-section type") {
CHECK(corsika::epos::getEposXSCode(Code::Electron) == 0);
CHECK(corsika::epos::getEposXSCode(Code::K0Long) == 3);
CHECK(corsika::epos::getEposXSCode(Code::SigmaPlus) == 1);
CHECK(corsika::epos::getEposXSCode(Code::K0Long) == 0);
CHECK(corsika::epos::getEposXSCode(Code::SigmaPlus) == 0);
CHECK(corsika::epos::getEposXSCode(Code::KMinus) == 3);
CHECK(corsika::epos::getEposXSCode(Code::PiMinus) == 2);
CHECK(corsika::epos::getEposXSCode(Code::Proton) == 1);
......@@ -144,23 +144,23 @@ TEST_CASE("EposInterface", "[processes]") {
// eposlhc accepts protons or nuclei with 4<=A<=18 as targets
CHECK_FALSE(model.isValidTarget(Code::Electron));
CHECK(model.isValidTarget(Code::Hydrogen));
CHECK_FALSE(model.isValidTarget(Code::Deuterium));
//CHECK_FALSE(model.isValidTarget(Code::Deuterium));
//CHECK_FALSE(model.isValidTarget(Code::Helium3));
CHECK(model.isValidTarget(Code::Helium));
CHECK_FALSE(model.isValidTarget(Code::Helium3));
CHECK_FALSE(model.isValidTarget(Code::Iron));
CHECK(model.isValidTarget(Code::Oxygen));
// hydrogen target == proton target == neutron target
auto const [xs_prod_pp, xs_ela_pp] =
model.getCrossSection(Code::Proton, Code::Proton, 100_GeV);
auto const [xs_prod_pn, xs_ela_pn] =
model.getCrossSection(Code::Proton, Code::Neutron, 100_GeV);
auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] =
model.getCrossSection(Code::Proton, 1, 1, Code::Hydrogen, 1, 1, 100_GeV);
CHECK(xs_prod_pp == xs_prod_pHydrogen);
CHECK(xs_prod_pp == xs_prod_pn);
CHECK(xs_ela_pp == xs_ela_pHydrogen);
CHECK(xs_ela_pn == xs_ela_pHydrogen);
// auto const [xs_prod_pp, xs_ela_pp] =
// model.getCrossSection(Code::Proton, Code::Proton, 100_GeV);
// auto const [xs_prod_pn, xs_ela_pn] =
// model.getCrossSection(Code::Proton, Code::Neutron, 100_GeV);
// auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] =
// model.getCrossSection(Code::Proton, 1, 1, Code::Hydrogen, 1, 1, 100_GeV);
// CHECK(xs_prod_pp == xs_prod_pHydrogen);
// CHECK(xs_prod_pp == xs_prod_pn);
// CHECK(xs_ela_pp == xs_ela_pHydrogen);
// CHECK(xs_ela_pn == xs_ela_pHydrogen);
}
SECTION("InteractionInterface - hadron cross sections") {
......@@ -191,7 +191,13 @@ TEST_CASE("EposInterface", "[processes]") {
auto const [xs_prod, xs_ela] =
model.getCrossSection(Code::Proton, Code::Oxygen, 100_GeV);
CHECK(xs_prod / 1_mb == Approx(330.7).margin(3.1));
CHECK(xs_prod / 1_mb == Approx(327.7).margin(5.1));
auto const [xs_prod2, xs_ela2] = model.getCrossSection(
Code::Nitrogen, Nitrogen::nucleus_A, Nitrogen::nucleus_Z, Code::Oxygen,
Oxygen::nucleus_A, Oxygen::nucleus_Z, 10_GeV);
CHECK(xs_prod2 / 1_mb == Approx(1076.7).margin(3.1));
}
SECTION("InteractionInterface - low energy") {
......@@ -248,7 +254,7 @@ TEST_CASE("EposInterface", "[processes]") {
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
CHECK(length / 1_g * 1_cm * 1_cm == Approx(75.2).margin(0.1));
CHECK(length / 1_g * 1_cm * 1_cm == Approx(12.8).margin(0.1));
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum