From 47a4ae403b8284c7008d73265f09d3b7b9e25234 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 19 May 2021 12:12:01 +0100
Subject: [PATCH] fixed cross section types, read tables instead of calculate
anew
---
corsika/detail/modules/epos/Interaction.inl | 184 +++++++++++++-------
corsika/modules/epos/Interaction.hpp | 25 ++-
corsika/modules/epos/ParticleConversion.hpp | 6 +-
modules/epos/epos.hpp | 13 +-
src/modules/epos/code_generator.py | 2 +-
tests/modules/testEpos.cpp | 56 +++---
6 files changed, 193 insertions(+), 93 deletions(-)
diff --git a/corsika/detail/modules/epos/Interaction.inl b/corsika/detail/modules/epos/Interaction.inl
index 79cc7cd0a..be56c26dc 100644
--- a/corsika/detail/modules/epos/Interaction.inl
+++ b/corsika/detail/modules/epos/Interaction.inl
@@ -232,7 +232,7 @@ namespace corsika::epos {
"Beam={}, "
"BeamA={}, "
"BeamZ={}, "
- "Target={}"
+ "Target={}, "
"TargetA={}, "
"TargetZ={} ",
idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
@@ -241,8 +241,10 @@ namespace corsika::epos {
::epos::hadr25_.idprojin = convertToEposRaw(Code::Proton);
::epos::nucl1_.laproj = iBeamZ; // get_nucleus_Z(idBeam);
::epos::nucl1_.maproj = iBeamA; // get_nucleus_A(idBeam);
+ ::epos::had10_.iclpro = corsika::epos::getEposXSCode(Code::Proton);
} else {
::epos::hadr25_.idprojin = convertToEposRaw(idBeam);
+ ::epos::had10_.iclpro = corsika::epos::getEposXSCode(idBeam);
::epos::nucl1_.laproj = -1;
::epos::nucl1_.maproj = 1;
}
@@ -251,14 +253,17 @@ namespace corsika::epos {
::epos::hadr25_.idtargin = convertToEposRaw(Code::Proton);
::epos::nucl1_.matarg = iTargetA; // get_nucleus_A(idTarget);
::epos::nucl1_.latarg = iTargetZ; // get_nucleus_Z(idTarget);
+ ::epos::had10_.icltar = corsika::epos::getEposXSCode(Code::Proton);
} else if (idTarget == Code::Proton || idTarget == Code::Hydrogen) {
::epos::hadr25_.idtargin = convertToEposRaw(Code::Proton);
::epos::nucl1_.matarg = 1;
::epos::nucl1_.latarg = -1;
+ ::epos::had10_.icltar = corsika::epos::getEposXSCode(Code::Proton);
} else if (idTarget == Code::Neutron) {
::epos::hadr25_.idtargin = convertToEposRaw(Code::Neutron);
::epos::nucl1_.matarg = 1;
::epos::nucl1_.latarg = -1;
+ ::epos::had10_.icltar = corsika::epos::getEposXSCode(Code::Proton);
} else {
throw std::runtime_error("Epos: configure_particles: target outside range!");
}
@@ -267,33 +272,36 @@ namespace corsika::epos {
"Id beam={}, "
"Z beam={}, "
"A beam={}, "
+ "XS beam={}, "
"Id target={}, "
"Z target={}, "
- "A target={}",
+ "A target={}, "
+ "XS target={} ",
::epos::hadr25_.idprojin, ::epos::nucl1_.laproj,
- ::epos::nucl1_.maproj, ::epos::hadr25_.idtargin,
- ::epos::nucl1_.latarg, ::epos::nucl1_.matarg);
+ ::epos::nucl1_.maproj, ::epos::had10_.iclpro,
+ ::epos::hadr25_.idtargin, ::epos::nucl1_.latarg,
+ ::epos::nucl1_.matarg, ::epos::had10_.icltar);
}
inline Interaction::~Interaction() { CORSIKA_LOGGER_DEBUG(logger_, "n={} ", count_); }
- inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
- Interaction::getCrossSection(corsika::Code const BeamId, corsika::Code const TargetId,
- const corsika::HEPEnergyType EnergyCOM) const {
- if (!is_nucleus(BeamId))
- return getCrossSection(BeamId, 1, 1, TargetId, get_nucleus_A(TargetId),
- get_nucleus_Z(TargetId), EnergyCOM);
- else
- throw("nuclear projecile, call getCrossSection with : BeamId, BeamA, BeamZ, ...");
- }
+ // inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
+ // Interaction::getCrossSection(corsika::Code const BeamId, corsika::Code const TargetId,
+ // const corsika::HEPEnergyType EnergyCOM) const {
+ // if (!is_nucleus(BeamId))
+ // return getCrossSection(BeamId, 1, 1, TargetId, get_nucleus_A(TargetId),
+ // get_nucleus_Z(TargetId), EnergyCOM);
+ // else
+ // throw("nuclear projecile, call getCrossSection with : BeamId, BeamA, BeamZ, ...");
+ // }
inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
- Interaction::getCrossSection(corsika::Code const BeamId, int const BeamA,
+ Interaction::calcCrossSectionCoM(corsika::Code const BeamId, int const BeamA,
int const BeamZ, corsika::Code const TargetId,
int const TargetA, int const TargetZ,
const corsika::HEPEnergyType EnergyCOM) const {
CORSIKA_LOGGER_DEBUG(logger_,
- "getCrossSection: input:"
+ "calcCrossSection: input:"
" beamId={}, beamA={}, beamZ={}"
" target={}, targetA={}, targetZ={}"
" Ecm={:4.3f} GeV,",
@@ -304,29 +312,29 @@ namespace corsika::epos {
BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
if (!iBeam)
throw std::runtime_error(
- "getCrossSection: interaction of beam hadron not defined in "
+ "calcCrossSectionCoM: interaction of beam hadron not defined in "
"Epos!");
CORSIKA_LOGGER_TRACE(logger_,
"projectile cross section type={} "
- "(0: cannot interact, 1:baryon, 2:pion, 3:kaon, 4:nucleus)",
+ "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
iBeam);
- // reset beam particle // (1: proton-like, 2: pion-like, 3:kaon-like, 4:nucleus)
+ // reset beam particle // (1: pion-like, 2: proton-like, 3:kaon-like)
if (iBeam == 1)
- initialize_event_CoM(Code::Proton, BeamA, BeamZ, TargetId, TargetA, TargetZ,
- EnergyCOM);
- else if (iBeam == 2)
initialize_event_CoM(Code::PiPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM);
+ else if (iBeam == 2)
+ initialize_event_CoM(Code::Proton, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+ EnergyCOM);
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 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 "
+ "calcCrossSectionCoM: interaction of beam hadron not defined in "
"Epos!");
//::epos::xsigma_();
@@ -349,33 +357,88 @@ namespace corsika::epos {
//::epos::crseaaepos_(sigTot, sigProd, sigCut, sigEla);
CORSIKA_LOGGER_DEBUG(logger_,
- "getCrossSection: output:"
+ "calcCrossSectionCoM: output:"
" sigProd={} mb,"
" sigEla={} mb",
sigProd, sigEla);
return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
+ }
+
+ inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
+ Interaction::readCrossSectionTableLab(Code const BeamId, int const BeamA,
+ int const BeamZ, Code const TargetId,
+ HEPEnergyType const EnergyLab) const {
+ CORSIKA_LOGGER_DEBUG(logger_,
+ "readCrossSectionTableLab: input: "
+ "beamId={}, "
+ "beamA={}, "
+ "beamZ={} "
+ "targetId={}, "
+ "ELab={:4.3f} GeV,",
+ BeamId, BeamA, BeamZ, TargetId, EnergyLab / 1_GeV);
// read cross section from epos internal tables
+ int Abeam;
+ float Ekin = -1;
+
+ if (is_nucleus(BeamId)) {
+ Abeam = BeamA;
+ // kinetic energy per nucleon
+ Ekin = (EnergyLab / Abeam - constants::nucleonMass) / 1_GeV;
+ } else {
+ ::epos::hadr2_.idproj = convertToEposRaw(BeamId);
+ int const iBeam = corsika::epos::getEposXSCode(
+ BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
+ CORSIKA_LOGGER_TRACE(logger_,
+ "projectile cross section type={} "
+ "(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
+ iBeam);
+
+ ::epos::had10_.iclpro = iBeam;
+ Abeam = 1;
+ Ekin = (EnergyLab - get_mass(BeamId)) / 1_GeV;
+ }
+ if (Ekin < 0) {
+ CORSIKA_LOGGER_ERROR(logger_,
+ "Negative kinetic energy!"
+ "Ekin={}",
+ Ekin);
+ throw std::runtime_error("Epos cross section failed! Negative kinetic energy!");
+ }
+
+ int Atarget;
+ if (is_nucleus(TargetId))
+ Atarget = get_nucleus_A(TargetId);
+ else
+ Atarget = 1;
+
+ int iMode = 3; // 0: air, >0 not air
+
+ CORSIKA_LOGGER_DEBUG(logger_,
+ "inside Epos "
+ "beamId={}, beamXS={}",
+ ::epos::hadr2_.idproj, ::epos::had10_.iclpro);
+
+ float sigProdEpos = ::epos::eposcrse_(Ekin, Abeam, Atarget, iMode);
+ float sigElaEpos = ::epos::eposelacrse_(Ekin, Abeam, Atarget, iMode);
- // int Abeam;
- // int iBeamId;
- // if(is_nucleus(BeamId)){
- // Abeam = get_nucleus_A(BeamId);
- // iBeamId = 1; //convertToEposRaw(BeamId);
- // } else {
- // ::epos::hadr2_.idproj = convertToEposRaw(BeamId);
- // Abeam = 1;
- // }
- // int Atarget;
- // if(is_nucleus(TargetId))
- // Atarget = get_nucleus_A(TargetId);
- // else
-
- // float Ekin = (EnergyLab-get_mass(BeamId)) / 1_GeV;
- // float sigProdEpos = ::epos::eposcrse_(Ekin,iBeamId);
- // float sigElaEpos = ::epos::eposelacrse_();
- // return std::make_tuple(sigProdEpos * 1_mb, sigElaEpos * 1_mb);
+ return std::make_tuple(sigProdEpos * 1_mb, sigElaEpos * 1_mb);
+ }
+
+ inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
+ Interaction::getCrossSectionLab(corsika::Code const BeamId, int const BeamA,
+ int const BeamZ, corsika::Code const TargetId,
+ int const TargetA, int const TargetZ,
+ const corsika::HEPEnergyType EnergyLab) const {
+ CORSIKA_LOGGER_DEBUG(logger_,
+ "getCrossSectionLab: input:"
+ " beamId={}, beamA={}, beamZ={}"
+ " target={}, targetA={}, targetZ={}"
+ " ELab={:4.3f} GeV,",
+ BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ,
+ EnergyLab / 1_GeV);
+ return readCrossSectionTableLab(BeamId, BeamA, BeamZ, TargetId, EnergyLab);
}
template <>
@@ -417,13 +480,13 @@ namespace corsika::epos {
// total momentum and energy
HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass;
- MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
- pTotLab += pLab;
- pTotLab += pTarget;
- auto const pTotLabNorm = pTotLab.getNorm();
+ // MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
+ // pTotLab += pLab;
+ // pTotLab += pTarget;
+ // auto const pTotLabNorm = pTotLab.getNorm();
// calculate cm. energy
- const HEPEnergyType ECoM = sqrt(
- (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
+ // const HEPEnergyType ECoM = sqrt(
+ // (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
auto const* currentNode = projectile.getNode();
const auto& mediumComposition =
@@ -431,9 +494,9 @@ namespace corsika::epos {
si::CrossSectionType weightedProdCrossSection = mediumComposition.getWeightedSum(
[=](corsika::Code targetID) -> si::CrossSectionType {
- return std::get<0>(this->getCrossSection(corsikaBeamId, beamA, beamZ,
- targetID, get_nucleus_A(targetID),
- get_nucleus_Z(targetID), ECoM));
+ return std::get<0>(this->getCrossSectionLab(corsikaBeamId, beamA, beamZ,
+ targetID, get_nucleus_A(targetID),
+ get_nucleus_Z(targetID), Elab));
});
CORSIKA_LOGGER_DEBUG(logger_, "InteractionLength: weighted CrossSection (mb): {} ",
@@ -485,11 +548,11 @@ namespace corsika::epos {
HEPEnergyType const projectileMomentumLabPerNucleon = projectileMomentum / beamA;
// define target
- auto const eTargetLab = 0_GeV + constants::nucleonMass;
- HEPEnergyType const Etot = eProjectileLab + eTargetLab;
+ // auto const eTargetLab = 0_GeV + constants::nucleonMass;
+ // HEPEnergyType const Etot = eProjectileLab + eTargetLab;
// invariant mass, i.e. cm. energy
- HEPEnergyType const Ecm =
- sqrt((Etot + projectileMomentum) * (Etot - projectileMomentum));
+ // HEPEnergyType const Ecm =
+ // sqrt((Etot + projectileMomentum) * (Etot - projectileMomentum));
// sample target mass number
auto const* currentNode = projectile.getNode();
@@ -506,9 +569,9 @@ namespace corsika::epos {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- [[maybe_unused]] auto const [sigProd, sigEla] =
- getCrossSection(corsikaBeamId, beamA, beamZ, targetId,
- get_nucleus_A(targetId), get_nucleus_Z(targetId), Ecm);
+ [[maybe_unused]] auto const [sigProd, sigEla] = getCrossSectionLab(
+ corsikaBeamId, beamA, beamZ, targetId, get_nucleus_A(targetId),
+ get_nucleus_Z(targetId), eProjectileLab);
cross_section_of_components[i] = sigProd;
}
@@ -535,7 +598,6 @@ namespace corsika::epos {
::epos::afinal_();
if (epos_listing_) {
- //std::string nam = "EPOSLHC&";
char nam[9] = "EPOSLHC&";
::epos::alistf_(nam);
}
diff --git a/corsika/modules/epos/Interaction.hpp b/corsika/modules/epos/Interaction.hpp
index 0f7c1f6ba..647fbec75 100644
--- a/corsika/modules/epos/Interaction.hpp
+++ b/corsika/modules/epos/Interaction.hpp
@@ -27,13 +27,30 @@ namespace corsika::epos {
//! returns production and elastic cross section for hadrons in epos. Inputs are:
//! CorsikaId of beam particle, CorsikaId of target particle, center-of-mass energy.
- //! Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
- std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
+ //! Allowed targets are: nuclei or single nucleons (p,n,hydrogen). This routine
+ //! calculates the cross sections from scratch. Very slow!
+ std::tuple<CrossSectionType, CrossSectionType> calcCrossSectionCoM(
Code const, int const, int const, Code const, int const, int const,
HEPEnergyType const) const;
- std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
- Code const, Code const, HEPEnergyType const) const;
+ //! returns production and elastic cross section for hadrons in epos by reading
+ //! pre-calculated tables from epos.
+ std::tuple<CrossSectionType, CrossSectionType> readCrossSectionTableLab(
+ Code const, int const, int const, Code const, HEPEnergyType const) const;
+
+ //! returns production and elastic cross section. Allowed configurations are
+ //! hadron-nucleon, hadron-nucleus and nucleus-nucleus. Inputs are particle id's mass
+ //! and charge numbers and total energy in the lab.
+ std::tuple<CrossSectionType, CrossSectionType> getCrossSectionLab(
+ Code const, int const, int const, Code const, int const, int const,
+ HEPEnergyType const) const;
+
+ // std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
+ // Code const, int const, int const, Code const, int const, int const,
+ // HEPEnergyType const) const;
+
+ // std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
+ // Code const, Code const, HEPEnergyType const) const;
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
diff --git a/corsika/modules/epos/ParticleConversion.hpp b/corsika/modules/epos/ParticleConversion.hpp
index 252ebea38..e50c80848 100644
--- a/corsika/modules/epos/ParticleConversion.hpp
+++ b/corsika/modules/epos/ParticleConversion.hpp
@@ -25,10 +25,10 @@ namespace corsika::epos {
*/
enum class EposXSClass : int8_t {
CannotInteract = 0,
- Baryon = 1,
- Pion = 2,
+ Baryon = 2,
+ Pion = 1,
Kaon = 3,
- Nucleus = 4
+ Charm = 4,
};
using EposXSClassIntType = std::underlying_type<EposXSClass>::type;
diff --git a/modules/epos/epos.hpp b/modules/epos/epos.hpp
index d8b77d359..a5ad0c9ff 100644
--- a/modules/epos/epos.hpp
+++ b/modules/epos/epos.hpp
@@ -52,8 +52,8 @@ namespace epos {
void conini_();
void psaini_();
- void idspin_(int&, int&, int&, int&);
- void iclass_(int&, int&);
+ //void idspin_(int&, int&, int&, int&);
+ // void iclass_(int&, int&);
void emsini_(double&, int&, int&);
void paramini_(int&);
void xsigma_();
@@ -489,6 +489,15 @@ namespace epos {
int nody[mxnody];
} nodcy_;
+
+ // integer iclpro,icltar,iclegy
+ // common/had10/iclpro,icltar,iclegy
+ extern struct {
+ int iclpro;
+ int icltar;
+ int iclegy;
+ } had10_;
+
/**
Small helper class to provide a data-directory name in the format eposlhc expects
*/
diff --git a/src/modules/epos/code_generator.py b/src/modules/epos/code_generator.py
index 970495cc0..74d2e05e2 100755
--- a/src/modules/epos/code_generator.py
+++ b/src/modules/epos/code_generator.py
@@ -57,7 +57,7 @@ def set_default_epos_definition(particle_db):
xsType = "CannotInteract"
hadronType = "UndefinedType"
if (pData['isNucleus']):
- xsType = "Nucleus"
+ xsType = "Baryon"
hadronType = "NucleusType"
pData['epos_xsType'] = xsType
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index 68b31673d..75f74acb3 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -119,6 +119,11 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
return sum;
}
+auto sqs2elab(HEPEnergyType const sqs, HEPEnergyType const ma,
+ HEPEnergyType const mb){
+ return (sqs * sqs - ma * ma - mb * mb) / 2. / mb;
+}
+
TEST_CASE("EposInterface", "[processes]") {
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
@@ -148,17 +153,17 @@ TEST_CASE("EposInterface", "[processes]") {
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);
+ // hydrogen target == proton target == neutron target
+ auto const [xs_prod_pp, xs_ela_pp] =
+ model.getCrossSectionLab(Code::Proton, 1, 1, Code::Proton, 1, 1, 100_GeV);
+ auto const [xs_prod_pn, xs_ela_pn] =
+ model.getCrossSectionLab(Code::Proton, 1, 1, Code::Neutron, 1, 0, 100_GeV);
+ auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] =
+ model.getCrossSectionLab(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") {
@@ -166,19 +171,22 @@ TEST_CASE("EposInterface", "[processes]") {
Interaction model;
// eposlhc accepts protons or nuclei with 4<=A<=18 as targets
- // p-p at 7TeV
+ // p-p at 7TeV around 70mb according to LHC
auto const [xs_prod, xs_ela] =
- model.getCrossSection(Code::Proton, Code::Proton, 7_TeV);
+ model.getCrossSectionLab(Code::Proton, 1, 1, Code::Proton, 1, 1,
+ sqs2elab(7_TeV, Proton::mass, Proton::mass));
CHECK(xs_prod / 1_mb == Approx(70.7).margin(2.1));
// pi-n at 7TeV
auto const [xs_prod1, xs_ela1] =
- model.getCrossSection(Code::PiPlus, Code::Neutron, 7_TeV);
+ model.getCrossSectionLab(Code::PiPlus, 0, 0, Code::Neutron, 1, 0,
+ sqs2elab(7_TeV, PiPlus::mass, Neutron::mass));
CHECK(xs_prod1 / 1_mb == Approx(52.7).margin(2.1));
// k-p at 7TeV
auto const [xs_prod2, xs_ela2] =
- model.getCrossSection(Code::KPlus, Code::Proton, 7_TeV);
+ model.getCrossSectionLab(Code::KPlus, 0, 0, Code::Proton, 1, 1,
+ sqs2elab(7_TeV, KPlus::mass, Proton::mass));
CHECK(xs_prod2 / 1_mb == Approx(45.7).margin(2.1));
}
@@ -187,15 +195,20 @@ TEST_CASE("EposInterface", "[processes]") {
Interaction model;
// eposlhc accepts protons or nuclei with 4<=A<=18 as targets
- auto const [xs_prod, xs_ela] =
- model.getCrossSection(Code::Proton, Code::Oxygen, 100_GeV);
- CHECK(xs_prod / 1_mb == Approx(327.7).margin(5.1));
+ auto const [xs_prod, xs_ela] = model.getCrossSectionLab(
+ Code::Proton, 1, 1, Code::Oxygen, Oxygen::nucleus_A, Oxygen::nucleus_Z, 100_GeV);
+ CHECK(xs_prod / 1_mb == Approx(287.0).margin(5.1));
- auto const [xs_prod2, xs_ela2] = model.getCrossSection(
+ auto const [xs_prod2, xs_ela2] = model.getCrossSectionLab(
Code::Nitrogen, Nitrogen::nucleus_A, Nitrogen::nucleus_Z, Code::Oxygen,
- Oxygen::nucleus_A, Oxygen::nucleus_Z, 10_GeV);
+ Oxygen::nucleus_A, Oxygen::nucleus_Z, 400_GeV);
CHECK(xs_prod2 / 1_mb == Approx(1076.7).margin(3.1));
+ // nuclear stack extension, particle "Nucleus"
+ auto const [xs_prod3, xs_ela3] = model.getCrossSectionLab(
+ Code::Nucleus, Nitrogen::nucleus_A, Nitrogen::nucleus_Z, Code::Oxygen,
+ Oxygen::nucleus_A, Oxygen::nucleus_Z, 400_GeV);
+ CHECK(xs_prod2 / xs_prod3 == 1);
}
SECTION("InteractionInterface - low energy") {
@@ -252,8 +265,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(12.8).margin(0.1));
+ CHECK(length / 1_g * 1_cm * 1_cm == Approx(12.8).margin(4.1));
}
-
}
--
GitLab