From ec68cdffb4ef1627b2e93d8d72a001f392135bd0 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sun, 31 Oct 2021 08:14:17 +0100
Subject: [PATCH] also adapted EPOS
---
corsika/detail/modules/epos/Interaction.inl | 404 +++++++-----------
.../modules/qgsjetII/InteractionModel.inl | 22 +-
corsika/modules/epos/Interaction.hpp | 57 ++-
corsika/modules/qgsjetII/InteractionModel.hpp | 5 +-
tests/modules/testEpos.cpp | 147 ++++---
tests/modules/testQGSJetII.cpp | 5 +-
6 files changed, 305 insertions(+), 335 deletions(-)
diff --git a/corsika/detail/modules/epos/Interaction.inl b/corsika/detail/modules/epos/Interaction.inl
index 1e69cff4e..b29e8c904 100644
--- a/corsika/detail/modules/epos/Interaction.inl
+++ b/corsika/detail/modules/epos/Interaction.inl
@@ -25,13 +25,10 @@
#include <string>
#include <tuple>
-using namespace corsika;
-using SetupParticle = setup::Stack::stack_iterator_type;
-
namespace corsika::epos {
- inline Interaction::Interaction(std::string const& dataPath,
- bool const epos_printout_on)
+ inline InteractionModel::InteractionModel(std::string const& dataPath,
+ bool const epos_printout_on)
: data_path_(dataPath)
, epos_listing_(epos_printout_on) {
if (dataPath == "") {
@@ -46,7 +43,7 @@ namespace corsika::epos {
setParticlesStable();
}
- inline void Interaction::setParticlesStable() const {
+ inline void InteractionModel::setParticlesStable() const {
CORSIKA_LOGGER_DEBUG(logger_,
"set all particles known to CORSIKA stable inside EPOS..");
for (auto& p : get_all_particles()) {
@@ -59,11 +56,22 @@ namespace corsika::epos {
}
}
- inline bool Interaction::isValidTarget(Code const TargetId) const {
- return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxTargetMassNumber_);
+ inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+ HEPEnergyType const sqrtS) const {
+ //! eposlhc only accepts nuclei with X<=A<=Y as targets, or protons aka Hydrogen or
+ //! neutrons (p,n == nucleon)
+ if (!is_nucleus(targetId) && targetId != Code::Neutron && targetId != Code::Proton) {
+ return false;
+ }
+ if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= maxTargetMassNumber_)) {
+ return false;
+ }
+ if ((minEnergyCoM_ > sqrtS) || (sqrtS > maxEnergyCoM_)) { return false; }
+ if (!epos::canInteract(projectileId)) { return false; }
+ return true;
}
- inline void Interaction::initialize() const {
+ inline void InteractionModel::initialize() const {
CORSIKA_LOGGER_DEBUG(logger_, "initializing...");
@@ -137,10 +145,10 @@ namespace corsika::epos {
Argon::nucleus_A, Argon::nucleus_Z, 100_GeV);
}
- inline void Interaction::initializeEventCoM(Code const idBeam, int const iBeamA,
- int const iBeamZ, Code const idTarget,
- int const iTargetA, int const iTargetZ,
- HEPEnergyType const Ecm) const {
+ inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
+ int const iBeamZ, Code const idTarget,
+ int const iTargetA, int const iTargetZ,
+ HEPEnergyType const Ecm) const {
CORSIKA_LOGGER_TRACE(logger_,
"initialize event in CoM frame!"
" Ecm={}",
@@ -163,10 +171,10 @@ namespace corsika::epos {
::epos::ainit_();
}
- inline void Interaction::initializeEventLab(Code const idBeam, int const iBeamA,
- int const iBeamZ, Code const idTarget,
- int const iTargetA, int const iTargetZ,
- HEPEnergyType const Plab) const {
+ inline void InteractionModel::initializeEventLab(Code const idBeam, int const iBeamA,
+ int const iBeamZ, Code const idTarget,
+ int const iTargetA, int const iTargetZ,
+ HEPEnergyType const Plab) const {
CORSIKA_LOGGER_TRACE(logger_,
"initialize event in lab. frame!"
" Plab per nuc={} GeV",
@@ -191,10 +199,10 @@ namespace corsika::epos {
::epos::ainit_();
}
- inline void Interaction::configureParticles(Code const idBeam, int const iBeamA,
- int const iBeamZ, Code const idTarget,
- int const iTargetA,
- int const iTargetZ) const {
+ inline void InteractionModel::configureParticles(Code const idBeam, int const iBeamA,
+ int const iBeamZ, Code const idTarget,
+ int const iTargetA,
+ int const iTargetZ) const {
CORSIKA_LOGGER_TRACE(logger_,
"setting "
"Beam={}, "
@@ -227,9 +235,8 @@ namespace corsika::epos {
::epos::hadr25_.idtargin = convertToEposRaw(Code::Neutron);
::epos::nucl1_.matarg = 1;
::epos::nucl1_.latarg = -1;
- } else {
- throw std::runtime_error("Epos: target outside range!");
}
+
CORSIKA_LOGGER_TRACE(logger_,
"inside EPOS: "
"Id beam={}, "
@@ -246,11 +253,15 @@ namespace corsika::epos {
::epos::nucl1_.matarg, ::epos::had10_.icltar);
}
- inline Interaction::~Interaction() { CORSIKA_LOGGER_DEBUG(logger_, "n={} ", count_); }
+ inline InteractionModel::~InteractionModel() {
+ CORSIKA_LOGGER_DEBUG(logger_, "n={} ", count_);
+ }
- inline std::tuple<CrossSectionType, CrossSectionType> Interaction::calcCrossSectionCoM(
- Code const BeamId, int const BeamA, int const BeamZ, Code const TargetId,
- int const TargetA, int const TargetZ, const HEPEnergyType EnergyCOM) const {
+ inline std::tuple<CrossSectionType, CrossSectionType>
+ InteractionModel::calcCrossSectionCoM(Code const BeamId, int const BeamA,
+ int const BeamZ, Code const TargetId,
+ int const TargetA, int const TargetZ,
+ const HEPEnergyType EnergyCOM) const {
CORSIKA_LOGGER_DEBUG(logger_,
"calcCrossSection: input:"
" beamId={}, beamA={}, beamZ={}"
@@ -259,12 +270,8 @@ namespace corsika::epos {
BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM / 1_GeV);
- const int iBeam = corsika::epos::getEposXSCode(
+ const int iBeam = epos::getEposXSCode(
BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
- if (!iBeam)
- throw std::runtime_error(
- "calcCrossSectionCoM: interaction of beam hadron not defined in "
- "Epos!");
CORSIKA_LOGGER_TRACE(logger_,
"projectile cross section type={} "
@@ -280,10 +287,6 @@ namespace corsika::epos {
else if (iBeam == 3)
initializeEventCoM(Code::KPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM);
- else
- throw std::runtime_error(
- "calcCrossSectionCoM: interaction of beam hadron not defined in "
- "Epos!");
double sigProd, sigEla = 0;
float sigTot1, sigProd1, sigCut1 = 0;
@@ -306,10 +309,10 @@ namespace corsika::epos {
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 {
+ inline std::tuple<CrossSectionType, CrossSectionType>
+ InteractionModel::readCrossSectionTableLab(Code const BeamId, int const BeamA,
+ int const BeamZ, Code const TargetId,
+ HEPEnergyType const EnergyLab) const {
CORSIKA_LOGGER_DEBUG(logger_,
"readCrossSectionTableLab: input: "
"beamId={}, "
@@ -329,7 +332,7 @@ namespace corsika::epos {
Ekin = (EnergyLab / Abeam - constants::nucleonMass) / 1_GeV;
} else {
::epos::hadr2_.idproj = convertToEposRaw(BeamId);
- int const iBeam = corsika::epos::getEposXSCode(
+ int const iBeam = epos::getEposXSCode(
BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
CORSIKA_LOGGER_TRACE(logger_,
"projectile cross section type={} "
@@ -340,13 +343,6 @@ namespace corsika::epos {
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 = 1;
if (is_nucleus(TargetId)) { Atarget = get_nucleus_A(TargetId); }
@@ -366,219 +362,145 @@ namespace corsika::epos {
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 {
+ inline std::tuple<CrossSectionType, CrossSectionType>
+ InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ auto const sqrtS = sqrt(sqrtS2);
+
+ if (!isValid(projectileId, targetId, sqrtS)) {
+ return {CrossSectionType::zero(), CrossSectionType::zero()};
+ }
+ HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+ static_pow<2>(get_mass(targetId))) /
+ (2 * get_mass(targetId));
+ int beamA = 1;
+ int beamZ = 1;
+ if (is_nucleus(projectileId)) {
+ beamA = get_nucleus_A(projectileId);
+ beamZ = get_nucleus_Z(projectileId);
+ }
+
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);
+ " target={}"
+ " ELab={:4.3f} GeV, sqrtS={}",
+ projectileId, beamA, beamZ, targetId, Elab / 1_GeV,
+ sqrtS / 1_GeV);
+ return readCrossSectionTableLab(projectileId, beamA, beamZ, targetId, Elab);
}
- template <>
- inline corsika::GrammageType Interaction::getInteractionLength(
- SetupParticle const& projectile) const {
-
- const corsika::Code corsikaBeamId = projectile.getPID();
- const bool kInteraction = corsika::epos::canInteract(corsikaBeamId);
- CORSIKA_LOGGER_DEBUG(logger_,
- "InteractionLength: input: \n"
- " energy: {} GeV "
- " beam can interact: {} "
- " beam pid: {}",
- projectile.getEnergy() / 1_GeV, kInteraction,
- projectile.getPID());
-
- if (kInteraction) {
-
- // define projectile nuclei
- int beamA = 1;
- int beamZ = 1;
- if (is_nucleus(corsikaBeamId)) {
- beamA = get_nucleus_A(corsikaBeamId);
- beamZ = get_nucleus_Z(corsikaBeamId);
- }
+ template <typename TSecondaryView>
+ inline void InteractionModel::doInteraction(TSecondaryView& view,
+ Code const projectileId,
+ Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) {
+
+ count_ = count_ + 1;
+
+ // define projectile
+ // define projectile, in lab frame
+ auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ auto const sqrtS = sqrt(sqrtS2);
+ if (!isValid(projectileId, targetId, sqrtS)) {
+ throw std::runtime_error("invalid projectiel/target/energy combination.");
+ }
+ HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+ static_pow<2>(get_mass(targetId))) /
+ (2 * get_mass(targetId));
- // get target from environment
- MomentumVector const& pLab = projectile.getMomentum();
- CoordinateSystemPtr const& labCS = pLab.getCoordinateSystem();
+ // system of initial-state
+ COMBoost boost(projectileP4, targetP4);
- // assume target is at rest!!
- MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV});
+ auto const& originalCS = boost.getOriginalCS();
+ auto const& csPrime =
+ boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
- // total momentum and energy
- HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass;
+ HEPMomentumType const pLabMag =
+ sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId)));
+ MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
- auto const* currentNode = projectile.getNode();
- const auto& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
+ // internal EPOS lab system
+ COMBoost boostInternal({Elab, pLab}, get_mass(targetId));
- si::CrossSectionType weightedProdCrossSection = mediumComposition.getWeightedSum(
- [=](corsika::Code targetID) -> si::CrossSectionType {
- return std::get<0>(this->getCrossSectionLab(corsikaBeamId, beamA, beamZ,
- targetID, get_nucleus_A(targetID),
- get_nucleus_Z(targetID), Elab));
- });
+ CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: {} interaction, Elab={} ", projectileId,
+ Elab);
- CORSIKA_LOGGER_DEBUG(logger_, "InteractionLength: weighted CrossSection (mb): {} ",
- weightedProdCrossSection / 1_mb);
+ int beamA = 1;
+ int beamZ = 1;
+ if (is_nucleus(projectileId)) {
+ beamA = get_nucleus_A(projectileId);
+ beamZ = get_nucleus_Z(projectileId);
+ CORSIKA_LOGGER_DEBUG(logger_, "A={}, Z={} ", beamA, beamZ);
+ }
- // calculate interaction length in medium
- GrammageType const int_length = mediumComposition.getAverageMassNumber() *
- constants::u / weightedProdCrossSection;
- CORSIKA_LOGGER_DEBUG(logger_, "interaction length (g/cm2): {} ",
- int_length / (0.001_kg) * 1_cm * 1_cm);
+ HEPMomentumType const projectileMomentumLabPerNucleon = pLabMag / beamA;
- return int_length;
+ // // from corsika7 interface
+ // // NEXLNK-part
+ int targetA = 1;
+ int targetZ = 1;
+ if (is_nucleus(targetId)) {
+ targetA = get_nucleus_A(targetId);
+ targetZ = get_nucleus_Z(targetId);
}
+ initializeEventLab(projectileId, beamA, beamZ, targetId, targetA, targetZ,
+ projectileMomentumLabPerNucleon);
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
+ // create event
+ int iarg = 1;
+ ::epos::aepos_(iarg);
- template <typename TSecondaryView>
- inline void Interaction::doInteraction(TSecondaryView& view) {
+ ::epos::afinal_();
- auto const projectile = view.getProjectile();
- auto const corsikaBeamId = projectile.getPID();
-
- CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: {} interaction, Elab={} ",
- corsikaBeamId, projectile.getEnergy());
-
- if (corsika::epos::canInteract(corsikaBeamId)) {
- count_ = count_ + 1;
- // 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;
- int beamZ = 1;
- if (is_nucleus(corsikaBeamId)) {
- beamA = get_nucleus_A(corsikaBeamId);
- beamZ = get_nucleus_Z(corsikaBeamId);
- CORSIKA_LOGGER_DEBUG(logger_, "A={}, Z={} ", beamA, beamZ);
- }
-
- HEPEnergyType const projectileMomentumLabPerNucleon = projectileMomentum / beamA;
-
- // define target
-
- // 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] = getCrossSectionLab(
- corsikaBeamId, beamA, beamZ, targetId, get_nucleus_A(targetId),
- get_nucleus_Z(targetId), eProjectileLab);
- cross_section_of_components[i] = sigProd;
- }
-
- const auto targetCode =
- mediumComposition.sampleTarget(cross_section_of_components, RNG_);
- CORSIKA_LOGGER_DEBUG(logger_, "target selected: {} ", targetCode);
-
- // // from corsika7 interface
- // // NEXLNK-part
- int targetA = 1;
- int targetZ = 1;
- if (is_nucleus(targetCode)) {
- targetA = get_nucleus_A(targetCode);
- targetZ = get_nucleus_Z(targetCode);
- }
- initializeEventLab(corsikaBeamId, beamA, beamZ, targetCode, targetA, targetZ,
- projectileMomentumLabPerNucleon);
+ if (epos_listing_) {
+ char nam[9] = "EPOSLHC&";
+ ::epos::alistf_(nam, 9);
+ }
- // create event
- int iarg = 1;
- ::epos::aepos_(iarg);
+ // NSTORE-part
- ::epos::afinal_();
+ MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV;
- if (epos_listing_) {
- char nam[9] = "EPOSLHC&";
- ::epos::alistf_(nam, 9);
- }
+ // position and time of interaction, not used in QgsjetII
+ auto const projectile = view.getProjectile();
+ Point const pOrig = projectile.getPosition();
+ TimeType const tOrig = projectile.getTime();
+
+ // secondaries
+ EposStack es;
+ CORSIKA_LOGGER_DEBUG(logger_, "number of particles: {}", es.getSize());
+ for (auto& psec : es) {
+ if (!psec.isFinal()) continue;
+
+ auto momentum = psec.getMomentum(csPrime);
+ // this is not "CoM" here, but rather the system defined by projectile+target,
+ // which in Cascade-mode is already lab
+ auto const P4com = boostInternal.toCoM(FourVector{psec.getEnergy(), momentum});
+ auto const P4output = boost.fromCoM(P4com);
+ auto p3output = P4output.getSpaceLikeComponents();
+ p3output.rebase(originalCS); // transform back into standard lab frame
+
+ EposCode const eposId = psec.getPID();
+ Code const pid = epos::convertFromEpos(eposId);
+ CORSIKA_LOGGER_TRACE(logger_,
+ " id= {}"
+ " p= {}",
+ pid, p3output.getComponents() / 1_GeV);
- // NSTORE-part
-
- MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV;
-
- // secondaries
- EposStack es;
- CORSIKA_LOGGER_DEBUG(logger_, "number of particles: {}", es.getSize());
- for (auto& psec : es) {
- if (!psec.isFinal()) continue;
-
- auto momentum = psec.getMomentum(zAxisFrame);
-
- momentum.rebase(originalCS); // transform back into standard lab frame
-
- EposCode const eposId = psec.getPID();
- Code const pid = corsika::epos::convertFromEpos(eposId);
- CORSIKA_LOGGER_TRACE(logger_,
- " id= {}"
- " p= {}",
- pid, momentum.getComponents() / 1_GeV);
- if (!is_nucleus(pid)) {
- auto pnew = view.addSecondary(std::make_tuple(pid, momentum, pOrig, tOrig));
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
- } else {
- unsigned int A = 0;
- unsigned int Z = 0;
- if (pid == Code::Deuterium) {
- A = 2;
- Z = 1;
- } else if (pid == Code::Tritium) {
- A = 3;
- Z = 1;
- } else if (pid == Code::Helium) {
- A = 4;
- Z = 2;
- } else {
- Z = get_nucleus_Z(eposId);
- A = get_nucleus_Z(eposId);
- }
- auto pnew = view.addSecondary(
- std::make_tuple(get_nucleus_code(A, Z), momentum, pOrig, tOrig));
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
- }
- }
- CORSIKA_LOGGER_DEBUG(
- logger_,
- "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
- ", Elab_final={}"
- ", Plab_final={}",
- Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
- } else
- CORSIKA_LOGGER_WARN(
- logger_, "Projectile not configured for interaction! This is likely an error!");
+ auto pnew = view.addSecondary(std::make_tuple(pid, p3output, pOrig, tOrig));
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
+ }
+ CORSIKA_LOGGER_DEBUG(
+ logger_,
+ "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
+ ", Elab_final={}"
+ ", Plab_final={}",
+ Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
}
} // namespace corsika::epos
diff --git a/corsika/detail/modules/qgsjetII/InteractionModel.inl b/corsika/detail/modules/qgsjetII/InteractionModel.inl
index b65d5b3f7..8203faa49 100644
--- a/corsika/detail/modules/qgsjetII/InteractionModel.inl
+++ b/corsika/detail/modules/qgsjetII/InteractionModel.inl
@@ -41,8 +41,10 @@ namespace corsika::qgsjetII {
CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_);
}
- inline bool InteractionModel::isValid(Code const projectileId,
- Code const targetId) const {
+ inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+ HEPEnergyType const sqrtS) const {
+
+ if (sqrtS < sqrtSmin_) { return false; }
if (is_nucleus(targetId)) {
size_t iTarget = get_nucleus_A(targetId);
if (iTarget > int(maxMassNumber_) || iTarget <= 0) { return false; }
@@ -66,10 +68,11 @@ namespace corsika::qgsjetII {
if (!corsika::qgsjetII::canInteract(projectileId)) {
return CrossSectionType::zero();
}
- if (!isValid(projectileId, targetId)) { return CrossSectionType::zero(); }
// define projectile, in lab frame
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ auto const sqrtS = sqrt(sqrtS2);
+ if (!isValid(projectileId, targetId, sqrtS)) { return CrossSectionType::zero(); }
HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
@@ -101,15 +104,13 @@ namespace corsika::qgsjetII {
"doInteraction: {} interaction possible? {}",
projectileId, corsika::qgsjetII::canInteract(projectileId));
+ // define projectile, in lab frame
+ auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ auto const sqrtS = sqrt(sqrtS2);
if (!corsika::qgsjetII::canInteract(projectileId) ||
- !isValid(projectileId, targetId)) {
+ !isValid(projectileId, targetId, sqrtS)) {
throw std::runtime_error("invalid target/projectile/energy combination.");
}
-
- CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
-
- // define projectile, in lab frame
- auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
@@ -150,6 +151,8 @@ namespace corsika::qgsjetII {
qgini_(Elab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
qgconf_();
+ CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
+
// bookkeeping
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
@@ -168,7 +171,6 @@ namespace corsika::qgsjetII {
auto const& csPrime =
boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
- MomentumVector const projectileMomentum = projectileP4.getSpaceLikeComponents();
HEPMomentumType const pLabMag =
sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId)));
MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
diff --git a/corsika/modules/epos/Interaction.hpp b/corsika/modules/epos/Interaction.hpp
index 56900ca01..96306b8fa 100644
--- a/corsika/modules/epos/Interaction.hpp
+++ b/corsika/modules/epos/Interaction.hpp
@@ -10,20 +10,19 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <tuple>
namespace corsika::epos {
- class Interaction : public InteractionProcess<Interaction> {
- std::string data_path_;
- unsigned int count_ = 0;
- bool epos_listing_;
+ class InteractionModel {
public:
- Interaction(std::string const& dataPath = "", bool const epos_printout_on = false);
- ~Interaction();
+ InteractionModel(std::string const& dataPath = "",
+ bool const epos_printout_on = false);
+ ~InteractionModel();
//! returns production and elastic cross section for hadrons in epos. Inputs are:
//! CorsikaId of beam particle, CorsikaId of target particle, center-of-mass energy.
@@ -41,27 +40,39 @@ namespace corsika::epos {
//! 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;
-
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle const&) const;
+ std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+ Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
/**
- In this function EPOSLHC is called to produce one event. The
- event is copied into the shower lab frame.
+ * Checks validity of projectile, target and energy combination.
*/
- template <typename TSecondaries>
- void doInteraction(TSecondaries&);
+ bool isValid(Code const projectileId, Code const targetId,
+ HEPEnergyType const sqrtS) const;
- bool isValidCoMEnergy(HEPEnergyType const ecm) const {
- return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
+ /**
+ * Get the inelatic/production cross section.
+ *
+ * @param projectileId
+ * @param targetId
+ * @param projectileP4
+ * @param targetP4
+ * @return CrossSectionType
+ */
+ CrossSectionType getCrossSection(Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ return std::get<0>(
+ getCrossSectionInelEla(projectileId, targetId, projectileP4, targetP4));
}
- //! eposlhc only accepts nuclei with X<=A<=Y as targets, or protons aka Hydrogen or
- //! neutrons (p,n == nucleon)
- bool isValidTarget(Code const) const;
+ /**
+ * In this function EPOSLHC is called to produce one event. The
+ * event is copied into the shower lab frame.
+ */
+ template <typename TSecondaries>
+ void doInteraction(TSecondaries&, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
void initialize() const;
void initializeEventCoM(Code const, int const, int const, Code const, int const,
@@ -73,6 +84,10 @@ namespace corsika::epos {
void setParticlesStable() const;
private:
+ std::string data_path_;
+ unsigned int count_ = 0;
+ bool epos_listing_;
+
default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("epos");
std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_epos_Interaction");
HEPEnergyType const minEnergyCoM_ = 6 * 1e9 * electronvolt;
diff --git a/corsika/modules/qgsjetII/InteractionModel.hpp b/corsika/modules/qgsjetII/InteractionModel.hpp
index 45e9f8ef0..0e3896df5 100644
--- a/corsika/modules/qgsjetII/InteractionModel.hpp
+++ b/corsika/modules/qgsjetII/InteractionModel.hpp
@@ -2,7 +2,7 @@
* (c) Copyright 2020 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
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version ofp
* the license.
*/
@@ -33,7 +33,7 @@ namespace corsika::qgsjetII {
* @param beamId
* @param targetId
*/
- bool isValid(Code const beamId, Code const targetId) const;
+ bool isValid(Code const beamId, Code const targetId, HEPEnergyType const sqrtS) const;
/**
* Return the QGSJETII inelastic/production cross section.
@@ -70,6 +70,7 @@ namespace corsika::qgsjetII {
corsika::default_prng_type& rng_ =
corsika::RNGManager<>::getInstance().getRandomStream("qgsjet");
static size_t constexpr maxMassNumber_ = 208;
+ static HEPEnergyType constexpr sqrtSmin_ = 10_GeV;
};
} // namespace corsika::qgsjetII
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index 3eab3eaaf..40d8690b6 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -29,7 +29,7 @@
using namespace corsika;
using namespace corsika::epos;
-TEST_CASE("epos", "module,process") {
+TEST_CASE("Epos", "module,process") {
logging::set_level(logging::level::trace);
@@ -134,75 +134,99 @@ TEST_CASE("EposInterface", "modules") {
CHECK(rndm < 1);
}
- SECTION("InteractionInterface - valid targets") {
+ SECTION("InteractionInterface - isValid") {
- Interaction model;
- CHECK_FALSE(model.isValidTarget(Code::Electron));
- CHECK(model.isValidTarget(Code::Hydrogen));
- CHECK(model.isValidTarget(Code::Helium));
- CHECK_FALSE(model.isValidTarget(Code::Iron));
- CHECK(model.isValidTarget(Code::Oxygen));
+ InteractionModel model;
+
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 100_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Helium, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Iron, 100_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Oxygen, 100_GeV));
+ }
+
+ SECTION("InteractionInterface - getCrossSectionInelEla") {
+
+ InteractionModel model;
// 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);
+ auto const [xs_prod_pp, xs_ela_pp] = model.getCrossSectionInelEla(
+ Code::Proton, Code::Proton,
+ {sqrt(static_pow<2>(100_GeV) + static_pow<2>(Proton::mass)),
+ {cs, 100_GeV, 0_GeV, 0_GeV}},
+ {Proton::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
+
+ auto const [xs_prod_pn, xs_ela_pn] = model.getCrossSectionInelEla(
+ Code::Proton, Code::Neutron,
+ {sqrt(static_pow<2>(100_GeV) + static_pow<2>(Proton::mass)),
+ {cs, 100_GeV, 0_GeV, 0_GeV}},
+ {Neutron::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
+
+ auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] = model.getCrossSectionInelEla(
+ Code::Proton, Code::Hydrogen,
+ {sqrt(static_pow<2>(100_GeV) + static_pow<2>(Proton::mass)),
+ {cs, 100_GeV, 0_GeV, 0_GeV}},
+ {Hydrogen::mass, {cs, 0_GeV, 0_GeV, 0_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") {
+ SECTION("InteractionModelInterface - hadron cross sections") {
- Interaction model;
+ InteractionModel model;
// p-p at 7TeV around 70mb according to LHC
- auto const [xs_prod, xs_ela] =
- model.getCrossSectionLab(Code::Proton, 1, 1, Code::Proton, 1, 1,
- sqs2elab(7_TeV, Proton::mass, Proton::mass));
+ auto const xs_prod = model.getCrossSection(
+ Code::Proton, Code::Proton,
+ {3.5_TeV,
+ {cs, sqrt(static_pow<2>(3.5_TeV) - static_pow<2>(Proton::mass)), 0_GeV, 0_GeV}},
+ {3.5_TeV,
+ {cs, -sqrt(static_pow<2>(3.5_TeV) - static_pow<2>(Proton::mass)), 0_GeV,
+ 0_GeV}});
CHECK(xs_prod / 1_mb == Approx(70.7).margin(2.1));
- { [[maybe_unused]] auto const& dum_xs = xs_ela; }
// pi-n at 7TeV
- auto const [xs_prod1, xs_ela1] =
- model.getCrossSectionLab(Code::PiPlus, 0, 0, Code::Neutron, 1, 0,
- sqs2elab(7_TeV, PiPlus::mass, Neutron::mass));
+ auto const xs_prod1 = model.getCrossSection(
+ Code::PiPlus, Code::Neutron,
+ {3.5_TeV,
+ {cs, sqrt(static_pow<2>(3.5_TeV) - static_pow<2>(PiPlus::mass)), 0_GeV, 0_GeV}},
+ {3.5_TeV,
+ {cs, -sqrt(static_pow<2>(3.5_TeV) - static_pow<2>(Neutron::mass)), 0_GeV,
+ 0_GeV}});
CHECK(xs_prod1 / 1_mb == Approx(52.7).margin(2.1));
- { [[maybe_unused]] auto const& dum_xs = xs_ela1; }
// k-p at 7TeV
- auto const [xs_prod2, xs_ela2] =
- model.getCrossSectionLab(Code::KPlus, 0, 0, Code::Proton, 1, 1,
- sqs2elab(7_TeV, KPlus::mass, Proton::mass));
+ auto const xs_prod2 = model.getCrossSection(
+ Code::KPlus, Code::Proton,
+ {3.5_TeV,
+ {cs, sqrt(static_pow<2>(3.5_TeV) - static_pow<2>(KPlus::mass)), 0_GeV, 0_GeV}},
+ {3.5_TeV,
+ {cs, -sqrt(static_pow<2>(3.5_TeV) - static_pow<2>(Proton::mass)), 0_GeV,
+ 0_GeV}});
CHECK(xs_prod2 / 1_mb == Approx(45.7).margin(2.1));
- { [[maybe_unused]] auto const& dum_xs = xs_ela2; }
}
SECTION("InteractionInterface - nuclear cross sections") {
- Interaction model;
+ InteractionModel model;
- auto const [xs_prod, xs_ela] = model.getCrossSectionLab(
- Code::Proton, 1, 1, Code::Oxygen, Oxygen::nucleus_A, Oxygen::nucleus_Z, 100_GeV);
+ auto const xs_prod = model.getCrossSection(
+ Code::Proton, Code::Oxygen,
+ {100_GeV,
+ {cs, sqrt(static_pow<2>(100_GeV) - static_pow<2>(Proton::mass)), 0_GeV, 0_GeV}},
+ {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
CHECK(xs_prod / 1_mb == Approx(287.0).margin(5.1));
- { [[maybe_unused]] auto const& dum_xs = xs_ela; }
- auto const [xs_prod2, xs_ela2] = model.getCrossSectionLab(
- Code::Nitrogen, Nitrogen::nucleus_A, Nitrogen::nucleus_Z, Code::Oxygen,
- Oxygen::nucleus_A, Oxygen::nucleus_Z, 400_GeV);
+ auto const xs_prod2 = model.getCrossSection(
+ Code::Nitrogen, Code::Oxygen,
+ {400_GeV,
+ {cs, sqrt(static_pow<2>(400_GeV) - static_pow<2>(Nitrogen::mass)), 0_GeV,
+ 0_GeV}},
+ {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
CHECK(xs_prod2 / 1_mb == Approx(1076.7).margin(3.1));
- { [[maybe_unused]] auto const& dum_xs = xs_ela2; }
-
- // 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);
- { [[maybe_unused]] auto const& dum_xs = xs_ela3; }
}
SECTION("InteractionInterface - low energy") {
@@ -214,10 +238,11 @@ TEST_CASE("EposInterface", "modules") {
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr;
- auto particle = stack->first();
-
- Interaction model;
- model.doInteraction(view);
+ InteractionModel model;
+ model.doInteraction(
+ view, Code::Proton, Code::Oxygen,
+ {sqrt(static_pow<2>(P0) + static_pow<2>(Proton::mass)), {cs, P0, 0_GeV, 0_GeV}},
+ {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
auto const pSum = sumMomentum(view, cs);
@@ -230,26 +255,29 @@ TEST_CASE("EposInterface", "modules") {
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(93.3).margin(0.1));
+ // [[maybe_unused]] const GrammageType length =
+ // model.getInteractionLength(particle);
+ // CHECK(length / 1_g * 1_cm * 1_cm == Approx(93.3).margin(0.1));
}
SECTION("InteractionInterface - nuclear projectile") {
- const HEPEnergyType P0 = 10_TeV;
+ HEPEnergyType const P0 = 10_TeV;
+ Code const pid = get_nucleus_code(8, 4);
auto [stack, viewPtr] = setup::testing::setup_stack(
- get_nucleus_code(8, 4), P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ pid, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
MomentumVector plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr;
- auto particle = stack->first();
-
- Interaction model;
+ InteractionModel model;
#ifndef __clang__
- // This is very obscure since it fails for -O2, but for both clang and gcc ???
- model.doInteraction(view);
+ // @todo This is very obscure since it fails for -O2, but for both clang and gcc ???
+ model.doInteraction(
+ view, pid, Code::Oxygen,
+ {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), {cs, P0, 0_GeV, 0_GeV}},
+ {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
auto const pSum = sumMomentum(view, cs);
@@ -263,8 +291,9 @@ TEST_CASE("EposInterface", "modules") {
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
#endif
- [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
- CHECK(length / 1_g * 1_cm * 1_cm ==
- Approx(30).margin(20)); // this is no physics validation
+ // [[maybe_unused]] const GrammageType length =
+ // model.getInteractionLength(particle);
+ // CHECK(length / 1_g * 1_cm * 1_cm ==
+ // Approx(30).margin(20)); // this is no physics validation
}
}
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 40f544798..966928f9a 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -96,8 +96,9 @@ TEST_CASE("QgsjetII", "[processes]") {
corsika::qgsjetII::InteractionModel model;
- CHECK_FALSE(model.isValid(Code::Electron, Code::Proton));
- CHECK_FALSE(model.isValid(Code::Proton, Code::Electron));
+ CHECK_FALSE(model.isValid(Code::Electron, Code::Proton, 1_TeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 1_TeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Proton, 1_GeV));
}
}
--
GitLab