From d7ac9f8a97065295e7a1e95c57c52feddf789813 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 25 Oct 2021 13:30:01 +0200
Subject: [PATCH] qgsjetII new interface
---
.../detail/modules/qgsjetII/Interaction.inl | 397 ------------------
.../modules/qgsjetII/InteractionModel.inl | 308 ++++++++++++++
.../modules/sibyll/InteractionModel.inl | 35 +-
corsika/framework/utility/COMBoost.hpp | 32 +-
corsika/modules/LongitudinalProfile.hpp | 3 +-
corsika/modules/QGSJetII.hpp | 21 +-
corsika/modules/Sibyll.hpp | 4 +-
corsika/modules/qgsjetII/Interaction.hpp | 65 ---
corsika/modules/qgsjetII/InteractionModel.hpp | 77 ++++
.../modules/qgsjetII/ParticleConversion.hpp | 8 +-
corsika/modules/sibyll/InteractionModel.hpp | 4 +-
tests/modules/CMakeLists.txt | 2 +-
tests/modules/testQGSJetII.cpp | 103 +++--
tests/modules/testSibyll.cpp | 7 +-
14 files changed, 514 insertions(+), 552 deletions(-)
delete mode 100644 corsika/detail/modules/qgsjetII/Interaction.inl
create mode 100644 corsika/detail/modules/qgsjetII/InteractionModel.inl
delete mode 100644 corsika/modules/qgsjetII/Interaction.hpp
create mode 100644 corsika/modules/qgsjetII/InteractionModel.hpp
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
deleted file mode 100644
index c3147dbb8..000000000
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ /dev/null
@@ -1,397 +0,0 @@
-/*
- * (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
- * the license.
- */
-
-#include <corsika/modules/qgsjetII/Interaction.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/framework/geometry/QuantityVector.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
-#include <corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp>
-#include <corsika/modules/qgsjetII/QGSJetIIStack.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-
-#include <sstream>
-#include <string>
-#include <tuple>
-
-#include <qgsjet-II-04.hpp>
-
-namespace corsika::qgsjetII {
-
- inline Interaction::Interaction(boost::filesystem::path dataPath) {
- CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
-
- // initialize QgsjetII
- static bool initialized = false;
- if (!initialized) {
- qgset_();
- datadir DIR(dataPath.string() + "/");
- qgaini_(DIR.data);
- initialized = true;
- }
- }
-
- inline Interaction::~Interaction() {
- CORSIKA_LOG_DEBUG("QgsjetII::Interaction n= {}", count_);
- }
-
- inline CrossSectionType Interaction::getCrossSection(const Code beamId,
- const Code targetId,
- const HEPEnergyType Elab,
- const unsigned int Abeam,
- const unsigned int targetA) const {
- double sigProd = std::numeric_limits<double>::infinity();
-
- if (corsika::qgsjetII::canInteract(beamId)) {
-
- int const iBeam = static_cast<QgsjetIIXSClassIntType>(
- corsika::qgsjetII::getQgsjetIIXSCode(beamId));
- int iTarget = 1;
- if (is_nucleus(targetId)) {
- iTarget = targetA;
- if (iTarget > int(maxMassNumber_) || iTarget <= 0) {
- std::ostringstream txt;
- txt << "QgsjetII target outside range. Atarget=" << iTarget;
- throw std::runtime_error(txt.str().c_str());
- }
- }
- int iProjectile = 1;
- if (is_nucleus(beamId)) {
- iProjectile = Abeam;
- if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) {
- std::ostringstream txt;
- txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
- throw std::runtime_error(txt.str().c_str());
- }
- }
-
- CORSIKA_LOG_DEBUG(
- "QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
- " iProjectile= {} iTarget= {}",
- Elab / 1_GeV, iBeam, iProjectile, iTarget);
- sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
- CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
- }
-
- return sigProd * 1_mb;
- }
-
- template <typename TParticle>
- inline GrammageType Interaction::getInteractionLength(const TParticle& particle) const {
-
- // coordinate system, get global frame of reference
- CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
-
- const Code corsikaBeamId = particle.getPID();
-
- // beam particles for qgsjetII : 1, 2, 3 for p, pi, k
- // read from cross section code table
- const bool kInteraction = corsika::qgsjetII::canInteract(corsikaBeamId);
-
- // FOR NOW: assume target is at rest
- MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
-
- // total momentum and energy
- HEPEnergyType const Elab = particle.getEnergy();
-
- CORSIKA_LOG_DEBUG(
- "Interaction: LambdaInt: \n"
- " input energy: {} GeV"
- " beam can interact: {}"
- " beam pid: {}",
- particle.getEnergy() / 1_GeV, kInteraction, corsikaBeamId);
-
- if (kInteraction) {
-
- int Abeam = 0;
- if (is_nucleus(corsikaBeamId)) Abeam = get_nucleus_A(corsikaBeamId);
-
- // get target from environment
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
-
- auto const* currentNode = particle.getNode();
- const auto& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
-
- CrossSectionType weightedProdCrossSection =
- mediumComposition.getWeightedSum([=](Code targetID) -> CrossSectionType {
- int targetA = 0;
- if (is_nucleus(targetID)) targetA = get_nucleus_A(targetID);
- return getCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA);
- });
-
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- "IntLength: weighted CrossSection (mb): {}",
- weightedProdCrossSection / 1_mb);
-
- // calculate interaction length in medium
- GrammageType const int_length = mediumComposition.getAverageMassNumber() *
- constants::u / weightedProdCrossSection;
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- "interaction length (g/cm2): {}",
- int_length / (0.001_kg) * 1_cm * 1_cm);
-
- return int_length;
- }
-
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
-
- /**
- In this function QGSJETII is called to produce one event. The
- event is copied (and boosted) into the shower lab frame.
- */
-
- template <typename TView>
- inline void Interaction::doInteraction(TView& view) {
-
- auto const projectile = view.getProjectile();
- auto const corsikaBeamId = projectile.getPID();
- CORSIKA_LOG_DEBUG(
- "ProcessQgsjetII: "
- "doInteraction: {} interaction possible? {}",
- corsikaBeamId, corsika::qgsjetII::canInteract(corsikaBeamId));
-
- if (!corsika::qgsjetII::canInteract(corsikaBeamId)) return;
-
- CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
-
- // position and time of interaction, not used in QgsjetII
- Point const pOrig = projectile.getPosition();
- TimeType const tOrig = projectile.getTime();
-
- // define target
- // for QgsjetII is always a single nucleon
- // FOR NOW: target is always at rest
- auto const targetEnergyLab = 0_GeV + constants::nucleonMass;
- auto const targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
- FourVector const PtargLab(targetEnergyLab, targetMomentumLab);
-
- // define projectile
- HEPEnergyType const projectileEnergyLab = projectile.getEnergy();
- auto const projectileMomentumLab = projectile.getMomentum();
-
- int beamA = 0;
- if (is_nucleus(corsikaBeamId)) beamA = get_nucleus_A(corsikaBeamId);
-
- HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
-
- CORSIKA_LOG_DEBUG(
- "ebeam lab: {} GeV "
- "pbeam lab: {} GeV ",
- projectileEnergyLab / 1_GeV, projectileMomentumLab.getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG(
- "etarget lab: {} GeV "
- "ptarget lab: {} GeV ",
- targetEnergyLab / 1_GeV, targetMomentumLab.getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("position of interaction: {}", pOrig.getCoordinates());
- CORSIKA_LOG_DEBUG("time: {} ", tOrig);
-
- // 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
- */
- 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];
- int targetA = 0;
- if (is_nucleus(targetId)) targetA = get_nucleus_A(targetId);
- const auto sigProd =
- getCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA);
- cross_section_of_components[i] = sigProd;
- }
-
- const auto targetCode =
- mediumComposition.sampleTarget(cross_section_of_components, rng_);
-
- int targetMassNumber = 1; // proton
- if (is_nucleus(targetCode)) { // nucleus
- targetMassNumber = get_nucleus_A(targetCode);
- if (targetMassNumber > int(maxMassNumber_))
- throw std::runtime_error(
- "QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
- // allowed path here
- } else {
- if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
- throw std::runtime_error(
- "QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
- // here
- }
- CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetCode, targetMassNumber);
-
- int projectileMassNumber = 1; // "1" means "hadron"
- QgsjetIIHadronType qgsjet_hadron_type =
- qgsjetII::getQgsjetIIHadronType(corsikaBeamId);
- if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
- projectileMassNumber = get_nucleus_A(corsikaBeamId);
- if (projectileMassNumber > int(maxMassNumber_))
- throw std::runtime_error(
- "QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
- // allowed path here
- std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
- QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
- std::uniform_int_distribution select(0, 1);
- qgsjet_hadron_type = nucleons[select(rng_)];
- } else {
- // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
- if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
- qgsjet_hadron_type = alternate_;
- alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
- ? QgsjetIIHadronType::PiMinusType
- : QgsjetIIHadronType::PiPlusType);
- }
- }
-
- // beam id for qgsjetII
- int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
- if (!is_nucleus(corsikaBeamId)) {
- kBeam = corsika::qgsjetII::convertToQgsjetIIRaw(corsikaBeamId);
- // from conex
- if (kBeam == 0) { // replace pi0 or rho0 with pi+/pi-
- static int select = 1;
- kBeam = select;
- select *= -1;
- }
- // replace lambda by neutron
- if (kBeam == 6)
- kBeam = 3;
- else if (kBeam == -6)
- kBeam = -3;
- // else if (abs(kBeam)>6) -> throw
- }
-
- count_++;
- int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
- CORSIKA_LOG_DEBUG(
- "qgsjet_hadron_type_int={} projectileMassNumber={} targetMassNumber={}",
- qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
- qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
- targetMassNumber);
- qgconf_();
-
- // bookkeeping
- MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV;
-
- // to read the secondaries
- // define rotation to and from CoM frame
- // CoM frame definition in QgsjetII projectile: +z
- auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
- CoordinateSystemPtr const zAxisFrame =
- make_rotationToZ(originalCS, projectileMomentumLab);
-
- // fragments
- QGSJetIIFragmentsStack qfs;
- for (auto& fragm : qfs) {
- int const A = fragm.getFragmentSize();
- if (A == 1) { // nucleon
- std::uniform_real_distribution<double> select;
- Code idFragm = Code::Proton;
- if (select(rng_) > 0.5) { idFragm = Code::Neutron; }
-
- const HEPMassType nucleonMass = get_mass(idFragm);
- // no pT, frgments just go forward
- auto momentum =
- Vector(zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
- 0.0_GeV, 0.0_GeV,
- sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
- (projectileEnergyLabPerNucleon - nucleonMass))});
-
- momentum.rebase(originalCS); // transform back into standard lab frame
- CORSIKA_LOG_DEBUG(
- "secondary fragment> id= {}"
- " p={}",
- idFragm, momentum.getComponents());
- auto pnew = view.addSecondary(std::make_tuple(idFragm, momentum, pOrig, tOrig));
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
-
- } else { // nucleus, A>1
-
- int Z = 0;
- switch (A) {
- case 2: // deuterium
- Z = 1;
- break;
- case 3: // tritium
- Z = 1;
- break;
- case 4: // helium
- Z = 2;
- break;
- default: // nucleus
- {
- Z = int(A / 2.15 + 0.7);
- }
- }
-
- HEPMassType const nucleusMass = Proton::mass * Z + Neutron::mass * (A - Z);
- // no pT, frgments just go forward
- auto momentum = Vector(
- zAxisFrame, QuantityVector<hepmomentum_d>{
- 0.0_GeV, 0.0_GeV,
- sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
- (projectileEnergyLabPerNucleon * A - nucleusMass))});
-
- momentum.rebase(originalCS); // transform back into standard lab frame
- CORSIKA_LOG_DEBUG(
- "secondary fragment> id={}"
- " p={}"
- " A={}"
- " Z={}",
- get_nucleus_code(A, Z), momentum.getComponents(), A, Z);
-
- auto pnew = view.addSecondary(
- std::make_tuple(get_nucleus_code(A, Z), momentum, pOrig, tOrig));
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
- }
- }
-
- // secondaries
- QGSJetIIStack qs;
- for (auto& psec : qs) {
-
- auto momentum = psec.getMomentum(zAxisFrame);
-
- momentum.rebase(originalCS); // transform back into standard lab frame
- CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}",
- corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
- momentum.getComponents());
- auto pnew = view.addSecondary(std::make_tuple(
- corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), 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={}"
- ", N_wounded,targ={}"
- ", N_wounded,proj={}"
- ", N_fragm,proj={}",
- Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents(),
- QGSJetIIFragmentsStackData::getWoundedNucleonsTarget(),
- QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile(), qfs.getSize());
- }
-} // namespace corsika::qgsjetII
diff --git a/corsika/detail/modules/qgsjetII/InteractionModel.inl b/corsika/detail/modules/qgsjetII/InteractionModel.inl
new file mode 100644
index 000000000..3d5209d8e
--- /dev/null
+++ b/corsika/detail/modules/qgsjetII/InteractionModel.inl
@@ -0,0 +1,308 @@
+/*
+ * (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
+ * the license.
+ */
+
+#include <corsika/modules/qgsjetII/InteractionModel.hpp>
+
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+#include <corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp>
+#include <corsika/modules/qgsjetII/QGSJetIIStack.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+
+#include <sstream>
+#include <tuple>
+
+#include <qgsjet-II-04.hpp>
+
+namespace corsika::qgsjetII {
+
+ inline InteractionModel::InteractionModel(boost::filesystem::path const dataPath) {
+ CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
+
+ // initialize QgsjetII
+ static bool initialized = false;
+ if (!initialized) {
+ qgset_();
+ datadir DIR(dataPath.string() + "/");
+ qgaini_(DIR.data);
+ initialized = true;
+ }
+ }
+
+ inline InteractionModel::~InteractionModel() {
+ CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_);
+ }
+
+ inline void InteractionModel::isValid(Code const projectileId,
+ Code const targetId) const {
+ if (is_nucleus(targetId)) {
+ size_t iTarget = get_nucleus_A(targetId);
+ if (iTarget > int(maxMassNumber_) || iTarget <= 0) {
+ std::ostringstream txt;
+ txt << "QgsjetII target outside range. Atarget=" << iTarget;
+ throw std::runtime_error(txt.str().c_str());
+ }
+ } else if (targetId != Proton::code) {
+ std::ostringstream txt;
+ txt << "QgsjetII not valid target=" << targetId;
+ throw std::runtime_error(txt.str().c_str());
+ }
+
+ if (is_nucleus(projectileId)) {
+ size_t iProjectile = get_nucleus_A(projectileId);
+ if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) {
+ std::ostringstream txt;
+ txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
+ throw std::runtime_error(txt.str().c_str());
+ }
+ } else if (!is_hadron(projectileId)) {
+ std::ostringstream txt;
+ txt << "QgsjetII projectile can only be hadrons. Is=" << projectileId;
+ throw std::runtime_error(txt.str().c_str());
+ }
+ }
+
+ inline CrossSectionType InteractionModel::getCrossSection(
+ Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+
+ double sigProd = std::numeric_limits<double>::infinity();
+
+ if (!corsika::qgsjetII::canInteract(projectileId)) { return sigProd * 1_mb; }
+
+ isValid(projectileId, targetId); // throws
+
+ // 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));
+
+ int const iBeam = static_cast<QgsjetIIXSClassIntType>(
+ corsika::qgsjetII::getQgsjetIIXSCode(projectileId));
+ int iTarget = 1;
+ if (is_nucleus(targetId)) { iTarget = get_nucleus_A(targetId); }
+ int iProjectile = 1;
+ if (is_nucleus(projectileId)) { iProjectile = get_nucleus_A(projectileId); }
+
+ CORSIKA_LOG_DEBUG(
+ "QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
+ " iProjectile= {} iTarget= {}",
+ Elab / 1_GeV, iBeam, iProjectile, iTarget);
+ sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
+ CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
+
+ return sigProd * 1_mb;
+ }
+
+ template <typename TSecondaries>
+ inline void InteractionModel::doInteraction(TSecondaries& view, Code const projectileId,
+ Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) {
+
+ CORSIKA_LOG_DEBUG(
+ "ProcessQgsjetII: "
+ "doInteraction: {} interaction possible? {}",
+ projectileId, corsika::qgsjetII::canInteract(projectileId));
+
+ if (!corsika::qgsjetII::canInteract(projectileId)) return;
+
+ isValid(projectileId, targetId); // throws
+
+ 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));
+
+ int beamA = 0;
+ if (is_nucleus(projectileId)) { beamA = get_nucleus_A(projectileId); }
+
+ CORSIKA_LOG_DEBUG("ebeam lab: {} GeV ", Elab / 1_GeV);
+
+ int targetMassNumber = 1; // proton
+ if (is_nucleus(targetId)) { // nucleus
+ targetMassNumber = get_nucleus_A(targetId);
+ }
+ CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetId, targetMassNumber);
+
+ // select QGSJetII internal projectile type
+ int projectileMassNumber = 1; // "1" means "hadron"
+ QgsjetIIHadronType qgsjet_hadron_type = qgsjetII::getQgsjetIIHadronType(projectileId);
+ if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
+ projectileMassNumber = get_nucleus_A(projectileId);
+ std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
+ QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
+ std::uniform_int_distribution select(0, 1);
+ qgsjet_hadron_type = nucleons[select(rng_)];
+ } else if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
+ // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
+ qgsjet_hadron_type = alternate_;
+ alternate_ =
+ (alternate_ == QgsjetIIHadronType::PiPlusType ? QgsjetIIHadronType::PiMinusType
+ : QgsjetIIHadronType::PiPlusType);
+ }
+
+ count_++;
+ int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
+ CORSIKA_LOG_DEBUG(
+ "qgsjet_hadron_type_int={} projectileMassNumber={} targetMassNumber={}",
+ qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
+ qgini_(Elab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
+ qgconf_();
+
+ // bookkeeping
+ MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV;
+
+ // to read the secondaries
+ // define rotation to and from CoM frame
+ // CoM frame definition in QgsjetII projectile: +z
+
+ // QGSJetII, both, in input and output only considers the lab frame with a target at
+ // rest.
+
+ // system of initial-state
+ COMBoost boost(projectileP4, targetP4);
+
+ auto const& originalCS = boost.getOriginalCS();
+ 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});
+
+ // internal QGSJetII system
+ COMBoost boostInternal({Elab, pLab}, get_mass(targetId));
+
+ // position and time of interaction, not used in QgsjetII
+ auto const projectile = view.getProjectile();
+ Point const pOrig = projectile.getPosition();
+ TimeType const tOrig = projectile.getTime();
+
+ // fragments
+ QGSJetIIFragmentsStack qfs;
+ for (auto& fragm : qfs) {
+ int const A = fragm.getFragmentSize();
+ if (A == 1) { // nucleon
+ std::uniform_real_distribution<double> select;
+ Code idFragm = Code::Proton;
+ if (select(rng_) > 0.5) { idFragm = Code::Neutron; }
+
+ const HEPMassType nucleonMass = get_mass(idFragm);
+ // no pT, fragments just go forward
+ HEPEnergyType const projectileEnergyLabPerNucleon = Elab / beamA;
+ MomentumVector momentum{csPrime,
+ {0.0_GeV, 0.0_GeV,
+ sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
+ (projectileEnergyLabPerNucleon - nucleonMass))}};
+
+ // 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{projectileEnergyLabPerNucleon, momentum});
+ auto const P4output = boost.fromCoM(P4com);
+ auto p3output = P4output.getSpaceLikeComponents();
+ p3output.rebase(originalCS); // transform back into standard lab frame
+
+ CORSIKA_LOG_DEBUG(
+ "secondary fragment> id= {}"
+ " p={}",
+ idFragm, p3output.getComponents());
+ auto pnew = view.addSecondary(std::make_tuple(idFragm, p3output, pOrig, tOrig));
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
+
+ } else { // nucleus, A>1
+
+ int Z = 0;
+ switch (A) {
+ case 2: // deuterium
+ Z = 1;
+ break;
+ case 3: // tritium
+ Z = 1;
+ break;
+ case 4: // helium
+ Z = 2;
+ break;
+ default: // nucleus
+ {
+ Z = int(A / 2.15 + 0.7);
+ }
+ }
+
+ HEPMassType const nucleusMass = Proton::mass * Z + Neutron::mass * (A - Z);
+ // no pT, frgments just go forward
+ HEPEnergyType const projectileEnergyLabPerNucleon = Elab / beamA;
+ MomentumVector momentum{
+ csPrime,
+ {0.0_GeV, 0.0_GeV,
+ sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
+ (projectileEnergyLabPerNucleon * A - nucleusMass))}};
+
+ // 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{projectileEnergyLabPerNucleon * A, momentum});
+ auto const P4output = boost.fromCoM(P4com);
+ auto p3output = P4output.getSpaceLikeComponents();
+ p3output.rebase(originalCS); // transform back into standard lab frame
+
+ CORSIKA_LOG_DEBUG(
+ "secondary fragment> id={}"
+ " p={}"
+ " A={}"
+ " Z={}",
+ get_nucleus_code(A, Z), p3output.getComponents(), A, Z);
+
+ auto pnew = view.addSecondary(
+ std::make_tuple(get_nucleus_code(A, Z), p3output, pOrig, tOrig));
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
+ }
+ }
+
+ // secondaries
+ QGSJetIIStack qs;
+ for (auto& psec : qs) {
+
+ 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
+
+ CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}",
+ corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
+ p3output.getComponents());
+ auto pnew = view.addSecondary(std::make_tuple(
+ corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), p3output, 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={}"
+ ", N_wounded,targ={}"
+ ", N_wounded,proj={}"
+ ", N_fragm,proj={}",
+ Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents(),
+ QGSJetIIFragmentsStackData::getWoundedNucleonsTarget(),
+ QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile(), qfs.getSize());
+ }
+} // namespace corsika::qgsjetII
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index c0716be3a..af1a01133 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -156,8 +156,8 @@ namespace corsika::sibyll {
auto const tmp = psib.getMomentum().getComponents();
auto const pCoM = MomentumVector(csPrime, tmp);
HEPEnergyType const eCoM = psib.getEnergy();
- auto const Plab = boost.fromCoM(FourVector{eCoM, pCoM});
- auto const p3lab = Plab.getSpaceLikeComponents();
+ auto const P4lab = boost.fromCoM(FourVector{eCoM, pCoM});
+ auto const p3lab = P4lab.getSpaceLikeComponents();
// add to corsika stack
auto pnew = secondaries.addSecondary(std::make_tuple(
@@ -167,19 +167,20 @@ namespace corsika::sibyll {
Elab_final += pnew.getEnergy();
Ecm_final += psib.getEnergy();
}
- HEPEnergyType const Elab_initial =
- static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
- CORSIKA_LOG_DEBUG(
- "conservation (all GeV): "
- "sqrtSnn={}, sqrtSnn_final={}, "
- "Elab_initial={}, Elab_final={}, "
- "diff(%)={}, "
- "E in nucleons={}, "
- "Plab_final={} ",
- sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Elab_initial,
- Elab_final / 1_GeV, (Elab_final - Elab_initial) / Elab_initial * 100,
- constants::nucleonMass * get_nwounded() / 1_GeV,
- (Plab_final / 1_GeV).getComponents());
- } // namespace corsika::sibyll
-
+ { // just output
+ HEPEnergyType const Elab_initial =
+ static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
+ CORSIKA_LOG_DEBUG(
+ "conservation (all GeV): "
+ "sqrtSnn={}, sqrtSnn_final={}, "
+ "Elab_initial={}, Elab_final={}, "
+ "diff(%)={}, "
+ "E in nucleons={}, "
+ "Plab_final={} ",
+ sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Elab_initial,
+ Elab_final / 1_GeV, (Elab_final - Elab_initial) / Elab_initial * 100,
+ constants::nucleonMass * get_nwounded() / 1_GeV,
+ (Plab_final / 1_GeV).getComponents());
+ }
+ }
} // namespace corsika::sibyll
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index 8bba6b43d..98be9fc07 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -25,11 +25,11 @@ namespace corsika {
*/
/**
- * @class COMBoost
+ * @class COMBoost
* @ingroup Utilities
*
- * This utility class handles Lorentz boost between different
- * referenence frames, using FourVector.
+ * This utility class handles Lorentz boost (in one spatial direction)
+ * between different referenence frames, using FourVector.
*
* The class is initialized with projectile and optionally target
* energy/momentum data. During initialization, a rotation matrix is
@@ -47,16 +47,26 @@ namespace corsika {
class COMBoost {
public:
- //! construct a COMBoost given four-vector of projectile and mass of target (target at
- //! rest)
+ /**
+ * Construct a COMBoost given four-vector of projectile and mass of target (target at
+ * rest).
+ *
+ * The FourMomentum and mass define the lab system.
+ */
COMBoost(FourMomentum const& P4projectile, HEPEnergyType const massTarget);
- //! construct a COMBoost given two four-vectors of projectile target
- COMBoost(FourMomentum const& P4projectile, FourMomentum const& P4target);
-
- //! construct a COMBoost to boost into the rest frame given a 3-momentum and mass
+ /**
+ * Construct a COMBoost to boost into the rest frame given a 3-momentum and mass.
+ */
COMBoost(MomentumVector const& momentum, HEPEnergyType const mass);
+ /**
+ * Construct a COMBoost given two four-vectors of projectile target.
+ *
+ * The tow FourMomentum can define an arbitrary system.
+ */
+ COMBoost(FourMomentum const& P4projectile, FourMomentum const& P4target);
+
//! transforms a 4-momentum from lab frame to the center-of-mass frame
template <typename FourVector>
FourVector toCoM(FourVector const& p4) const;
@@ -65,10 +75,10 @@ namespace corsika {
template <typename FourVector>
FourVector fromCoM(FourVector const& p4) const;
- //! returns the rotated coordinate system
+ //! returns the rotated coordinate system: +z is projectile direction
CoordinateSystemPtr getRotatedCS() const;
- //! returns the original coordinate system
+ //! returns the original coordinate system of the projectile (lab)
CoordinateSystemPtr getOriginalCS() const;
protected:
diff --git a/corsika/modules/LongitudinalProfile.hpp b/corsika/modules/LongitudinalProfile.hpp
index 1d84462b0..b5d199fb8 100644
--- a/corsika/modules/LongitudinalProfile.hpp
+++ b/corsika/modules/LongitudinalProfile.hpp
@@ -31,8 +31,7 @@ namespace corsika {
* simulation into a projected grammage range and counts for
* different particle species when they cross dX (default: 10g/cm2)
* boundaries.
- *
- **/
+ */
class LongitudinalProfile : public ContinuousProcess<LongitudinalProfile> {
diff --git a/corsika/modules/QGSJetII.hpp b/corsika/modules/QGSJetII.hpp
index 9540c4ad1..dd7422442 100644
--- a/corsika/modules/QGSJetII.hpp
+++ b/corsika/modules/QGSJetII.hpp
@@ -8,4 +8,23 @@
#pragma once
-#include <corsika/modules/qgsjetII/Interaction.hpp>
+#include <corsika/modules/qgsjetII/InteractionModel.hpp>
+
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+/**
+ * @file QGSJetII.hpp
+ *
+ * Includes all the parts of the QGSJetII model. Defines the InteractionProcess<TModel>
+ * classes needed for the ProcessSequence.
+ */
+
+namespace corsika::qgsjetII {
+ /**
+ * @brief qgsjetII::Interaction is the process for ProcessSequence.
+ *
+ * The qgsjetII::InteractionModel is wrapped as an InteractionProcess here in order
+ * to provide all the functions for ProcessSequence.
+ */
+ class Interaction : public InteractionModel, public InteractionProcess<Interaction> {};
+} // namespace corsika::qgsjetII
diff --git a/corsika/modules/Sibyll.hpp b/corsika/modules/Sibyll.hpp
index 401766122..9737e260d 100644
--- a/corsika/modules/Sibyll.hpp
+++ b/corsika/modules/Sibyll.hpp
@@ -18,7 +18,7 @@
/**
* @file Sibyll.hpp
*
- * Includes all the parts of the Sibyll model. Defines the InteractoinProcess<TModel>
+ * Includes all the parts of the Sibyll model. Defines the InteractionProcess<TModel>
* classes needed for the ProcessSequence.
*/
@@ -26,7 +26,7 @@ namespace corsika::sibyll {
/**
* @brief sibyll::Interaction is the process for ProcessSequence.
*
- * The sibyll::Model is wrapped as an InteractionProcess here in order
+ * The sibyll::InteractionModel is wrapped as an InteractionProcess here in order
* to provide all the functions for ProcessSequence.
*/
class Interaction : public InteractionModel, public InteractionProcess<Interaction> {};
diff --git a/corsika/modules/qgsjetII/Interaction.hpp b/corsika/modules/qgsjetII/Interaction.hpp
deleted file mode 100644
index 0030de2fc..000000000
--- a/corsika/modules/qgsjetII/Interaction.hpp
+++ /dev/null
@@ -1,65 +0,0 @@
-/*
- * (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
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
-#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
-
-#include <boost/filesystem/path.hpp>
-
-#include <qgsjet-II-04.hpp>
-#include <string>
-
-namespace corsika::qgsjetII {
-
- class Interaction : public corsika::InteractionProcess<Interaction> {
-
- public:
- Interaction(boost::filesystem::path dataPath = corsika_data("QGSJetII"));
- ~Interaction();
-
- bool wasInitialized() { return initialized_; }
- unsigned int getMaxTargetMassNumber() const { return maxMassNumber_; }
- bool isValidTarget(corsika::Code TargetId) const {
- return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxMassNumber_);
- }
-
- CrossSectionType getCrossSection(const Code, const Code, const HEPEnergyType,
- const unsigned int Abeam = 0,
- const unsigned int Atarget = 0) const;
-
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle const&) const;
-
- /**
- In this function QGSJETII is called to produce one event. The
- event is copied (and boosted) into the shower lab frame.
- */
-
- template <typename TSecondaryView>
- void doInteraction(TSecondaryView&);
-
- private:
- int count_ = 0;
- bool initialized_ = false;
- QgsjetIIHadronType alternate_ =
- QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
-
- corsika::default_prng_type& rng_ =
- corsika::RNGManager<>::getInstance().getRandomStream("qgsjet");
- static unsigned int constexpr maxMassNumber_ = 208;
- };
-
-} // namespace corsika::qgsjetII
-
-#include <corsika/detail/modules/qgsjetII/Interaction.inl>
diff --git a/corsika/modules/qgsjetII/InteractionModel.hpp b/corsika/modules/qgsjetII/InteractionModel.hpp
new file mode 100644
index 000000000..4ec1bf6a2
--- /dev/null
+++ b/corsika/modules/qgsjetII/InteractionModel.hpp
@@ -0,0 +1,77 @@
+/*
+ * (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
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+#include <qgsjet-II-04.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <boost/filesystem/path.hpp>
+
+namespace corsika::qgsjetII {
+
+ class InteractionModel {
+
+ public:
+ InteractionModel(boost::filesystem::path dataPath = corsika_data("QGSJetII"));
+ ~InteractionModel();
+
+ /**
+ * Throws exception if invalid system is passed.
+ *
+ * @param beamId
+ * @param targetId
+ */
+ void isValid(Code const beamId, Code const targetId) const;
+
+ /**
+ * Return the QGSJETII inelastic/production cross section.
+ *
+ * This cross section must correspond to the process described in doInteraction.
+ * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+ *
+ * @param projectile is the Code of the projectile
+ * @param target is the Code of the target
+ * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+ * @param Aprojectile is the mass number of the projectils, if it is a nucleus
+ * @param Atarget is the mass number of the target, if it is a nucleus
+ *
+ * @return inelastic cross section.
+ */
+ CrossSectionType getCrossSection(Code const projectile, Code const target,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
+
+ /**
+ * In this function QGSJETII is called to produce one event.
+ *
+ * The event is copied (and boosted) into the shower lab frame.
+ */
+ template <typename TSecondaries>
+ void doInteraction(TSecondaries&, Code const projectile, Code const target,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+ private:
+ int count_ = 0;
+ QgsjetIIHadronType alternate_ =
+ QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
+
+ corsika::default_prng_type& rng_ =
+ corsika::RNGManager<>::getInstance().getRandomStream("qgsjet");
+ static size_t constexpr maxMassNumber_ = 208;
+ };
+
+} // namespace corsika::qgsjetII
+
+#include <corsika/detail/modules/qgsjetII/InteractionModel.inl>
diff --git a/corsika/modules/qgsjetII/ParticleConversion.hpp b/corsika/modules/qgsjetII/ParticleConversion.hpp
index 3b4f689b8..0a3ff8c08 100644
--- a/corsika/modules/qgsjetII/ParticleConversion.hpp
+++ b/corsika/modules/qgsjetII/ParticleConversion.hpp
@@ -15,13 +15,13 @@
namespace corsika::qgsjetII {
/**
- These are the possible secondaries produced by QGSJetII
+ * These are the possible secondaries produced by QGSJetII.
*/
enum class QgsjetIICode : int8_t;
using QgsjetIICodeIntType = std::underlying_type<QgsjetIICode>::type;
/**
- These are the possible projectile for which QGSJetII knwos cross section
+ * These are the possible projectile for which QGSJetII knwos cross section.
*/
enum class QgsjetIIXSClass : int8_t {
CannotInteract = 0,
@@ -32,7 +32,7 @@ namespace corsika::qgsjetII {
using QgsjetIIXSClassIntType = std::underlying_type<QgsjetIIXSClass>::type;
/**
- These are the only possible projectile types in QGSJetII
+ * These are the only possible projectile types in QGSJetII.
*/
enum class QgsjetIIHadronType : int8_t {
UndefinedType = 0,
@@ -58,7 +58,7 @@ namespace corsika::qgsjetII {
namespace corsika::qgsjetII {
- QgsjetIICode constexpr convertToQgsjetII(Code pCode) {
+ QgsjetIICode constexpr convertToQgsjetII(Code const pCode) {
return corsika2qgsjetII[static_cast<CodeIntType>(pCode)];
}
diff --git a/corsika/modules/sibyll/InteractionModel.hpp b/corsika/modules/sibyll/InteractionModel.hpp
index 2b4191368..91b4c98a9 100644
--- a/corsika/modules/sibyll/InteractionModel.hpp
+++ b/corsika/modules/sibyll/InteractionModel.hpp
@@ -48,7 +48,7 @@ namespace corsika::sibyll {
HEPEnergyType const sqrtSnn) const;
/**
- * @brief returns inelastic AND elastic cross sections.
+ * Returns inelastic AND elastic cross sections.
*
* These cross sections must correspond to the process described in doInteraction
* AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
@@ -68,7 +68,7 @@ namespace corsika::sibyll {
FourMomentum const& targetP4) const;
/**
- * @brief returns inelastic (production) cross section.
+ * Returns inelastic (production) cross section.
*
* This cross section must correspond to the process described in doInteraction.
* Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 496deb3eb..5d6ec770f 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -4,7 +4,7 @@ set (test_modules_sources
testTracking.cpp
## testExecTime.cpp --> needs to be fixed, see #326
testObservationPlane.cpp
- #testQGSJetII.cpp
+ testQGSJetII.cpp
#testPythia8.cpp
#testUrQMD.cpp
#testCONEX.cpp
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 84f1b4ebe..2a8cd57a8 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -6,8 +6,7 @@
* the license.
*/
-#include <corsika/modules/qgsjetII/Interaction.hpp>
-#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+#include <corsika/modules/QGSJetII.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
@@ -49,6 +48,7 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("QgsjetII", "[processes]") {
logging::set_level(logging::level::info);
+ RNGManager<>::getInstance().registerRandomStream("qgsjet");
SECTION("Corsika -> QgsjetII") {
CHECK(corsika::qgsjetII::convertToQgsjetII(PiMinus::code) ==
@@ -91,6 +91,14 @@ TEST_CASE("QgsjetII", "[processes]") {
CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::PiMinus) ==
corsika::qgsjetII::QgsjetIIXSClass::LightMesons);
}
+
+ SECTION("valid") {
+
+ corsika::qgsjetII::InteractionModel model;
+
+ CHECK_THROWS(model.isValid(Code::Electron, Code::Proton));
+ CHECK_THROWS(model.isValid(Code::Proton, Code::Electron));
+ }
}
#include <corsika/framework/geometry/Point.hpp>
@@ -115,25 +123,23 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
logging::set_level(logging::level::info);
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+ auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
[[maybe_unused]] auto const& node_dummy = nodePtr;
- RNGManager<>::getInstance().registerRandomStream("qgsjet");
-
SECTION("InteractionInterface") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Proton, 110_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
- auto particle = stackPtr->first();
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
- corsika::qgsjetII::Interaction model;
- model.doInteraction(view);
- [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
-
- CHECK(length / (1_g / square(1_cm)) == Approx(93.04).margin(0.1));
+ corsika::qgsjetII::InteractionModel model;
+ model.doInteraction(view, Code::Proton, Code::Oxygen,
+ {sqrt(static_pow<2>(110_GeV) + static_pow<2>(Proton::mass)),
+ MomentumVector{cs, 110_GeV, 0_GeV, 0_GeV}},
+ {Oxygen::mass, MomentumVector{cs, {0_eV, 0_eV, 0_eV}}});
/* **********************************
As it turned out already two times (#291 and #307) that the detailed output of
@@ -152,24 +158,25 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
SECTION("InteractionInterface Nuclei") {
+ HEPEnergyType const P0 = 20100_GeV;
+ MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
+ Code const pid = get_nucleus_code(60, 30);
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- get_nucleus_code(60, 30), 20100_GeV,
- (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ pid, P0, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
- auto particle = stackPtr->first();
- auto projectile = secViewPtr->getProjectile();
- auto const projectileMomentum = projectile.getMomentum();
- corsika::qgsjetII::Interaction model;
- model.doInteraction(view); // this also should produce some fragments
+ HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid)));
+ FourMomentum const projectileP4(Elab, plab);
+ FourMomentum const targetP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+ view.clear();
+
+ corsika::qgsjetII::InteractionModel model;
+ model.doInteraction(view, pid, Code::Oxygen, projectileP4,
+ targetP4); // this also should produce some fragments
CHECK(view.getSize() == Approx(300).margin(150)); // this is not physics validation
int countFragments = 0;
for (auto const& sec : view) { countFragments += (is_nucleus(sec.getPID())); }
CHECK(countFragments == Approx(4).margin(2)); // this is not physics validation
- [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
-
- CHECK(length / (1_g / square(1_cm)) ==
- Approx(12).margin(2)); // this is not physics validation
}
SECTION("Heavy nuclei") {
@@ -178,38 +185,33 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
get_nucleus_code(1000, 1000), 1100_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
- auto particle = stackPtr->first();
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
- corsika::qgsjetII::Interaction model;
+ corsika::qgsjetII::InteractionModel model;
+ FourMomentum const aP4(100_GeV, {cs, 99_GeV, 0_GeV, 0_GeV});
+ FourMomentum const bP4(1_TeV, {cs, 0.9_TeV, 0_GeV, 0_GeV});
+
+ CHECK_THROWS(model.getCrossSection(get_nucleus_code(10, 5),
+ get_nucleus_code(1000, 500), aP4, bP4));
+ CHECK_THROWS(model.getCrossSection(Code::Nucleus, Code::Nucleus, aP4, bP4));
CHECK_THROWS(
- model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 10., 1000.));
- CHECK_THROWS(
- model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 1000., 10.));
- CHECK_THROWS(model.doInteraction(view));
- CHECK_THROWS(model.getInteractionLength(particle));
+ model.doInteraction(view, get_nucleus_code(1000, 500), Code::Oxygen, aP4, bP4));
}
SECTION("Allowed Particles") {
- { // electron
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Electron, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
- [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
- auto particle = stackPtr->first();
- corsika::qgsjetII::Interaction model;
- GrammageType const length = model.getInteractionLength(particle);
- CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
- }
+
{ // pi0 is internally converted into pi+/pi-
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Pi0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
- corsika::qgsjetII::Interaction model;
- model.doInteraction(view);
+ corsika::qgsjetII::InteractionModel model;
+ model.doInteraction(view, Code::Pi0, Code::Oxygen,
+ {sqrt(static_pow<2>(1_TeV) + static_pow<2>(Pi0::mass)),
+ MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
+ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(10).margin(8)); // this is not physics validation
}
{ // rho0 is internally converted into pi-/pi+
@@ -217,8 +219,11 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
Code::Rho0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
- corsika::qgsjetII::Interaction model;
- model.doInteraction(view);
+ corsika::qgsjetII::InteractionModel model;
+ model.doInteraction(view, Code::Rho0, Code::Oxygen,
+ {sqrt(static_pow<2>(1_TeV) + static_pow<2>(Rho0::mass)),
+ MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
+ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(25).margin(20)); // this is not physics validation
}
{ // Lambda is internally converted into neutron
@@ -227,8 +232,11 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
*csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
- corsika::qgsjetII::Interaction model;
- model.doInteraction(view);
+ corsika::qgsjetII::InteractionModel model;
+ model.doInteraction(view, Code::Lambda0, Code::Oxygen,
+ {sqrt(static_pow<2>(100_GeV) + static_pow<2>(Lambda0::mass)),
+ MomentumVector{cs, 100_GeV, 0_GeV, 0_GeV}},
+ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(25).margin(20)); // this is not physics validation
}
{ // AntiLambda is internally converted into anti neutron
@@ -237,8 +245,11 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
*csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
- corsika::qgsjetII::Interaction model;
- model.doInteraction(view);
+ corsika::qgsjetII::InteractionModel model;
+ model.doInteraction(view, Code::Lambda0Bar, Code::Oxygen,
+ {sqrt(static_pow<2>(1_TeV) + static_pow<2>(Lambda0Bar::mass)),
+ MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
+ {Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(70).margin(67)); // this is not physics validation
}
}
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index a08e33d8c..718a43601 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -88,7 +88,6 @@ TEST_CASE("Sibyll", "modules") {
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-
#include <corsika/framework/core/ParticleProperties.hpp>
#include <SetupTestEnvironment.hpp>
@@ -186,9 +185,9 @@ TEST_CASE("SibyllInterface", "modules") {
model.setVerbose(true);
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(Proton::mass));
FourMomentum const projectileP4(Elab, plab);
- FourMomentum const nucleonP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+ FourMomentum const nucleusP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
view.clear();
- model.doInteraction(view, Code::Proton, Code::Oxygen, projectileP4, nucleonP4);
+ model.doInteraction(view, Code::Proton, Code::Oxygen, projectileP4, nucleusP4);
auto const pSum = sumMomentum(view, cs);
/*
@@ -255,7 +254,7 @@ TEST_CASE("SibyllInterface", "modules") {
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] CrossSectionType const cx =
- model.getCrossSection(Code::Proton, Code::Oxygen, projectileP4, nucleonP4);
+ model.getCrossSection(Code::Proton, Code::Oxygen, projectileP4, nucleusP4);
CHECK(cx / 1_mb == Approx(300).margin(1));
// CHECK(view.getEntries() == 9); //! \todo: this was 20 before refactory-2020: check
// "also sibyll not stable wrt. to compiler
--
GitLab