diff --git a/corsika/detail/modules/epos/Interaction.inl b/corsika/detail/modules/epos/Interaction.inl
index c3053ba3c52083a8e1e3e005f3564c7836ae271a..dc4b0df99ab194587631f6f0d24dc58181523927 100644
--- a/corsika/detail/modules/epos/Interaction.inl
+++ b/corsika/detail/modules/epos/Interaction.inl
@@ -13,8 +13,7 @@
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
-//#include <corsika/modules/sibyll/ParticleConversion.hpp>
-//#include <corsika/modules/sibyll/SibStack.hpp>
+#include <corsika/modules/epos/EposStack.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
@@ -46,7 +45,6 @@ namespace corsika::epos {
}
}
-
inline void Interaction::initialize_eposlhc_c7() {
// corsika7 ini
@@ -57,19 +55,19 @@ namespace corsika::epos {
//::epos::prnt1_.ish = 3; // debug level in epos
//::epos::prnt1_.ifch = 6; // output unit
//::epos::prnt1_.iecho = 1;
-
+
// dummy set seeds for random number generator in epos. need to fool epos checks...
// we will use external generator
- ::epos::cseed_.seedi = 1;
+ ::epos::cseed_.seedi = 1;
::epos::cseed_.seedj = 1;
- ::epos::cseed_.seedc = 1;
+ ::epos::cseed_.seedc = 1;
::epos::enrgy_.egymin = 6.;
::epos::enrgy_.egymax = 2.e6;
-
+
::epos::lhcparameters_();
-
- ::epos::hadr6_.isigma = 0;// do not show cross section
+
+ ::epos::hadr6_.isigma = 0; // do not show cross section
::epos::hadr6_.isetcs = 3; /* !option to obtain pomeron parameters
! 0.....determine parameters but do not use Kfit
! 1.....determine parameters and use Kfit
@@ -83,11 +81,11 @@ namespace corsika::epos {
*/
::epos::cjinti_.ionudi =
1; // !include quasi elastic events but strict calculation of xs
- ::epos::cjinti_.iorsce = 0; // !color exchange turned on(1) or off(0)
- ::epos::cjinti_.iorsdf = 3; // !droplet formation turned on(>0) or off(0)
- ::epos::cjinti_.iorshh = 0; // !other hadron-hadron int. turned on(1) or off(0)
+ ::epos::cjinti_.iorsce = 0; // !color exchange turned on(1) or off(0)
+ ::epos::cjinti_.iorsdf = 3; // !droplet formation turned on(>0) or off(0)
+ ::epos::cjinti_.iorshh = 0; // !other hadron-hadron int. turned on(1) or off(0)
- ::epos::othe1_.istore = 0; // do not produce epos output file
+ ::epos::othe1_.istore = 0; // do not produce epos output file
::epos::nucl6_.infragm = 0; // keep free nucleons in fragmentation
::epos::othe2_.iframe = 12; // lab frame, target at rest
@@ -96,7 +94,7 @@ namespace corsika::epos {
::epos::datadir BASE(data_path_);
strcpy(::epos::fname_.fnnx, BASE.data);
::epos::nfname_.nfnnx = BASE.length;
-
+
::epos::datadir TL(data_path_ + "epos.initl");
strcpy(::epos::fname_.fnii, TL.data);
::epos::nfname_.nfnii = TL.length;
@@ -104,15 +102,15 @@ namespace corsika::epos {
::epos::datadir EV(data_path_ + "epos.iniev");
strcpy(::epos::fname_.fnie, EV.data);
::epos::nfname_.nfnie = EV.length;
-
+
::epos::datadir RJ(data_path_ + "epos.inirj"); // lhcparameters adds ".lhc"
strcpy(::epos::fname_.fnrj, RJ.data);
::epos::nfname_.nfnrj = RJ.length;
-
+
::epos::datadir CS(data_path_ + "epos.inics"); // lhcparameters adds ".lhc"
strcpy(::epos::fname_.fncs, CS.data);
::epos::nfname_.nfncs = CS.length;
-
+
//::epos::fname_.fnid="/home/felix/ngcorsika/corsika/modules/data/EPOS/epos.inidi";
// EPOPAR input ../epos/epos.param !initialization input file for epos
// EPOPAR fname inics ../epos/epos.inics !initialization input file for epos
@@ -134,9 +132,8 @@ namespace corsika::epos {
::epos::nucl1_.latarg = 1;
::epos::hadr1_.pnll = 200.;
::epos::lept1_.engy = -1.;
-
- ::epos::ainit_();
+ ::epos::ainit_();
}
inline Interaction::~Interaction() {
@@ -160,154 +157,195 @@ namespace corsika::epos {
auto const projectile = view.getProjectile();
auto const corsikaBeamId = projectile.getPID();
- if (corsika::is_nucleus(corsikaBeamId)) {
- // nuclei handled by different process, this should not happen
- throw std::runtime_error("Nuclear projectile are not handled by EPOSLHC!");
- }
-
- // position and time of interaction, not used in Epos
- Point const pOrig = projectile.getPosition();
- TimeType const tOrig = projectile.getTime();
-
- // define projectile
- HEPEnergyType const eProjectileLab = projectile.getEnergy();
- auto const pProjectileLab = projectile.getMomentum();
- CoordinateSystemPtr const& originalCS = pProjectileLab.getCoordinateSystem();
-
- CORSIKA_LOG_DEBUG(
- "ProcessEPOS: "
- "DoInteraction: {} interaction ",
- corsikaBeamId);
-
- // define target
- // for Sibyll is always a single nucleon
- // FOR NOW: target is always at rest
- const auto eTargetLab = 0_GeV + constants::nucleonMass;
- const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
- const FourVector PtargLab(eTargetLab, pTargetLab);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: ebeam lab: {} GeV"
- "Interaction: pbeam lab: {} GeV",
- eProjectileLab / 1_GeV, pProjectileLab.getComponents());
- CORSIKA_LOG_DEBUG(
- "Interaction: etarget lab: {} GeV "
- "Interaction: ptarget lab: {} GeV",
- eTargetLab / 1_GeV, pTargetLab.getComponents() / 1_GeV);
-
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
-
- // define target kinematics in lab frame
- // define boost to and from CoM frame
- // CoM frame definition in Sibyll projectile: +z
- COMBoost const boost(PprojLab, constants::nucleonMass);
- auto const& csPrime = boost.getRotatedCS();
-
- // just for show:
- // boost projecticle
- [[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
-
- // boost target
- [[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: ebeam CoM: {} GeV "
- "Interaction: pbeam CoM: {} GeV ",
- PprojCoM.getTimeLikeComponent() / 1_GeV,
- PprojCoM.getSpaceLikeComponents().getComponents(csPrime) / 1_GeV);
- CORSIKA_LOG_DEBUG(
- "Interaction: etarget CoM: {} GeV "
- "Interaction: ptarget CoM: {} GeV ",
- PtargCoM.getTimeLikeComponent() / 1_GeV,
- PtargCoM.getSpaceLikeComponents().getComponents(csPrime) / 1_GeV);
-
- CORSIKA_LOG_DEBUG("Interaction: position of interaction: {} ",
- pOrig.getCoordinates());
- CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig);
-
- HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = projectile.getMomentum();
- // invariant mass, i.e. cm. energy
- HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.getSquaredNorm());
-
- // sample target mass number
- auto const* currentNode = projectile.getNode();
- auto const& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
- // get cross sections for target materials
- /*
- Here we read the cross section from the interaction model again,
- should be passed from getInteractionLength if possible
- */
- //#warning reading interaction cross section again, should not be necessary
- auto const& compVec = mediumComposition.getComponents();
- std::vector<CrossSectionType> cross_section_of_components(compVec.size());
-
- for (size_t i = 0; i < compVec.size(); ++i) {
- auto const targetId = compVec[i];
- [[maybe_unused]] auto const [sigProd, sigEla] =
- getCrossSection(corsikaBeamId, targetId, Ecm);
- cross_section_of_components[i] = sigProd;
- }
+ if (corsika::epos::canInteract(corsikaBeamId)) {
+
+ // position and time of interaction, not used in Epos
+ Point const pOrig = projectile.getPosition();
+ TimeType const tOrig = projectile.getTime();
+
+ // define projectile
+ HEPEnergyType const eProjectileLab = projectile.getEnergy();
+ auto const pProjectileLab = projectile.getMomentum();
+ auto const projectileMomentum = pProjectileLab.getNorm();
+ CoordinateSystemPtr const& originalCS = pProjectileLab.getCoordinateSystem();
+
+ // epos frame with z along the projectile direction
+ CoordinateSystemPtr const zAxisFrame = make_rotationToZ(originalCS, pProjectileLab);
+
+ int beamA = 1;
+ if (is_nucleus(corsikaBeamId)) beamA = projectile.getNuclearA();
+
+ HEPEnergyType const projectileMomentumLabPerNucleon = projectileMomentum / beamA;
+
+ CORSIKA_LOG_DEBUG(
+ "ProcessEPOS: "
+ "DoInteraction: {} interaction ",
+ corsikaBeamId);
+
+ // define target
+ // for Sibyll is always a single nucleon
+ // FOR NOW: target is always at rest
+ const auto eTargetLab = 0_GeV + constants::nucleonMass;
+ const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
+ const FourVector PtargLab(eTargetLab, pTargetLab);
+
+ CORSIKA_LOG_DEBUG(
+ "Interaction: ebeam lab: {} GeV"
+ "Interaction: pbeam lab: {} GeV",
+ eProjectileLab / 1_GeV, pProjectileLab.getComponents());
+ CORSIKA_LOG_DEBUG(
+ "Interaction: etarget lab: {} GeV "
+ "Interaction: ptarget lab: {} GeV",
+ eTargetLab / 1_GeV, pTargetLab.getComponents() / 1_GeV);
+
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Sibyll projectile: +z
+ COMBoost const boost(PprojLab, constants::nucleonMass);
+ auto const& csPrime = boost.getRotatedCS();
+
+ // just for show:
+ // boost projecticle
+ [[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
+
+ // boost target
+ [[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
+
+ CORSIKA_LOG_DEBUG(
+ "Interaction: ebeam CoM: {} GeV "
+ "Interaction: pbeam CoM: {} GeV ",
+ PprojCoM.getTimeLikeComponent() / 1_GeV,
+ PprojCoM.getSpaceLikeComponents().getComponents(csPrime) / 1_GeV);
+ CORSIKA_LOG_DEBUG(
+ "Interaction: etarget CoM: {} GeV "
+ "Interaction: ptarget CoM: {} GeV ",
+ PtargCoM.getTimeLikeComponent() / 1_GeV,
+ PtargCoM.getSpaceLikeComponents().getComponents(csPrime) / 1_GeV);
+
+ CORSIKA_LOG_DEBUG("Interaction: position of interaction: {} ",
+ pOrig.getCoordinates());
+ CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig);
+
+ HEPEnergyType Etot = eProjectileLab + eTargetLab;
+ MomentumVector Ptot = projectile.getMomentum();
+ // invariant mass, i.e. cm. energy
+ HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.getSquaredNorm());
+
+ // sample target mass number
+ auto const* currentNode = projectile.getNode();
+ auto const& mediumComposition =
+ currentNode->getModelProperties().getNuclearComposition();
+ // get cross sections for target materials
+ /*
+ Here we read the cross section from the interaction model again,
+ should be passed from getInteractionLength if possible
+ */
+ //#warning reading interaction cross section again, should not be necessary
+ auto const& compVec = mediumComposition.getComponents();
+ std::vector<CrossSectionType> cross_section_of_components(compVec.size());
+
+ for (size_t i = 0; i < compVec.size(); ++i) {
+ auto const targetId = compVec[i];
+ [[maybe_unused]] auto const [sigProd, sigEla] =
+ getCrossSection(corsikaBeamId, targetId, Ecm);
+ cross_section_of_components[i] = sigProd;
+ }
- const auto targetCode =
- mediumComposition.sampleTarget(cross_section_of_components, RNG_);
- CORSIKA_LOG_DEBUG("Interaction: target selected: {} ", targetCode);
-
- // from corsika7 interface
- // NEXLNK-part
-
- //projectile
- // if(is_nucleus(corsikaBeamId)){
- // ::epos::hadr25_.idprojin =1120;
- // ::epos::nucl1_.laproj = projectile.get_nucleus_Z(); // Z
- // ::epos::nucl1_.maproj = projectile.get_nucleus_A();; // A
- // } else {
- ::epos::hadr25_.idprojin = convertToEposRaw(corsikaBeamId); //1120 ; // id "NEXUS code"
- ::epos::nucl1_.laproj = -1; // Z (-1 for hadron)
- ::epos::nucl1_.maproj = 1; // A
- //}
-
- // target
- int targetMassNumber = 1; // proton
- if (is_nucleus(targetCode)) { // nucleus
- targetMassNumber = get_nucleus_A(targetCode);
- if (targetMassNumber > maxTargetMassNumber_)
- throw std::runtime_error("Epos target mass outside range.");
- } else {
- if (targetCode != Proton::code || targetCode != Neutron::code)
- throw std::runtime_error("Epos target not possible.");
- // proton or neutron target
- ::epos::hadr25_.idtargin = convertToEposRaw(targetCode);
- if (targetCode == Proton::code)
- ::epos::nucl1_.latarg = 1; // Z
- else
- ::epos::nucl1_.latarg = -1; // Z (-1 with id 1220 for neutron)
- ::epos::nucl1_.matarg = 1; // A
+ const auto targetCode =
+ mediumComposition.sampleTarget(cross_section_of_components, RNG_);
+ CORSIKA_LOG_DEBUG("Interaction: target selected: {} ", targetCode);
+
+ // from corsika7 interface
+ // NEXLNK-part
+
+ // projectile
+ // if(is_nucleus(corsikaBeamId)){
+ // ::epos::hadr25_.idprojin =1120;
+ // ::epos::nucl1_.laproj = projectile.get_nucleus_Z(); // Z
+ // ::epos::nucl1_.maproj = projectile.get_nucleus_A();; // A
+ // } else {
+ ::epos::hadr25_.idprojin =
+ convertToEposRaw(corsikaBeamId); // 1120 ; // id "NEXUS code"
+ ::epos::nucl1_.laproj = -1; // Z (-1 for hadron)
+ ::epos::nucl1_.maproj = 1; // A
+ //}
+
+ // target
+ int targetMassNumber = 1; // proton
+ if (is_nucleus(targetCode)) { // nucleus
+ targetMassNumber = get_nucleus_A(targetCode);
+ if (targetMassNumber > maxTargetMassNumber_)
+ throw std::runtime_error("Epos target mass outside range.");
+ } else {
+ if (targetCode != Proton::code || targetCode != Neutron::code)
+ throw std::runtime_error("Epos target not possible.");
+ // proton or neutron target
+ ::epos::hadr25_.idtargin = convertToEposRaw(targetCode);
+ if (targetCode == Proton::code)
+ ::epos::nucl1_.latarg = 1; // Z
+ else
+ ::epos::nucl1_.latarg = -1; // Z (-1 with id 1220 for neutron)
+ ::epos::nucl1_.matarg = 1; // A
+ }
+ CORSIKA_LOG_DEBUG("Interaction: target epos code/A: {}", targetMassNumber);
+
+ // hadron-nucleon momentum
+ ::epos::hadr1_.pnll = float(projectileMomentumLabPerNucleon / 1_GeV); // float(200);
+
+ // C SET ENGY NEGATIVE TO FORCE CALCULATION IN LAB FRAME
+ ::epos::lept1_.engy = -1.;
+ ::epos::enrgy_.ecms = -1.;
+ ::epos::enrgy_.elab = -1.;
+ ::epos::enrgy_.ekin = -1.;
+
+ // C INTIALIZE ENERGY AND PARTICLE DEPENDENT PORTION OF EPOS/NEXUS
+ // C AT THE FIRST CALL: READ ALSO DATA SETS
+ ::epos::ainit_();
+
+ // create event
+ int iarg = 1;
+ ::epos::aepos_(iarg);
+
+ ::epos::afinal_();
+
+ // NSTORE-part
+
+ MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV;
+
+ // secondaries
+ EposStack es;
+ CORSIKA_LOG_DEBUG("npart: {}", es.getSize());
+ for (auto& psec : es) {
+ if (psec.hasDecayed()) continue;
+ 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(
+ "secondary fragment> id= {}"
+ " p= {}",
+ pid, momentum.getComponents() / 1_GeV);
+ auto pnew =
+ view.addSecondary(std::make_tuple(pid, energy, momentum, pOrig, tOrig));
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
+ }
+ CORSIKA_LOG_DEBUG(
+ "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
+ "Elab_final="
+ ", Plab_final={}",
+ Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
}
- CORSIKA_LOG_DEBUG("Interaction: target epos code/A: {}", targetMassNumber);
- // hadron-nucleon momentum
- ::epos::hadr1_.pnll = float(200);
-
- // C SET ENGY NEGATIVE TO FORCE CALCULATION IN LAB FRAME
- ::epos::lept1_.engy = -1.;
- ::epos::enrgy_.ecms = -1.;
- ::epos::enrgy_.elab = -1.;
- ::epos::enrgy_.ekin = -1.;
-
- //C INTIALIZE ENERGY AND PARTICLE DEPENDENT PORTION OF EPOS/NEXUS
- //C AT THE FIRST CALL: READ ALSO DATA SETS
- ::epos::ainit_();
-
- // create event
- int iarg = 1;
- ::epos::aepos_(iarg);
-
- ::epos::afinal_();
- // NSTORE-part
- std::cout << "npart: " << ::epos::cptl_.nptl << std::endl;
-
/*
maproj=maprojxs
laproj=laprojxs
@@ -365,5 +403,4 @@ c xsbmaxim=dble(bmaxim)
*/
}
-
} // namespace corsika::epos
diff --git a/corsika/modules/epos/EposStack.hpp b/corsika/modules/epos/EposStack.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f40f71a616cf5544575f905206ee5a3eed6a7a6a
--- /dev/null
+++ b/corsika/modules/epos/EposStack.hpp
@@ -0,0 +1,174 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+#include <corsika/modules/epos/ParticleConversion.hpp>
+
+#include <epos.hpp>
+
+namespace corsika::epos {
+
+ typedef corsika::Vector<hepmomentum_d> MomentumVector;
+
+ class EposStackData {
+
+ public:
+ void dump() const {}
+
+ void clear() { ::epos::cptl_.nptl = 0; }
+ unsigned int getSize() const { return ::epos::cptl_.nptl; }
+ unsigned int getCapacity() const { return ::epos::mxptl; }
+
+ 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;
+ }
+ void setMass(const unsigned int i, const HEPMassType v) {
+ ::epos::cptl_.pptl[4][i] = 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;
+ }
+ 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;
+ }
+ HEPEnergyType getMass(const unsigned int i) const {
+ return ::epos::cptl_.pptl[4][i] * 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};
+ 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};
+ return MomentumVector(CS, components);
+ }
+
+ void copy(const unsigned int i1, const unsigned int i2) {
+ ::epos::cptl_.idptl[i2] = ::epos::cptl_.idptl[i1];
+ ::epos::cptl_.iorptl[i2] = ::epos::cptl_.iorptl[i1];
+ ::epos::cptl_.jorptl[i2] = ::epos::cptl_.jorptl[i1];
+ ::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];
+ 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];
+ }
+ 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];
+ }
+ }
+
+ void swap(const unsigned int i1, const unsigned int i2) {
+ std::swap(::epos::cptl_.idptl[i2], ::epos::cptl_.idptl[i1]);
+ std::swap(::epos::cptl_.iorptl[i2], ::epos::cptl_.iorptl[i1]);
+ std::swap(::epos::cptl_.jorptl[i2], ::epos::cptl_.jorptl[i1]);
+ 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]);
+ 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]);
+ }
+ 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]);
+ }
+ }
+
+ void incrementSize() { ::epos::cptl_.nptl++; }
+ void decrementSize() {
+ if (::epos::cptl_.nptl > 0) { ::epos::cptl_.nptl--; }
+ }
+ };
+
+ template <typename StackIteratorInterface>
+ class ParticleInterface : public corsika::ParticleBase<StackIteratorInterface> {
+
+ using corsika::ParticleBase<StackIteratorInterface>::getStackData;
+ using corsika::ParticleBase<StackIteratorInterface>::getIndex;
+
+ public:
+ void setParticleData(const int vID, // corsika::epos::EposCode vID,
+ const HEPEnergyType vE, const MomentumVector& vP,
+ const HEPMassType vM) {
+ setPID(vID);
+ setEnergy(vE);
+ setMomentum(vP);
+ setMass(vM);
+ setState(0);
+ }
+
+ void setParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
+ const int vID, // corsika::epos::EposCode vID,
+ const HEPEnergyType vE, const MomentumVector& vP,
+ const HEPMassType vM) {
+ setPID(vID);
+ setEnergy(vE);
+ setMomentum(vP);
+ setMass(vM);
+ setState(0);
+ }
+
+ void setEnergy(const HEPEnergyType v) { getStackData().setEnergy(getIndex(), v); }
+
+ HEPEnergyType getEnergy() const { return getStackData().getEnergy(getIndex()); }
+
+ bool hasDecayed() const { return getStackData().getState(getIndex()) == 0; }
+
+ void setMass(const HEPMassType v) { getStackData().setMass(getIndex(), v); }
+
+ HEPEnergyType getMass() const { return getStackData().getMass(getIndex()); }
+
+ void setPID(const int v) { getStackData().setId(getIndex(), v); }
+
+ corsika::epos::EposCode getPID() const {
+ return static_cast<corsika::epos::EposCode>(getStackData().getId(getIndex()));
+ }
+
+ void setState(const int v) { getStackData().setState(getIndex(), v); }
+
+ corsika::epos::EposCode getState() const {
+ return static_cast<corsika::epos::EposCode>(getStackData().getState(getIndex()));
+ }
+
+ MomentumVector getMomentum() const { return getStackData().getMomentum(getIndex()); }
+
+ MomentumVector getMomentum(const CoordinateSystemPtr& CS) const {
+ return getStackData().getMomentum(getIndex(), CS);
+ }
+
+ void setMomentum(const MomentumVector& v) {
+ getStackData().setMomentum(getIndex(), v);
+ }
+ };
+
+ typedef corsika::Stack<EposStackData, ParticleInterface> EposStack;
+
+} // end namespace corsika::epos
diff --git a/modules/epos/epos.hpp b/modules/epos/epos.hpp
index a454c4a6e954f5f9a2732dd0bd35402c2eba9e41..383fae4317f198e6632dc7bbd7ad7a8bbe1e9f37 100644
--- a/modules/epos/epos.hpp
+++ b/modules/epos/epos.hpp
@@ -76,7 +76,7 @@ namespace epos {
// get mass for id
void idmass_(int&, double&);
// convert id from one format to another
- //int idtrafo_(char[3], char[3], int&);
+ int idtrafo_(char[3], char[3], int&);
// conex routine
// void cxidmass_(int&, int&);
@@ -126,7 +126,6 @@ namespace epos {
int modelxs;
} xsappli_;
-
extern struct {
int nevent;
int nfull;
@@ -139,15 +138,15 @@ namespace epos {
int iframexs;
} xsevent_;
- // common/metr1/iospec,iocova,iopair,iozero,ioflac,iomom
- extern struct {
- int iospec;
- int iocova;
- int iopair;
- int iozero;
- int ioflac;
- int iomom;
- } metr1_;
+ // common/metr1/iospec,iocova,iopair,iozero,ioflac,iomom
+ extern struct {
+ int iospec;
+ int iocova;
+ int iopair;
+ int iozero;
+ int ioflac;
+ int iomom;
+ } metr1_;
extern struct {
int ifrade;
@@ -332,96 +331,105 @@ namespace epos {
// integer istore,istmax,irescl,ntrymx,nclean,iopdg,ioidch
// common/othe1/istore,istmax,gaumx,irescl,ntrymx,nclean,iopdg,ioidch
- extern struct {
- int istore;
- int istmax;
- int gaumx;
- int irescl;
- int ntrymx;
- int nclean;
- int iopdg;
- int ioidch;
- } othe1_;
-
- // 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
- // * ,fngrv,fncp,fnnx,fncs,fndr,fnhpf
- extern struct {
- char fnch[500];
- char fnhi[500];
- char fndt[500];
- char fnii[500];
- char fnid[500];
- char fnie[500];
- char fnrj[500];
- char fnmt[500];
- char fngrv[500];
- char fncp[500];
- char fnnx[500];
- char fncs[500];
- char fndr[500];
- char fnhpf[500];
-
- } fname_;
-
- // integer nfnch,nfnhi,nfndt,nfnii,nfnid,nfnie,nfnrj,nfnmt
- // *,nfngrv,nfncp,nfnnx,nfncs,nfndr,nfnhpf
- // common/nfname/nfnch,nfnhi,nfndt,nfnii,nfnid,nfnie,nfnrj,nfnmt
- // *,nfngrv,nfncp,nfnnx,nfncs,nfndr,nfnhpf
- extern struct {
- int nfnch;
- int nfnhi;
- int nfndt;
- int nfnii;
- int nfnid;
- int nfnie;
- int nfnrj;
- int nfnmt;
- int nfngrv;
- int nfncp;
- int nfnnx;
- int nfncs;
- int nfndr;
- int nfnhpf;
- } nfname_;
-
- // integer iprmpt,ish,ishsub,irandm,irewch,iecho,modsho,idensi
- // common/prnt1/iprmpt,ish,ishsub,irandm,irewch,iecho,modsho,idensi
- extern struct {
- int iprmpt;
- int ish;
- int ishsub;
- int irandm;
- int irewch;
- int iecho;
- int modsho;
- int idensi;
- } prnt1_;
- unsigned int constexpr mmry = 1;
- unsigned int constexpr mxptl = 200000 / mmry;
- // real pptl,tivptl,xorptl
- // integer nptl,iorptl,idptl,istptl,ifrptl,jorptl,ibptl,ityptl
- // common/cptl/nptl,pptl(5,mxptl),iorptl(mxptl),idptl(mxptl)
- extern struct {
- int nptl;
- float pptl[mxptl][5];
- int iorptl[mxptl];
- int idptl[mxptl];
- } cptl_;
-
- /**
- Small helper class to provide a data-directory name in the format eposlhc expects
- */
- class datadir {
- private:
- datadir operator=(const std::string& dir);
- datadir operator=(const datadir&);
-
- public:
- datadir(const std::string& dir);
- char data[500];
- int length;
+ extern struct {
+ int istore;
+ int istmax;
+ int gaumx;
+ int irescl;
+ int ntrymx;
+ int nclean;
+ int iopdg;
+ int ioidch;
+ } othe1_;
+
+ // 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
+ // * ,fngrv,fncp,fnnx,fncs,fndr,fnhpf
+ extern struct {
+ char fnch[500];
+ char fnhi[500];
+ char fndt[500];
+ char fnii[500];
+ char fnid[500];
+ char fnie[500];
+ char fnrj[500];
+ char fnmt[500];
+ char fngrv[500];
+ char fncp[500];
+ char fnnx[500];
+ char fncs[500];
+ char fndr[500];
+ char fnhpf[500];
+
+ } fname_;
+
+ // integer nfnch,nfnhi,nfndt,nfnii,nfnid,nfnie,nfnrj,nfnmt
+ // *,nfngrv,nfncp,nfnnx,nfncs,nfndr,nfnhpf
+ // common/nfname/nfnch,nfnhi,nfndt,nfnii,nfnid,nfnie,nfnrj,nfnmt
+ // *,nfngrv,nfncp,nfnnx,nfncs,nfndr,nfnhpf
+ extern struct {
+ int nfnch;
+ int nfnhi;
+ int nfndt;
+ int nfnii;
+ int nfnid;
+ int nfnie;
+ int nfnrj;
+ int nfnmt;
+ int nfngrv;
+ int nfncp;
+ int nfnnx;
+ int nfncs;
+ int nfndr;
+ int nfnhpf;
+ } nfname_;
+
+ // integer iprmpt,ish,ishsub,irandm,irewch,iecho,modsho,idensi
+ // common/prnt1/iprmpt,ish,ishsub,irandm,irewch,iecho,modsho,idensi
+ extern struct {
+ int iprmpt;
+ int ish;
+ int ishsub;
+ int irandm;
+ int irewch;
+ int iecho;
+ int modsho;
+ int idensi;
+ } prnt1_;
+ unsigned int constexpr mmry = 1;
+ unsigned int constexpr mxptl = 200000 / mmry;
+ // real pptl,tivptl,xorptl
+ // integer nptl,iorptl,idptl,istptl,ifrptl,jorptl,ibptl,ityptl
+ // common/cptl/nptl,pptl(5,mxptl),iorptl(mxptl),idptl(mxptl)
+ // *,istptl(mxptl),tivptl(2,mxptl),ifrptl(2,mxptl),jorptl(mxptl)
+ // *,xorptl(4,mxptl),ibptl(4,mxptl),ityptl(mxptl)
+ extern struct {
+ int nptl;
+ float pptl[mxptl][5];
+ int iorptl[mxptl];
+ int idptl[mxptl];
+ int istptl[mxptl];
+ float tivptl[mxptl][2];
+ int ifrptl[mxptl][2];
+ int jorptl[mxptl];
+ float xorptl[mxptl][4];
+ int ibptl[mxptl][4];
+ int ityptl[mxptl];
+ } cptl_;
+
+ /**
+ Small helper class to provide a data-directory name in the format eposlhc expects
+ */
+ class datadir {
+ private:
+ datadir operator=(const std::string& dir);
+ datadir operator=(const datadir&);
+
+ public:
+ datadir(const std::string& dir);
+ char data[500];
+ int length;
};
}
} // namespace epos