IAP GITLAB

Skip to content
Snippets Groups Projects
Commit ec68cdff authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

also adapted EPOS

parent 3704a385
No related branches found
No related tags found
1 merge request!387Hadronic model interface overhaul
Hide whitespace changes
Inline Side-by-side
corsika/detail/modules/epos/Interaction.inl
+ 163
− 241
View file @ ec68cdff
...@@ -25,13 +25,10 @@ ...@@ -25,13 +25,10 @@
#include <string> #include <string>
#include <tuple> #include <tuple>
using namespace corsika;
using SetupParticle = setup::Stack::stack_iterator_type;
namespace corsika::epos { namespace corsika::epos {
inline Interaction::Interaction(std::string const& dataPath, inline InteractionModel::InteractionModel(std::string const& dataPath,
bool const epos_printout_on) bool const epos_printout_on)
: data_path_(dataPath) : data_path_(dataPath)
, epos_listing_(epos_printout_on) { , epos_listing_(epos_printout_on) {
if (dataPath == "") { if (dataPath == "") {
...@@ -46,7 +43,7 @@ namespace corsika::epos { ...@@ -46,7 +43,7 @@ namespace corsika::epos {
setParticlesStable(); setParticlesStable();
} }
inline void Interaction::setParticlesStable() const { inline void InteractionModel::setParticlesStable() const {
CORSIKA_LOGGER_DEBUG(logger_, CORSIKA_LOGGER_DEBUG(logger_,
"set all particles known to CORSIKA stable inside EPOS.."); "set all particles known to CORSIKA stable inside EPOS..");
for (auto& p : get_all_particles()) { for (auto& p : get_all_particles()) {
...@@ -59,11 +56,22 @@ namespace corsika::epos { ...@@ -59,11 +56,22 @@ namespace corsika::epos {
} }
} }
inline bool Interaction::isValidTarget(Code const TargetId) const { inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxTargetMassNumber_); 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..."); CORSIKA_LOGGER_DEBUG(logger_, "initializing...");
...@@ -137,10 +145,10 @@ namespace corsika::epos { ...@@ -137,10 +145,10 @@ namespace corsika::epos {
Argon::nucleus_A, Argon::nucleus_Z, 100_GeV); Argon::nucleus_A, Argon::nucleus_Z, 100_GeV);
} }
inline void Interaction::initializeEventCoM(Code const idBeam, int const iBeamA, inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
int const iBeamZ, Code const idTarget, int const iBeamZ, Code const idTarget,
int const iTargetA, int const iTargetZ, int const iTargetA, int const iTargetZ,
HEPEnergyType const Ecm) const { HEPEnergyType const Ecm) const {
CORSIKA_LOGGER_TRACE(logger_, CORSIKA_LOGGER_TRACE(logger_,
"initialize event in CoM frame!" "initialize event in CoM frame!"
" Ecm={}", " Ecm={}",
...@@ -163,10 +171,10 @@ namespace corsika::epos { ...@@ -163,10 +171,10 @@ namespace corsika::epos {
::epos::ainit_(); ::epos::ainit_();
} }
inline void Interaction::initializeEventLab(Code const idBeam, int const iBeamA, inline void InteractionModel::initializeEventLab(Code const idBeam, int const iBeamA,
int const iBeamZ, Code const idTarget, int const iBeamZ, Code const idTarget,
int const iTargetA, int const iTargetZ, int const iTargetA, int const iTargetZ,
HEPEnergyType const Plab) const { HEPEnergyType const Plab) const {
CORSIKA_LOGGER_TRACE(logger_, CORSIKA_LOGGER_TRACE(logger_,
"initialize event in lab. frame!" "initialize event in lab. frame!"
" Plab per nuc={} GeV", " Plab per nuc={} GeV",
...@@ -191,10 +199,10 @@ namespace corsika::epos { ...@@ -191,10 +199,10 @@ namespace corsika::epos {
::epos::ainit_(); ::epos::ainit_();
} }
inline void Interaction::configureParticles(Code const idBeam, int const iBeamA, inline void InteractionModel::configureParticles(Code const idBeam, int const iBeamA,
int const iBeamZ, Code const idTarget, int const iBeamZ, Code const idTarget,
int const iTargetA, int const iTargetA,
int const iTargetZ) const { int const iTargetZ) const {
CORSIKA_LOGGER_TRACE(logger_, CORSIKA_LOGGER_TRACE(logger_,
"setting " "setting "
"Beam={}, " "Beam={}, "
...@@ -227,9 +235,8 @@ namespace corsika::epos { ...@@ -227,9 +235,8 @@ namespace corsika::epos {
::epos::hadr25_.idtargin = convertToEposRaw(Code::Neutron); ::epos::hadr25_.idtargin = convertToEposRaw(Code::Neutron);
::epos::nucl1_.matarg = 1; ::epos::nucl1_.matarg = 1;
::epos::nucl1_.latarg = -1; ::epos::nucl1_.latarg = -1;
} else {
throw std::runtime_error("Epos: target outside range!");
} }
CORSIKA_LOGGER_TRACE(logger_, CORSIKA_LOGGER_TRACE(logger_,
"inside EPOS: " "inside EPOS: "
"Id beam={}, " "Id beam={}, "
...@@ -246,11 +253,15 @@ namespace corsika::epos { ...@@ -246,11 +253,15 @@ namespace corsika::epos {
::epos::nucl1_.matarg, ::epos::had10_.icltar); ::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( inline std::tuple<CrossSectionType, CrossSectionType>
Code const BeamId, int const BeamA, int const BeamZ, Code const TargetId, InteractionModel::calcCrossSectionCoM(Code const BeamId, int const BeamA,
int const TargetA, int const TargetZ, const HEPEnergyType EnergyCOM) const { int const BeamZ, Code const TargetId,
int const TargetA, int const TargetZ,
const HEPEnergyType EnergyCOM) const {
CORSIKA_LOGGER_DEBUG(logger_, CORSIKA_LOGGER_DEBUG(logger_,
"calcCrossSection: input:" "calcCrossSection: input:"
" beamId={}, beamA={}, beamZ={}" " beamId={}, beamA={}, beamZ={}"
...@@ -259,12 +270,8 @@ namespace corsika::epos { ...@@ -259,12 +270,8 @@ namespace corsika::epos {
BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ, BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM / 1_GeV); 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) 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_, CORSIKA_LOGGER_TRACE(logger_,
"projectile cross section type={} " "projectile cross section type={} "
...@@ -280,10 +287,6 @@ namespace corsika::epos { ...@@ -280,10 +287,6 @@ namespace corsika::epos {
else if (iBeam == 3) else if (iBeam == 3)
initializeEventCoM(Code::KPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ, initializeEventCoM(Code::KPlus, BeamA, BeamZ, TargetId, TargetA, TargetZ,
EnergyCOM); EnergyCOM);
else
throw std::runtime_error(
"calcCrossSectionCoM: interaction of beam hadron not defined in "
"Epos!");
double sigProd, sigEla = 0; double sigProd, sigEla = 0;
float sigTot1, sigProd1, sigCut1 = 0; float sigTot1, sigProd1, sigCut1 = 0;
...@@ -306,10 +309,10 @@ namespace corsika::epos { ...@@ -306,10 +309,10 @@ namespace corsika::epos {
return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb); return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
} }
inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType> inline std::tuple<CrossSectionType, CrossSectionType>
Interaction::readCrossSectionTableLab(Code const BeamId, int const BeamA, InteractionModel::readCrossSectionTableLab(Code const BeamId, int const BeamA,
int const BeamZ, Code const TargetId, int const BeamZ, Code const TargetId,
HEPEnergyType const EnergyLab) const { HEPEnergyType const EnergyLab) const {
CORSIKA_LOGGER_DEBUG(logger_, CORSIKA_LOGGER_DEBUG(logger_,
"readCrossSectionTableLab: input: " "readCrossSectionTableLab: input: "
"beamId={}, " "beamId={}, "
...@@ -329,7 +332,7 @@ namespace corsika::epos { ...@@ -329,7 +332,7 @@ namespace corsika::epos {
Ekin = (EnergyLab / Abeam - constants::nucleonMass) / 1_GeV; Ekin = (EnergyLab / Abeam - constants::nucleonMass) / 1_GeV;
} else { } else {
::epos::hadr2_.idproj = convertToEposRaw(BeamId); ::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) BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
CORSIKA_LOGGER_TRACE(logger_, CORSIKA_LOGGER_TRACE(logger_,
"projectile cross section type={} " "projectile cross section type={} "
...@@ -340,13 +343,6 @@ namespace corsika::epos { ...@@ -340,13 +343,6 @@ namespace corsika::epos {
Abeam = 1; Abeam = 1;
Ekin = (EnergyLab - get_mass(BeamId)) / 1_GeV; 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; int Atarget = 1;
if (is_nucleus(TargetId)) { Atarget = get_nucleus_A(TargetId); } if (is_nucleus(TargetId)) { Atarget = get_nucleus_A(TargetId); }
...@@ -366,219 +362,145 @@ namespace corsika::epos { ...@@ -366,219 +362,145 @@ namespace corsika::epos {
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> inline std::tuple<CrossSectionType, CrossSectionType>
Interaction::getCrossSectionLab(corsika::Code const BeamId, int const BeamA, InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
int const BeamZ, corsika::Code const TargetId, FourMomentum const& projectileP4,
int const TargetA, int const TargetZ, FourMomentum const& targetP4) const {
const corsika::HEPEnergyType EnergyLab) 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_, CORSIKA_LOGGER_DEBUG(logger_,
"getCrossSectionLab: input:" "getCrossSectionLab: input:"
" beamId={}, beamA={}, beamZ={}" " beamId={}, beamA={}, beamZ={}"
" target={}, targetA={}, targetZ={}" " target={}"
" ELab={:4.3f} GeV,", " ELab={:4.3f} GeV, sqrtS={}",
BeamId, BeamA, BeamZ, TargetId, TargetA, TargetZ, projectileId, beamA, beamZ, targetId, Elab / 1_GeV,
EnergyLab / 1_GeV); sqrtS / 1_GeV);
return readCrossSectionTableLab(BeamId, BeamA, BeamZ, TargetId, EnergyLab); return readCrossSectionTableLab(projectileId, beamA, beamZ, targetId, Elab);
} }
template <> template <typename TSecondaryView>
inline corsika::GrammageType Interaction::getInteractionLength( inline void InteractionModel::doInteraction(TSecondaryView& view,
SetupParticle const& projectile) const { Code const projectileId,
Code const targetId,
const corsika::Code corsikaBeamId = projectile.getPID(); FourMomentum const& projectileP4,
const bool kInteraction = corsika::epos::canInteract(corsikaBeamId); FourMomentum const& targetP4) {
CORSIKA_LOGGER_DEBUG(logger_,
"InteractionLength: input: \n" count_ = count_ + 1;
" energy: {} GeV "
" beam can interact: {} " // define projectile
" beam pid: {}", // define projectile, in lab frame
projectile.getEnergy() / 1_GeV, kInteraction, auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
projectile.getPID()); auto const sqrtS = sqrt(sqrtS2);
if (!isValid(projectileId, targetId, sqrtS)) {
if (kInteraction) { throw std::runtime_error("invalid projectiel/target/energy combination.");
}
// define projectile nuclei HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
int beamA = 1; static_pow<2>(get_mass(targetId))) /
int beamZ = 1; (2 * get_mass(targetId));
if (is_nucleus(corsikaBeamId)) {
beamA = get_nucleus_A(corsikaBeamId);
beamZ = get_nucleus_Z(corsikaBeamId);
}
// get target from environment // system of initial-state
MomentumVector const& pLab = projectile.getMomentum(); COMBoost boost(projectileP4, targetP4);
CoordinateSystemPtr const& labCS = pLab.getCoordinateSystem();
// assume target is at rest!! auto const& originalCS = boost.getOriginalCS();
MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV}); auto const& csPrime =
boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
// total momentum and energy HEPMomentumType const pLabMag =
HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass; sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId)));
MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
auto const* currentNode = projectile.getNode(); // internal EPOS lab system
const auto& mediumComposition = COMBoost boostInternal({Elab, pLab}, get_mass(targetId));
currentNode->getModelProperties().getNuclearComposition();
si::CrossSectionType weightedProdCrossSection = mediumComposition.getWeightedSum( CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: {} interaction, Elab={} ", projectileId,
[=](corsika::Code targetID) -> si::CrossSectionType { Elab);
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): {} ", int beamA = 1;
weightedProdCrossSection / 1_mb); 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 HEPMomentumType const projectileMomentumLabPerNucleon = pLabMag / beamA;
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);
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> ::epos::afinal_();
inline void Interaction::doInteraction(TSecondaryView& view) {
auto const projectile = view.getProjectile(); if (epos_listing_) {
auto const corsikaBeamId = projectile.getPID(); char nam[9] = "EPOSLHC&";
::epos::alistf_(nam, 9);
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);
// create event // NSTORE-part
int iarg = 1;
::epos::aepos_(iarg);
::epos::afinal_(); MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
if (epos_listing_) { // position and time of interaction, not used in QgsjetII
char nam[9] = "EPOSLHC&"; auto const projectile = view.getProjectile();
::epos::alistf_(nam, 9); 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 auto pnew = view.addSecondary(std::make_tuple(pid, p3output, pOrig, tOrig));
Plab_final += pnew.getMomentum();
MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV}); Elab_final += pnew.getEnergy();
HEPEnergyType Elab_final = 0_GeV; }
CORSIKA_LOGGER_DEBUG(
// secondaries logger_,
EposStack es; "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
CORSIKA_LOGGER_DEBUG(logger_, "number of particles: {}", es.getSize()); ", Elab_final={}"
for (auto& psec : es) { ", Plab_final={}",
if (!psec.isFinal()) continue; Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
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!");
} }
} // namespace corsika::epos } // namespace corsika::epos
corsika/detail/modules/qgsjetII/InteractionModel.inl
+ 12
− 10
View file @ ec68cdff
...@@ -41,8 +41,10 @@ namespace corsika::qgsjetII { ...@@ -41,8 +41,10 @@ namespace corsika::qgsjetII {
CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_); CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_);
} }
inline bool InteractionModel::isValid(Code const projectileId, inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
Code const targetId) const { HEPEnergyType const sqrtS) const {
if (sqrtS < sqrtSmin_) { return false; }
if (is_nucleus(targetId)) { if (is_nucleus(targetId)) {
size_t iTarget = get_nucleus_A(targetId); size_t iTarget = get_nucleus_A(targetId);
if (iTarget > int(maxMassNumber_) || iTarget <= 0) { return false; } if (iTarget > int(maxMassNumber_) || iTarget <= 0) { return false; }
...@@ -66,10 +68,11 @@ namespace corsika::qgsjetII { ...@@ -66,10 +68,11 @@ namespace corsika::qgsjetII {
if (!corsika::qgsjetII::canInteract(projectileId)) { if (!corsika::qgsjetII::canInteract(projectileId)) {
return CrossSectionType::zero(); return CrossSectionType::zero();
} }
if (!isValid(projectileId, targetId)) { return CrossSectionType::zero(); }
// define projectile, in lab frame // define projectile, in lab frame
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr(); 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)) - HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) / static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId)); (2 * get_mass(targetId));
...@@ -101,15 +104,13 @@ namespace corsika::qgsjetII { ...@@ -101,15 +104,13 @@ namespace corsika::qgsjetII {
"doInteraction: {} interaction possible? {}", "doInteraction: {} interaction possible? {}",
projectileId, corsika::qgsjetII::canInteract(projectileId)); 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) || if (!corsika::qgsjetII::canInteract(projectileId) ||
!isValid(projectileId, targetId)) { !isValid(projectileId, targetId, sqrtS)) {
throw std::runtime_error("invalid target/projectile/energy combination."); 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)) - HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) / static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId)); (2 * get_mass(targetId));
...@@ -150,6 +151,8 @@ namespace corsika::qgsjetII { ...@@ -150,6 +151,8 @@ namespace corsika::qgsjetII {
qgini_(Elab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber); qgini_(Elab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
qgconf_(); qgconf_();
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
// bookkeeping // bookkeeping
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV}); MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV; HEPEnergyType Elab_final = 0_GeV;
...@@ -168,7 +171,6 @@ namespace corsika::qgsjetII { ...@@ -168,7 +171,6 @@ namespace corsika::qgsjetII {
auto const& csPrime = auto const& csPrime =
boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade) boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
MomentumVector const projectileMomentum = projectileP4.getSpaceLikeComponents();
HEPMomentumType const pLabMag = HEPMomentumType const pLabMag =
sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId))); sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId)));
MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag}); MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
......
corsika/modules/epos/Interaction.hpp
+ 36
− 21
View file @ ec68cdff
...@@ -10,20 +10,19 @@ ...@@ -10,20 +10,19 @@
#include <corsika/framework/core/ParticleProperties.hpp> #include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/random/RNGManager.hpp> #include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/process/InteractionProcess.hpp> #include <corsika/framework/process/InteractionProcess.hpp>
#include <tuple> #include <tuple>
namespace corsika::epos { namespace corsika::epos {
class Interaction : public InteractionProcess<Interaction> { class InteractionModel {
std::string data_path_;
unsigned int count_ = 0;
bool epos_listing_;
public: public:
Interaction(std::string const& dataPath = "", bool const epos_printout_on = false); InteractionModel(std::string const& dataPath = "",
~Interaction(); bool const epos_printout_on = false);
~InteractionModel();
//! returns production and elastic cross section for hadrons in epos. Inputs are: //! returns production and elastic cross section for hadrons in epos. Inputs are:
//! CorsikaId of beam particle, CorsikaId of target particle, center-of-mass energy. //! CorsikaId of beam particle, CorsikaId of target particle, center-of-mass energy.
...@@ -41,27 +40,39 @@ namespace corsika::epos { ...@@ -41,27 +40,39 @@ namespace corsika::epos {
//! returns production and elastic cross section. Allowed configurations are //! returns production and elastic cross section. Allowed configurations are
//! hadron-nucleon, hadron-nucleus and nucleus-nucleus. Inputs are particle id's mass //! hadron-nucleon, hadron-nucleus and nucleus-nucleus. Inputs are particle id's mass
//! and charge numbers and total energy in the lab. //! and charge numbers and total energy in the lab.
std::tuple<CrossSectionType, CrossSectionType> getCrossSectionLab( std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
Code const, int const, int const, Code const, int const, int const, Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
HEPEnergyType const) const; FourMomentum const& targetP4) const;
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
/** /**
In this function EPOSLHC is called to produce one event. The * Checks validity of projectile, target and energy combination.
event is copied into the shower lab frame.
*/ */
template <typename TSecondaries> bool isValid(Code const projectileId, Code const targetId,
void doInteraction(TSecondaries&); 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) * In this function EPOSLHC is called to produce one event. The
bool isValidTarget(Code const) const; * 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 initialize() const;
void initializeEventCoM(Code const, int const, int const, Code const, int const, void initializeEventCoM(Code const, int const, int const, Code const, int const,
...@@ -73,6 +84,10 @@ namespace corsika::epos { ...@@ -73,6 +84,10 @@ namespace corsika::epos {
void setParticlesStable() const; void setParticlesStable() const;
private: private:
std::string data_path_;
unsigned int count_ = 0;
bool epos_listing_;
default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("epos"); default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("epos");
std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_epos_Interaction"); std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_epos_Interaction");
HEPEnergyType const minEnergyCoM_ = 6 * 1e9 * electronvolt; HEPEnergyType const minEnergyCoM_ = 6 * 1e9 * electronvolt;
......
corsika/modules/qgsjetII/InteractionModel.hpp
+ 3
− 2
View file @ ec68cdff
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * 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. * the license.
*/ */
...@@ -33,7 +33,7 @@ namespace corsika::qgsjetII { ...@@ -33,7 +33,7 @@ namespace corsika::qgsjetII {
* @param beamId * @param beamId
* @param targetId * @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. * Return the QGSJETII inelastic/production cross section.
...@@ -70,6 +70,7 @@ namespace corsika::qgsjetII { ...@@ -70,6 +70,7 @@ namespace corsika::qgsjetII {
corsika::default_prng_type& rng_ = corsika::default_prng_type& rng_ =
corsika::RNGManager<>::getInstance().getRandomStream("qgsjet"); corsika::RNGManager<>::getInstance().getRandomStream("qgsjet");
static size_t constexpr maxMassNumber_ = 208; static size_t constexpr maxMassNumber_ = 208;
static HEPEnergyType constexpr sqrtSmin_ = 10_GeV;
}; };
} // namespace corsika::qgsjetII } // namespace corsika::qgsjetII
......
tests/modules/testEpos.cpp
+ 88
− 59
View file @ ec68cdff
...@@ -29,7 +29,7 @@ ...@@ -29,7 +29,7 @@
using namespace corsika; using namespace corsika;
using namespace corsika::epos; using namespace corsika::epos;
TEST_CASE("epos", "module,process") { TEST_CASE("Epos", "module,process") {
logging::set_level(logging::level::trace); logging::set_level(logging::level::trace);
...@@ -134,75 +134,99 @@ TEST_CASE("EposInterface", "modules") { ...@@ -134,75 +134,99 @@ TEST_CASE("EposInterface", "modules") {
CHECK(rndm < 1); CHECK(rndm < 1);
} }
SECTION("InteractionInterface - valid targets") { SECTION("InteractionInterface - isValid") {
Interaction model; InteractionModel model;
CHECK_FALSE(model.isValidTarget(Code::Electron));
CHECK(model.isValidTarget(Code::Hydrogen)); CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 100_GeV));
CHECK(model.isValidTarget(Code::Helium)); CHECK(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV));
CHECK_FALSE(model.isValidTarget(Code::Iron)); CHECK(model.isValid(Code::Proton, Code::Helium, 100_GeV));
CHECK(model.isValidTarget(Code::Oxygen)); 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 // hydrogen target == proton target == neutron target
auto const [xs_prod_pp, xs_ela_pp] = auto const [xs_prod_pp, xs_ela_pp] = model.getCrossSectionInelEla(
model.getCrossSectionLab(Code::Proton, 1, 1, Code::Proton, 1, 1, 100_GeV); Code::Proton, Code::Proton,
auto const [xs_prod_pn, xs_ela_pn] = {sqrt(static_pow<2>(100_GeV) + static_pow<2>(Proton::mass)),
model.getCrossSectionLab(Code::Proton, 1, 1, Code::Neutron, 1, 0, 100_GeV); {cs, 100_GeV, 0_GeV, 0_GeV}},
auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] = {Proton::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
model.getCrossSectionLab(Code::Proton, 1, 1, Code::Hydrogen, 1, 1, 100_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_pHydrogen);
CHECK(xs_prod_pp == xs_prod_pn); CHECK(xs_prod_pp == xs_prod_pn);
CHECK(xs_ela_pp == xs_ela_pHydrogen); CHECK(xs_ela_pp == xs_ela_pHydrogen);
CHECK(xs_ela_pn == 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 // p-p at 7TeV around 70mb according to LHC
auto const [xs_prod, xs_ela] = auto const xs_prod = model.getCrossSection(
model.getCrossSectionLab(Code::Proton, 1, 1, Code::Proton, 1, 1, Code::Proton, Code::Proton,
sqs2elab(7_TeV, Proton::mass, Proton::mass)); {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)); CHECK(xs_prod / 1_mb == Approx(70.7).margin(2.1));
{ [[maybe_unused]] auto const& dum_xs = xs_ela; }
// pi-n at 7TeV // pi-n at 7TeV
auto const [xs_prod1, xs_ela1] = auto const xs_prod1 = model.getCrossSection(
model.getCrossSectionLab(Code::PiPlus, 0, 0, Code::Neutron, 1, 0, Code::PiPlus, Code::Neutron,
sqs2elab(7_TeV, PiPlus::mass, Neutron::mass)); {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)); CHECK(xs_prod1 / 1_mb == Approx(52.7).margin(2.1));
{ [[maybe_unused]] auto const& dum_xs = xs_ela1; }
// k-p at 7TeV // k-p at 7TeV
auto const [xs_prod2, xs_ela2] = auto const xs_prod2 = model.getCrossSection(
model.getCrossSectionLab(Code::KPlus, 0, 0, Code::Proton, 1, 1, Code::KPlus, Code::Proton,
sqs2elab(7_TeV, KPlus::mass, Proton::mass)); {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)); CHECK(xs_prod2 / 1_mb == Approx(45.7).margin(2.1));
{ [[maybe_unused]] auto const& dum_xs = xs_ela2; }
} }
SECTION("InteractionInterface - nuclear cross sections") { SECTION("InteractionInterface - nuclear cross sections") {
Interaction model; InteractionModel model;
auto const [xs_prod, xs_ela] = model.getCrossSectionLab( auto const xs_prod = model.getCrossSection(
Code::Proton, 1, 1, Code::Oxygen, Oxygen::nucleus_A, Oxygen::nucleus_Z, 100_GeV); 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)); 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( auto const xs_prod2 = model.getCrossSection(
Code::Nitrogen, Nitrogen::nucleus_A, Nitrogen::nucleus_Z, Code::Oxygen, Code::Nitrogen, Code::Oxygen,
Oxygen::nucleus_A, Oxygen::nucleus_Z, 400_GeV); {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)); 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") { SECTION("InteractionInterface - low energy") {
...@@ -214,10 +238,11 @@ TEST_CASE("EposInterface", "modules") { ...@@ -214,10 +238,11 @@ TEST_CASE("EposInterface", "modules") {
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr; setup::StackView& view = *viewPtr;
auto particle = stack->first(); InteractionModel model;
model.doInteraction(
Interaction model; view, Code::Proton, Code::Oxygen,
model.doInteraction(view); {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); auto const pSum = sumMomentum(view, cs);
...@@ -230,26 +255,29 @@ TEST_CASE("EposInterface", "modules") { ...@@ -230,26 +255,29 @@ TEST_CASE("EposInterface", "modules") {
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV)); Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05)); CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle); // [[maybe_unused]] const GrammageType length =
CHECK(length / 1_g * 1_cm * 1_cm == Approx(93.3).margin(0.1)); // model.getInteractionLength(particle);
// CHECK(length / 1_g * 1_cm * 1_cm == Approx(93.3).margin(0.1));
} }
SECTION("InteractionInterface - nuclear projectile") { 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( 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 plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr; setup::StackView& view = *viewPtr;
auto particle = stack->first(); InteractionModel model;
Interaction model;
#ifndef __clang__ #ifndef __clang__
// This is very obscure since it fails for -O2, but for both clang and gcc ??? // @todo This is very obscure since it fails for -O2, but for both clang and gcc ???
model.doInteraction(view); 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); auto const pSum = sumMomentum(view, cs);
...@@ -263,8 +291,9 @@ TEST_CASE("EposInterface", "modules") { ...@@ -263,8 +291,9 @@ TEST_CASE("EposInterface", "modules") {
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV)); Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05)); CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
#endif #endif
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle); // [[maybe_unused]] const GrammageType length =
CHECK(length / 1_g * 1_cm * 1_cm == // model.getInteractionLength(particle);
Approx(30).margin(20)); // this is no physics validation // CHECK(length / 1_g * 1_cm * 1_cm ==
// Approx(30).margin(20)); // this is no physics validation
} }
} }
tests/modules/testQGSJetII.cpp
+ 3
− 2
View file @ ec68cdff
...@@ -96,8 +96,9 @@ TEST_CASE("QgsjetII", "[processes]") { ...@@ -96,8 +96,9 @@ TEST_CASE("QgsjetII", "[processes]") {
corsika::qgsjetII::InteractionModel model; corsika::qgsjetII::InteractionModel model;
CHECK_FALSE(model.isValid(Code::Electron, Code::Proton)); CHECK_FALSE(model.isValid(Code::Electron, Code::Proton, 1_TeV));
CHECK_FALSE(model.isValid(Code::Proton, Code::Electron)); CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 1_TeV));
CHECK_FALSE(model.isValid(Code::Proton, Code::Proton, 1_GeV));
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum