diff --git a/corsika/detail/modules/epos/Interaction.inl b/corsika/detail/modules/epos/Interaction.inl
index dc4b0df99ab194587631f6f0d24dc58181523927..c7724bad3359bc859af3232d24c59d5d7d48bb4d 100644
--- a/corsika/detail/modules/epos/Interaction.inl
+++ b/corsika/detail/modules/epos/Interaction.inl
@@ -52,8 +52,8 @@ namespace corsika::epos {
::epos::aaset_(iarg);
//::epos::atitle_();
- //::epos::prnt1_.ish = 3; // debug level in epos
- //::epos::prnt1_.ifch = 6; // output unit
+ ::epos::prnt1_.ish = 3; // debug level in epos
+ ::epos::files_.ifch = 6; // output unit
//::epos::prnt1_.iecho = 1;
// dummy set seeds for random number generator in epos. need to fool epos checks...
@@ -279,7 +279,7 @@ namespace corsika::epos {
if (targetMassNumber > maxTargetMassNumber_)
throw std::runtime_error("Epos target mass outside range.");
} else {
- if (targetCode != Proton::code || targetCode != Neutron::code)
+ if (targetCode != Proton::code && targetCode != Neutron::code)
throw std::runtime_error("Epos target not possible.");
// proton or neutron target
::epos::hadr25_.idtargin = convertToEposRaw(targetCode);
@@ -318,16 +318,26 @@ namespace corsika::epos {
// secondaries
EposStack es;
CORSIKA_LOG_DEBUG("npart: {}", es.getSize());
+ int i=-1;
for (auto& psec : es) {
- if (psec.hasDecayed()) continue;
+ ++i;
+ if (!psec.isFinal()) continue;
+
+ CORSIKA_LOG_DEBUG("EPOS: i, id, energy, mass, type, state: {} {} {} {} {}", i,
+ ::epos::cptl_.idptl[i], ::epos::cptl_.pptl[3][i],
+ ::epos::cptl_.pptl[4][i], ::epos::cptl_.ityptl[i],
+ ::epos::cptl_.istptl[i]);
+ CORSIKA_LOG_DEBUG("id, energy: {} {}", psec.getPID(), psec.getEnergy());
+ int id = abs(static_cast<int>(psec.getPID()));
+ CORSIKA_LOG_DEBUG("epos id to pdg: {}",
+ ::epos::idtrafo_("nxs", "pdg", id));
+
+
auto momentum = psec.getMomentum(zAxisFrame);
auto const energy = psec.getEnergy();
momentum.rebase(originalCS); // transform back into standard lab frame
- CORSIKA_LOG_DEBUG("id, energy: {} {}", psec.getPID(), psec.getEnergy()/1_GeV);
- int id = abs(static_cast<int>(psec.getPID()));
- CORSIKA_LOG_DEBUG("epos id to pdg:",
- ::epos::idtrafo_("nxs", "pdg", id));
+
auto const pid = corsika::epos::convertFromEpos(psec.getPID());
CORSIKA_LOG_DEBUG(
@@ -341,7 +351,7 @@ namespace corsika::epos {
}
CORSIKA_LOG_DEBUG(
"conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
- "Elab_final="
+ ", Elab_final={}"
", Plab_final={}",
Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
}
diff --git a/corsika/modules/epos/EposStack.hpp b/corsika/modules/epos/EposStack.hpp
index f40f71a616cf5544575f905206ee5a3eed6a7a6a..3be163e629f789bf5cde25d41c858b94ad755909 100644
--- a/corsika/modules/epos/EposStack.hpp
+++ b/corsika/modules/epos/EposStack.hpp
@@ -31,38 +31,38 @@ namespace corsika::epos {
void setId(const unsigned int i, const int v) { ::epos::cptl_.idptl[i] = v; }
void setEnergy(const unsigned int i, const HEPEnergyType v) {
- ::epos::cptl_.pptl[3][i] = v / 1_GeV;
+ ::epos::cptl_.pptl[i][3] = v / 1_GeV;
}
void setMass(const unsigned int i, const HEPMassType v) {
- ::epos::cptl_.pptl[4][i] = v / 1_GeV;
+ ::epos::cptl_.pptl[i][4] = v / 1_GeV;
}
void setMomentum(const unsigned int i, const MomentumVector& v) {
auto tmp = v.getComponents();
- for (int idx = 0; idx < 3; ++idx) ::epos::cptl_.pptl[idx][i] = tmp[idx] / 1_GeV;
+ for (int idx = 0; idx < 3; ++idx) ::epos::cptl_.pptl[i][idx] = tmp[idx] / 1_GeV;
}
void setState(const unsigned int i, const int v) { ::epos::cptl_.istptl[i] = v; }
int getId(const unsigned int i) const { return ::epos::cptl_.idptl[i]; }
int getState(const unsigned int i) const { return ::epos::cptl_.istptl[i]; }
HEPEnergyType getEnergy(const int i) const {
- return ::epos::cptl_.pptl[3][i] * 1_GeV;
+ return ::epos::cptl_.pptl[i][3] * 1_GeV;
}
HEPEnergyType getMass(const unsigned int i) const {
- return ::epos::cptl_.pptl[4][i] * 1_GeV;
+ return ::epos::cptl_.pptl[i][4] * 1_GeV;
}
MomentumVector getMomentum(const unsigned int i) const {
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
- QuantityVector<hepmomentum_d> components = {::epos::cptl_.pptl[0][i] * 1_GeV,
- ::epos::cptl_.pptl[1][i] * 1_GeV,
- ::epos::cptl_.pptl[2][i] * 1_GeV};
+ QuantityVector<hepmomentum_d> components = {::epos::cptl_.pptl[i][0] * 1_GeV,
+ ::epos::cptl_.pptl[i][1] * 1_GeV,
+ ::epos::cptl_.pptl[i][2] * 1_GeV};
return MomentumVector(rootCS, components);
}
MomentumVector getMomentum(const unsigned int i,
const CoordinateSystemPtr& CS) const {
- QuantityVector<hepmomentum_d> components = {::epos::cptl_.pptl[0][i] * 1_GeV,
- ::epos::cptl_.pptl[1][i] * 1_GeV,
- ::epos::cptl_.pptl[2][i] * 1_GeV};
+ QuantityVector<hepmomentum_d> components = {::epos::cptl_.pptl[i][0] * 1_GeV,
+ ::epos::cptl_.pptl[i][1] * 1_GeV,
+ ::epos::cptl_.pptl[i][2] * 1_GeV};
return MomentumVector(CS, components);
}
@@ -73,14 +73,14 @@ namespace corsika::epos {
::epos::cptl_.istptl[i2] = ::epos::cptl_.istptl[i1];
::epos::cptl_.ityptl[i2] = ::epos::cptl_.ityptl[i1];
for (unsigned int i = 0; i < 5; ++i)
- ::epos::cptl_.pptl[i][i2] = ::epos::cptl_.pptl[i][i1];
+ ::epos::cptl_.pptl[i2][i] = ::epos::cptl_.pptl[i1][i];
for (unsigned int i = 0; i < 2; ++i) {
- ::epos::cptl_.tivptl[i][i2] = ::epos::cptl_.tivptl[i][i1];
- ::epos::cptl_.ifrptl[i][i2] = ::epos::cptl_.ifrptl[i][i1];
+ ::epos::cptl_.tivptl[i2][i] = ::epos::cptl_.tivptl[i1][i];
+ ::epos::cptl_.ifrptl[i2][i] = ::epos::cptl_.ifrptl[i1][i];
}
for (unsigned int i = 0; i < 4; ++i) {
- ::epos::cptl_.xorptl[i][i2] = ::epos::cptl_.xorptl[i][i1];
- ::epos::cptl_.ibptl[i][i2] = ::epos::cptl_.ibptl[i][i1];
+ ::epos::cptl_.xorptl[i2][i] = ::epos::cptl_.xorptl[i1][i];
+ ::epos::cptl_.ibptl[i2][i] = ::epos::cptl_.ibptl[i1][i];
}
}
@@ -91,14 +91,14 @@ namespace corsika::epos {
std::swap(::epos::cptl_.istptl[i2], ::epos::cptl_.istptl[i1]);
std::swap(::epos::cptl_.ityptl[i2], ::epos::cptl_.ityptl[i1]);
for (unsigned int i = 0; i < 5; ++i)
- std::swap(::epos::cptl_.pptl[i][i2], ::epos::cptl_.pptl[i][i1]);
+ std::swap(::epos::cptl_.pptl[i2][i], ::epos::cptl_.pptl[i1][i]);
for (unsigned int i = 0; i < 2; ++i) {
- std::swap(::epos::cptl_.tivptl[i][i2], ::epos::cptl_.tivptl[i][i1]);
- std::swap(::epos::cptl_.ifrptl[i][i2], ::epos::cptl_.ifrptl[i][i1]);
+ std::swap(::epos::cptl_.tivptl[i2][i], ::epos::cptl_.tivptl[i1][i]);
+ std::swap(::epos::cptl_.ifrptl[i2][i], ::epos::cptl_.ifrptl[i1][i]);
}
for (unsigned int i = 0; i < 4; ++i) {
- std::swap(::epos::cptl_.xorptl[i][i2], ::epos::cptl_.xorptl[i][i1]);
- std::swap(::epos::cptl_.ibptl[i][i2], ::epos::cptl_.ibptl[i][i1]);
+ std::swap(::epos::cptl_.xorptl[i2][i], ::epos::cptl_.xorptl[i1][i]);
+ std::swap(::epos::cptl_.ibptl[i2][i], ::epos::cptl_.ibptl[i1][i]);
}
}
@@ -140,7 +140,7 @@ namespace corsika::epos {
HEPEnergyType getEnergy() const { return getStackData().getEnergy(getIndex()); }
- bool hasDecayed() const { return getStackData().getState(getIndex()) == 0; }
+ bool isFinal() const { return getStackData().getState(getIndex()) == 0; }
void setMass(const HEPMassType v) { getStackData().setMass(getIndex(), v); }
diff --git a/modules/epos/epos.hpp b/modules/epos/epos.hpp
index 383fae4317f198e6632dc7bbd7ad7a8bbe1e9f37..68ead99755e44de8f052d4f6fc925ba05750aff3 100644
--- a/modules/epos/epos.hpp
+++ b/modules/epos/epos.hpp
@@ -342,6 +342,19 @@ namespace epos {
int ioidch;
} othe1_;
+ // integer ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr
+ // common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr
+ extern struct {
+ int ifop;
+ int ifmt;
+ int ifch;
+ int ifcx;
+ int ifhi;
+ int ifdt;
+ int ifcp;
+ int ifdr;
+ } files_;
+
// character*500 fnch,fnhi,fndt,fnii,fnid,fnie,fnrj,fnmt
// * ,fngrv,fncp,fnnx,fncs,fndr,fnhpf
// common/fname/ fnch, fnhi, fndt, fnii, fnid, fnie, fnrj, fnmt
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index 04e667fa7c1a0fe59499078ade94b0ec14ca597d..752e1f1c708ce2e34149f6fed2e19ee62968f30b 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -122,7 +122,7 @@ TEST_CASE("EposInterface", "[processes]") {
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
- auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
@@ -173,6 +173,16 @@ TEST_CASE("EposInterface", "[processes]") {
Interaction model;
model.doInteraction(view);
+ auto const pSum = sumMomentum(view, cs);
+
+ CHECK(pSum.getComponents(cs).getX() / P0 == Approx(1).margin(0.05));
+ CHECK(pSum.getComponents(cs).getY() / 1_GeV == Approx(0).margin(1e-4));
+ CHECK(pSum.getComponents(cs).getZ() / 1_GeV == Approx(0).margin(1e-4));
+
+ CHECK((pSum - plab).getNorm() / 1_GeV ==
+ Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
+ 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(88.7).margin(0.1));