Merge branch '238-add-qgsjetii' into 'master'

Resolve "Add qgsjetII"

Closes #238

See merge request AirShowerPhysics/corsika!167

.gitlab-ci.yml

@@ -19,7 +19,10 @@ config-u-18_04:
   stage: config
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
+    - git clone https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika-data.git    
     - mkdir build
     - cd build
     - cmake .. -DCMAKE_BUILD_TYPE=Debug -DWITH_PYTHIA=ON
@@ -27,6 +30,7 @@ config-u-18_04:
     expire_in: 1 day
     paths:
        - build
+       - corsika-data
 
 # job/stage to just prepare cmake
 config-clang-8:
@@ -34,7 +38,10 @@ config-clang-8:
   stage: config
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_BUILDS_DIR}/AirShowerPhysics/corsika/corsika-data/"
   script:
+    - git clone https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika-data.git
     - mkdir build
     - cd build
     - cmake .. -DCMAKE_BUILD_TYPE=Debug -DWITH_PYTHIA=ON
@@ -42,6 +49,7 @@ config-clang-8:
     expire_in: 1 day
     paths:
        - build
+       - corsika-data
 
 # normal pipeline for each commit
 build-test-u-18_04: 
@@ -51,6 +59,8 @@ build-test-u-18_04:
   stage: build_test
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
     - cd build
     - cmake --build . -- -j4
@@ -73,6 +83,8 @@ build-test-clang-8:
   stage: build_test
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
     - cd build
     - cmake --build . -- -j4
@@ -95,6 +107,8 @@ release-u-18_04:
   stage: optional
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
     - cd build
     - cmake .. -DCMAKE_BUILD_TYPE=Release
@@ -119,6 +133,8 @@ release-clang-8:
   stage: optional
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
     - cd build
     - cmake .. -DCMAKE_BUILD_TYPE=Release
@@ -143,6 +159,8 @@ release-clang-8:
   stage: optional
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
     - cd build
     - cmake .. -DCMAKE_BUILD_TYPE=Coverage
@@ -176,6 +194,8 @@ documentation:
   stage: optional
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
     - cd build
     - cmake --build . --target doxygen -- -j4
@@ -196,6 +216,8 @@ sanity:
   stage: optional
   tags:
     - corsika
+  variables:
+    CORSIKA_DATA: "${CI_PROJECT_DIR}/corsika-data/"
   script:
     - cd build
     - cmake .. -DWITH_CORSIKA_SANITIZERS_ENABLED=ON

Processes/CMakeLists.txt

@@ -4,6 +4,7 @@ add_subdirectory (NullModel)
 add_subdirectory (TrackingLine) 
 # hadron interaction models
 add_subdirectory (Sibyll)
+add_subdirectory (QGSJetII)
 if (PYTHIA8_FOUND)
   add_subdirectory (Pythia)
 endif (PYTHIA8_FOUND)

Processes/QGSJetII/CMakeLists.txt

+set(Python_ADDITIONAL_VERSIONS 3)
+find_package(PythonInterp 3 REQUIRED)
+
+if (NOT DEFINED ENV{CORSIKA_DATA})
+  message (WARNING "WARNING: download corsika-data repository from gitlab for needed data tables, and set CORSIKA_DATA path!")
+endif ()
+
+add_custom_command (
+  OUTPUT  ${PROJECT_BINARY_DIR}/Processes/QGSJetII/Generated.inc
+  COMMAND ${PROJECT_SOURCE_DIR}/Processes/QGSJetII/code_generator.py 
+          ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
+          ${PROJECT_SOURCE_DIR}/Processes/QGSJetII/qgsjet-II-04-codes.dat
+  DEPENDS code_generator.py
+          qgsjet-II-04-codes.dat
+          ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
+  WORKING_DIRECTORY
+          ${PROJECT_BINARY_DIR}/Processes/QGSJetII/
+  COMMENT "Generate conversion tables for particle codes QGSJetII <-> CORSIKA"
+  VERBATIM
+  )
+  
+set (
+  MODEL_SOURCES
+  ParticleConversion.cc
+  Interaction.cc
+  qgsjet-II-04.f
+  qgsjet-II-04.cc
+  )
+
+set (
+  MODEL_HEADERS
+  ParticleConversion.h
+  qgsjet-II-04.h
+  QGSJetIIStack.h
+  QGSJetIIFragmentsStack.h
+  Interaction.h
+  ${PROJECT_BINARY_DIR}/Processes/QGSJetII/Generated.inc
+  )
+
+set (
+  MODEL_NAMESPACE
+  corsika/process/qgsjetII
+  )
+
+add_library (ProcessQGSJetII STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessQGSJetII ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+# ....................................................
+# since Generated.inc is an automatically produced file in the build directory,
+# create a symbolic link into the source tree, so that it can be found and edited more easily
+# this is not needed for the build to succeed! .......
+add_custom_command (
+  OUTPUT  ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
+  COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/process/qgsjetII/Generated.inc ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
+  COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc"
+  )
+add_custom_target (SourceDirLinkQgsII DEPENDS ${PROJECT_BINARY_DIR}/Processes/QGSJetII/Generated.inc)
+add_dependencies (ProcessQGSJetII SourceDirLinkQgsII)
+# .....................................................
+
+
+
+set_target_properties (
+  ProcessQGSJetII
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+  )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+  ProcessQGSJetII
+  CORSIKAprocesssequence
+  CORSIKAparticles
+  CORSIKAutilities
+  CORSIKAunits
+  CORSIKAthirdparty
+  CORSIKAgeometry
+  CORSIKAenvironment
+  )
+
+target_include_directories (
+  ProcessQGSJetII 
+  INTERFACE 
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS ProcessQGSJetII
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+  )
+
+# install(
+#    FILES qgsdat-II-04 sectnu-II-04 DESTINATION etc/data
+# )
+
+
+# --------------------
+# code unit testing
+CORSIKA_ADD_TEST(testQGSJetII
+  SOURCES
+  testQGSJetII.cc
+  ${MODEL_HEADERS}
+)
+
+target_link_libraries (
+  testQGSJetII
+  ProcessQGSJetII
+  CORSIKAsetup
+  CORSIKArandom
+  CORSIKAgeometry
+  CORSIKAunits
+  CORSIKAtesting
+  )
+
+# also provide QGSJetII large data tables for testing in build tree (links)
+# -> changed that to use environment variable "${CORSIKA_DATA}/QGSJetII"
+# add_custom_command (
+#  OUTPUT  ${CMAKE_CURRENT_BINARY_DIR}/qgsdat-II-04 ${CMAKE_CURRENT_BINARY_DIR}/sectnu-II-04
+#  COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR}/qgsdat-II-04 ${CMAKE_CURRENT_BINARY_DIR}/qgsdat-II-04
+#  COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR}/sectnu-II-04 ${CMAKE_CURRENT_BINARY_DIR}/sectnu-II-04
+#  COMMENT "Generate link in source-dir: qgsjet-II-04"
+#  )
+
+# add_custom_target (QGSJetDataLink DEPENDS
+#  ${CMAKE_CURRENT_BINARY_DIR}/qgsdat-II-04
+#  ${CMAKE_CURRENT_BINARY_DIR}/sectnu-II-04
+#  )
+
+# add_dependencies (testQGSJetII QGSJetDataLink)

Processes/QGSJetII/Interaction.cc

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/process/qgsjetII/Interaction.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/QuantityVector.h>
+#include <corsika/geometry/FourVector.h>
+#include <corsika/process/qgsjetII/ParticleConversion.h>
+#include <corsika/process/qgsjetII/QGSJetIIFragmentsStack.h>
+#include <corsika/process/qgsjetII/QGSJetIIStack.h>
+#include <corsika/process/qgsjetII/qgsjet-II-04.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/utl/COMBoost.h>
+
+#include <sstream>
+#include <string>
+#include <tuple>
+
+using std::cout;
+using std::endl;
+using std::ostringstream;
+using std::string;
+using std::tuple;
+
+using namespace corsika;
+using namespace corsika::setup;
+using SetupParticle = setup::Stack::StackIterator;
+using SetupProjectile = setup::StackView::StackIterator;
+using Track = Trajectory;
+
+namespace corsika::process::qgsjetII {
+
+  Interaction::Interaction(const string& dataPath)
+      : data_path_(dataPath) {
+    if (dataPath == "") {
+      if (std::getenv("CORSIKA_DATA")) {
+        data_path_ = string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
+        cout << "Searching for QGSJetII data tables in " << data_path_ << endl;
+      }
+    }
+  }
+
+  Interaction::~Interaction() { cout << "QgsjetII::Interaction n=" << count_ << endl; }
+
+  void Interaction::Init() {
+
+    using random::RNGManager;
+
+    // initialize QgsjetII
+    if (!initialized_) {
+      qgset_();
+      datadir DIR(data_path_);
+      qgaini_(DIR.data);
+      initialized_ = true;
+    }
+  }
+
+  units::si::CrossSectionType Interaction::GetCrossSection(
+      const particles::Code beamId, const particles::Code targetId,
+      const units::si::HEPEnergyType Elab, const unsigned int Abeam,
+      const unsigned int targetA) const {
+    using namespace units::si;
+    double sigProd = std::numeric_limits<double>::infinity();
+
+    if (process::qgsjetII::CanInteract(beamId)) {
+
+      const int iBeam = process::qgsjetII::GetQgsjetIIXSCode(beamId);
+      int iTarget = 1;
+      if (particles::IsNucleus(targetId)) {
+        iTarget = targetA;
+        if (iTarget > maxMassNumber_ || iTarget <= 0) {
+          std::ostringstream txt;
+          txt << "QgsjetII target outside range. iTarget=" << iTarget;
+          throw std::runtime_error(txt.str().c_str());
+        }
+      }
+      int iProjectile = 1;
+      if (particles::IsNucleus(beamId)) {
+        iProjectile = Abeam;
+        if (iProjectile > maxMassNumber_ || iProjectile <= 0)
+          throw std::runtime_error("QgsjetII target outside range. ");
+      }
+
+      cout << "QgsjetII::GetCrossSection Elab=" << Elab << " iBeam=" << iBeam
+           << " iProjectile=" << iProjectile << " iTarget=" << iTarget << endl;
+      sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
+      cout << "QgsjetII::GetCrossSection sigProd=" << sigProd << endl;
+    }
+
+    return sigProd * 1_mb;
+  }
+
+  template <>
+  units::si::GrammageType Interaction::GetInteractionLength(
+      SetupParticle const& vP) const {
+
+    using namespace units;
+    using namespace units::si;
+    using namespace geometry;
+
+    // coordinate system, get global frame of reference
+    CoordinateSystem& rootCS =
+        RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+    const particles::Code corsikaBeamId = vP.GetPID();
+
+    // beam particles for qgsjetII : 1, 2, 3 for p, pi, k
+    // read from cross section code table
+    const bool kInteraction = process::qgsjetII::CanInteract(corsikaBeamId);
+
+    // FOR NOW: assume target is at rest
+    MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
+
+    // total momentum and energy
+    HEPEnergyType Elab = vP.GetEnergy();
+
+    cout << "Interaction: LambdaInt: \n"
+         << " input energy: " << vP.GetEnergy() / 1_GeV << endl
+         << " beam can interact:" << kInteraction << endl
+         << " beam pid:" << vP.GetPID() << endl;
+
+    if (kInteraction) {
+
+      int Abeam = 0;
+      if (particles::IsNucleus(vP.GetPID())) Abeam = vP.GetNuclearA();
+
+      // 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 = vP.GetNode();
+      const auto& mediumComposition =
+          currentNode->GetModelProperties().GetNuclearComposition();
+
+      si::CrossSectionType weightedProdCrossSection = mediumComposition.WeightedSum(
+          [=](particles::Code targetID) -> si::CrossSectionType {
+            int targetA = 0;
+            if (corsika::particles::IsNucleus(targetID))
+              targetA = particles::GetNucleusA(targetID);
+            return GetCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA);
+          });
+
+      cout << "Interaction: "
+           << "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
+           << endl;
+
+      // calculate interaction length in medium
+      GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
+                                      units::constants::u / weightedProdCrossSection;
+      cout << "Interaction: "
+           << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
+           << endl;
+
+      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 <>
+  process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
+
+    using namespace units;
+    using namespace utl;
+    using namespace units::si;
+    using namespace geometry;
+
+    const auto corsikaBeamId = vP.GetPID();
+    cout << "ProcessQgsjetII: "
+         << "DoInteraction: " << corsikaBeamId << " interaction? "
+         << process::qgsjetII::CanInteract(corsikaBeamId) << endl;
+
+    if (process::qgsjetII::CanInteract(corsikaBeamId)) {
+
+      const CoordinateSystem& rootCS =
+          RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+      // position and time of interaction, not used in QgsjetII
+      Point pOrig = vP.GetPosition();
+      TimeType tOrig = vP.GetTime();
+
+      // define target
+      // for QgsjetII is always a single nucleon
+      // FOR NOW: target is always at rest
+      const auto targetEnergyLab = 0_GeV + constants::nucleonMass;
+      const auto targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+      const FourVector PtargLab(targetEnergyLab, targetMomentumLab);
+
+      // define projectile
+      HEPEnergyType const projectileEnergyLab = vP.GetEnergy();
+      auto const projectileMomentumLab = vP.GetMomentum();
+
+      int beamA = 0;
+      if (particles::IsNucleus(corsikaBeamId)) beamA = vP.GetNuclearA();
+
+      cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << endl
+           << "Interaction: pbeam lab: " << projectileMomentumLab.GetComponents() / 1_GeV
+           << endl;
+      cout << "Interaction: etarget lab: " << targetEnergyLab / 1_GeV << endl
+           << "Interaction: ptarget lab: " << targetMomentumLab.GetComponents() / 1_GeV << endl;
+
+      cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
+      cout << "Interaction: time: " << tOrig << endl;
+
+      // sample target mass number
+      auto const* currentNode = vP.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<si::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 (corsika::particles::IsNucleus(targetId))
+          targetA = particles::GetNucleusA(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, fRNG);
+      cout << "Interaction: target selected: " << targetCode << endl;
+
+      int targetQgsCode = -1;
+      if (particles::IsNucleus(targetCode))
+        targetQgsCode = particles::GetNucleusA(targetCode);
+      if (targetCode == particles::Proton::GetCode()) targetQgsCode = 1;
+      cout << "Interaction: target qgsjetII code/A: " << targetQgsCode << endl;
+      if (targetQgsCode > maxMassNumber_ || targetQgsCode < 1)
+        throw std::runtime_error("QgsjetII target outside range.");
+
+      int projQgsCode = 1;
+      if (particles::IsNucleus(corsikaBeamId)) projQgsCode = vP.GetNuclearA();
+      cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " "
+           << corsikaBeamId << endl;
+      if (projQgsCode > maxMassNumber_ || projQgsCode < 1)
+        throw std::runtime_error("QgsjetII target outside range.");
+
+      // beam id for qgsjetII
+      int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
+      if (corsikaBeamId != particles::Code::Nucleus) {
+        kBeam = process::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
+      }
+
+      cout << "Interaction: "
+           << " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << endl;
+      count_++;
+      qgini_(projectileEnergyLab / 1_GeV, kBeam, projQgsCode, targetQgsCode);
+      // this is from CRMC, is this REALLY needed ???
+      qgini_(projectileEnergyLab / 1_GeV, kBeam, projQgsCode, targetQgsCode);
+      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();
+      geometry::CoordinateSystem const zAxisFrame =
+	originalCS.RotateToZ(projectileMomentumLab);
+      
+      // fragments
+      QGSJetIIFragmentsStack qfs;
+      for (auto& fragm : qfs) {
+        particles::Code idFragm = particles::Code::Nucleus;
+        int A = fragm.GetFragmentSize();
+        int Z = 0;
+        switch (A) {
+          case 1: { // proton/neutron
+            idFragm = particles::Code::Proton;
+
+	    auto momentum = geometry::Vector(
+					     zAxisFrame,
+					     corsika::geometry::QuantityVector<hepmomentum_d>{0.0_GeV, 0.0_GeV,
+						 sqrt((projectileEnergyLab + particles::Proton::GetMass()) *
+						      (projectileEnergyLab - particles::Proton::GetMass()))});
+	    
+	    auto const energy = sqrt(momentum.squaredNorm() + square(particles::GetMass(idFragm)));	    
+	    momentum.rebase(originalCS); // transform back into standard lab frame
+	    std::cout << "secondary fragment> id=" << idFragm << " p=" << momentum.GetComponents() << std::endl;
+            auto pnew = vP.AddSecondary(
+                tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+                      geometry::Point, units::si::TimeType>{
+		  idFragm, energy, momentum, pOrig, tOrig});
+            Plab_final += pnew.GetMomentum();
+            Elab_final += pnew.GetEnergy();
+          } break;
+          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);
+          }
+        }
+
+        if (idFragm == particles::Code::Nucleus) {
+	    auto momentum = geometry::Vector(
+					     zAxisFrame,
+					     geometry::QuantityVector<hepmomentum_d>{0.0_GeV, 0.0_GeV,
+						 sqrt((projectileEnergyLab + constants::nucleonMass * A) *
+						      (projectileEnergyLab - constants::nucleonMass * A))});
+	    
+	    auto const energy = sqrt(momentum.squaredNorm() + square(constants::nucleonMass*A));	    
+	    momentum.rebase(originalCS); // transform back into standard lab frame
+	    std::cout << "secondary fragment> id=" << idFragm << " p=" << momentum.GetComponents() << " A=" << A << " Z=" << Z << std::endl;
+            auto pnew = vP.AddSecondary(
+                tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+		geometry::Point, units::si::TimeType, unsigned short, unsigned short>{
+		  idFragm, energy, momentum, pOrig, tOrig, A, Z});
+            Plab_final += pnew.GetMomentum();
+            Elab_final += pnew.GetEnergy();
+        }
+      }
+
+      // secondaries
+      QGSJetIIStack qs;
+      for (auto& psec : qs) {
+
+        auto momentum = psec.GetMomentum(zAxisFrame);
+        auto const energy = psec.GetEnergy();
+
+	momentum.rebase(originalCS); // transform back into standard lab frame
+	std::cout << "secondary fragment> id=" << process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()) << " p=" << momentum.GetComponents() << std::endl;
+	auto pnew = vP.AddSecondary(
+				    tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
+				    geometry::Point, units::si::TimeType>{
+		  process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()), energy, momentum, pOrig, tOrig});
+	Plab_final += pnew.GetMomentum();
+	Elab_final += pnew.GetEnergy();
+      }
+      cout << "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/ << endl
+           << "Elab_final=" << Elab_final / 1_GeV
+           << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
+           << ", N_wounded,targ="
+           << QGSJetIIFragmentsStackData::GetWoundedNucleonsTarget()
+           << ", N_wounded,proj="
+           << QGSJetIIFragmentsStackData::GetWoundedNucleonsProjectile()
+           << ", N_fragm,proj=" << qfs.GetSize() << endl;
+    }
+    return process::EProcessReturn::eOk;
+  }
+
+} // namespace corsika::process::qgsjetII

Processes/QGSJetII/Interaction.h

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _corsika_process_qgsjetII_interaction_h_
+#define _corsika_process_qgsjetII_interaction_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <string>
+
+namespace corsika::process::qgsjetII {
+
+  class Interaction : public corsika::process::InteractionProcess<Interaction> {
+
+    std::string data_path_;
+    int count_ = 0;
+    bool initialized_ = false;
+
+  public:
+    Interaction(const std::string& dataPath = "");
+    ~Interaction();
+
+    void Init();
+
+    bool WasInitialized() { return initialized_; }
+    int GetMaxTargetMassNumber() const { return maxMassNumber_; }
+    bool IsValidTarget(corsika::particles::Code TargetId) const {
+      return (corsika::particles::GetNucleusA(TargetId) < maxMassNumber_) &&
+             corsika::particles::IsNucleus(TargetId);
+    }
+
+    corsika::units::si::CrossSectionType GetCrossSection(
+        const corsika::particles::Code, const corsika::particles::Code,
+        const corsika::units::si::HEPEnergyType, const unsigned int Abeam = 0,
+        const unsigned int Atarget = 0) const;
+
+    template <typename TParticle>
+    corsika::units::si::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 TProjectile>
+    corsika::process::EProcessReturn DoInteraction(TProjectile&);
+
+  private:
+    corsika::random::RNG& fRNG =
+        corsika::random::RNGManager::GetInstance().GetRandomStream("qgran");
+    const int maxMassNumber_ = 208;
+  };
+
+} // namespace corsika::process::qgsjetII
+
+#endif

Processes/QGSJetII/ParticleConversion.cc

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/particles/ParticleProperties.h>
+#include <corsika/process/qgsjetII/ParticleConversion.h>
+
+using namespace corsika::process::qgsjetII;

Processes/QGSJetII/ParticleConversion.h

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_processes_qgsjetII_particles_h_
+#define _include_processes_qgsjetII_particles_h_
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <string>
+
+namespace corsika::process::qgsjetII {
+
+  enum class QgsjetIICode : int8_t;
+  using QgsjetIICodeIntType = std::underlying_type<QgsjetIICode>::type;
+
+#include <corsika/process/qgsjetII/Generated.inc>
+
+  QgsjetIICode constexpr ConvertToQgsjetII(corsika::particles::Code pCode) {
+    return static_cast<QgsjetIICode>(
+        corsika2qgsjetII[static_cast<corsika::particles::CodeIntType>(pCode)]);
+  }
+
+  corsika::particles::Code constexpr ConvertFromQgsjetII(QgsjetIICode pCode) {
+    auto const pCodeInt = static_cast<QgsjetIICodeIntType>(pCode);
+    auto const corsikaCode = qgsjetII2corsika[pCodeInt - minQgsjetII];
+    if (corsikaCode == corsika::particles::Code::Unknown) {
+      throw std::runtime_error(std::string("QGSJETII/CORSIKA conversion of pCodeInt=")
+                                   .append(std::to_string(pCodeInt))
+                                   .append(" impossible"));
+    }
+    return corsikaCode;
+  }
+
+  int constexpr ConvertToQgsjetIIRaw(corsika::particles::Code pCode) {
+    return static_cast<int>(ConvertToQgsjetII(pCode));
+  }
+
+  int constexpr GetQgsjetIIXSCode(corsika::particles::Code pCode) {
+    if (pCode == corsika::particles::Code::Nucleus) return 2;
+    return corsika2qgsjetIIXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
+  }
+
+  bool constexpr CanInteract(corsika::particles::Code pCode) {
+    return (GetQgsjetIIXSCode(pCode) > 0) && (ConvertToQgsjetIIRaw(pCode) <= 5);
+  }
+
+} // namespace corsika::process::qgsjetII
+
+#endif

Processes/QGSJetII/QGSJetIIFragmentsStack.h

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_qgsjetIIfragmentsstack_h_
+#define _include_qgsjetIIfragmentsstack_h_
+
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/process/qgsjetII/ParticleConversion.h>
+#include <corsika/process/qgsjetII/qgsjet-II-04.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process::qgsjetII {
+
+  class QGSJetIIFragmentsStackData {
+
+  public:
+    void Init();
+    void Dump() const {}
+
+    void Clear() {
+      qgarr13_.nsf = 0;
+      qgarr55_.nwt = 0;
+    }
+    unsigned int GetSize() const { return qgarr13_.nsf; }
+    unsigned int GetCapacity() const { return iapmax; }
+
+    static unsigned int GetWoundedNucleonsTarget() { return qgarr55_.nwt; }
+    static unsigned int GetWoundedNucleonsProjectile() { return qgarr55_.nwp; }
+
+    int GetFragmentSize(const unsigned int i) const { return qgarr13_.iaf[i]; }
+    void SetFragmentSize(const unsigned int i, const int v) { qgarr13_.iaf[i] = v; }
+
+    void Copy(const unsigned int i1, const unsigned int i2) {
+      qgarr13_.iaf[i2] = qgarr13_.iaf[i1];
+    }
+
+    void Swap(const unsigned int i1, const unsigned int i2) {
+      std::swap(qgarr13_.iaf[i1], qgarr13_.iaf[i2]);
+    }
+
+    void IncrementSize() { qgarr13_.nsf++; }
+    void DecrementSize() {
+      if (qgarr13_.nsf > 0) { qgarr13_.nsf--; }
+    }
+  };
+
+  template <typename StackIteratorInterface>
+  class FragmentsInterface : public corsika::stack::ParticleBase<StackIteratorInterface> {
+
+    using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+    using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+  public:
+    void SetParticleData(const int vSize) { SetFragmentSize(vSize); }
+
+    void SetParticleData(FragmentsInterface<StackIteratorInterface>& /*parent*/,
+                         const int vSize) {
+      SetFragmentSize(vSize);
+    }
+
+    void SetFragmentSize(const int v) { GetStackData().SetFragmentSize(GetIndex(), v); }
+
+    double GetFragmentSize() const { return GetStackData().GetFragmentSize(GetIndex()); }
+  };
+
+  typedef corsika::stack::Stack<QGSJetIIFragmentsStackData, FragmentsInterface>
+      QGSJetIIFragmentsStack;
+
+} // end namespace corsika::process::qgsjetII
+
+#endif

Processes/QGSJetII/QGSJetIIStack.h

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_qgsjetIIstack_h_
+#define _include_qgsjetIIstack_h_
+
+#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/process/qgsjetII/ParticleConversion.h>
+#include <corsika/process/qgsjetII/qgsjet-II-04.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process::qgsjetII {
+
+  typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
+
+  class QGSJetIIStackData {
+
+  public:
+    void Init();
+    void Dump() const {}
+
+    void Clear() {
+      qgarr12_.nsp = 0;
+      qgarr13_.nsf = 0;
+      qgarr55_.nwt = 0;
+    }
+    unsigned int GetSize() const { return qgarr12_.nsp; }
+    unsigned int GetCapacity() const { return nptmax; }
+
+    void SetId(const unsigned int i, const int v) { qgarr14_.ich[i] = v; }
+    void SetEnergy(const unsigned int i, const corsika::units::si::HEPEnergyType v) {
+      using namespace corsika::units::si;
+      qgarr14_.esp[i][0] = v / 1_GeV;
+    }
+
+    void SetMomentum(const unsigned int i, const MomentumVector& v) {
+      using namespace corsika::units::si;
+      auto tmp = v.GetComponents();
+      qgarr14_.esp[i][2] = tmp[0] / 1_GeV;
+      qgarr14_.esp[i][3] = tmp[1] / 1_GeV;
+      qgarr14_.esp[i][1] = tmp[2] / 1_GeV;
+    }
+
+    int GetId(const unsigned int i) const { return qgarr14_.ich[i]; }
+    corsika::units::si::HEPEnergyType GetEnergy(const int i) const {
+      using namespace corsika::units::si;
+      return qgarr14_.esp[i][0] * 1_GeV;
+    }
+    MomentumVector GetMomentum(const unsigned int i, const corsika::geometry::CoordinateSystem& CS) const {
+      using namespace corsika::units::si;
+      geometry::QuantityVector<hepmomentum_d> components = {qgarr14_.esp[i][2] * 1_GeV,
+                                                  qgarr14_.esp[i][3] * 1_GeV,
+                                                  qgarr14_.esp[i][1] * 1_GeV};
+      return MomentumVector(CS, components);
+    }
+
+    void Copy(const unsigned int i1, const unsigned int i2) {
+      qgarr14_.ich[i2] = qgarr14_.ich[i1];
+      for (unsigned int i = 0; i < 4; ++i) qgarr14_.esp[i2][i] = qgarr14_.esp[i1][i];
+    }
+
+    void Swap(const unsigned int i1, const unsigned int i2) {
+      std::swap(qgarr14_.ich[i1], qgarr14_.ich[i2]);
+      for (unsigned int i = 0; i < 4; ++i)
+        std::swap(qgarr14_.esp[i1][i], qgarr14_.esp[i2][i]);
+    }
+
+    void IncrementSize() { qgarr12_.nsp++; }
+    void DecrementSize() {
+      if (qgarr12_.nsp > 0) { qgarr12_.nsp--; }
+    }
+  };
+
+  template <typename StackIteratorInterface>
+  class ParticleInterface : public corsika::stack::ParticleBase<StackIteratorInterface> {
+
+    using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+    using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+  public:
+    void SetParticleData(const int vID, 
+                         const corsika::units::si::HEPEnergyType vE,
+                         const MomentumVector& vP,
+                         const corsika::units::si::HEPMassType vM) {
+      SetPID(vID);
+      SetEnergy(vE);
+      SetMomentum(vP);
+    }
+
+    void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
+                         const int vID, 
+                         const corsika::units::si::HEPEnergyType vE,
+                         const MomentumVector& vP,
+                         const corsika::units::si::HEPMassType vM) {
+      SetPID(vID);
+      SetEnergy(vE);
+      SetMomentum(vP);
+    }
+
+    void SetEnergy(const corsika::units::si::HEPEnergyType v) {
+      GetStackData().SetEnergy(GetIndex(), v);
+    }
+
+    corsika::units::si::HEPEnergyType GetEnergy() const {
+      return GetStackData().GetEnergy(GetIndex());
+    }
+
+    void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
+
+    corsika::process::qgsjetII::QgsjetIICode GetPID() const {
+      return static_cast<corsika::process::qgsjetII::QgsjetIICode>(
+          GetStackData().GetId(GetIndex()));
+    }
+
+    MomentumVector GetMomentum(const corsika::geometry::CoordinateSystem& CS) const { return GetStackData().GetMomentum(GetIndex(), CS); }
+
+    void SetMomentum(const MomentumVector& v) {
+      GetStackData().SetMomentum(GetIndex(), v);
+    }
+  };
+
+  typedef corsika::stack::Stack<QGSJetIIStackData, ParticleInterface> QGSJetIIStack;
+
+} // end namespace corsika::process::qgsjetII
+
+#endif

Processes/QGSJetII/code_generator.py

+#!/usr/bin/env python3
+
+# (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+#
+# See file AUTHORS for a list of contributors.
+#
+# 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.
+
+
+import pickle, sys, itertools
+
+
+
+# loads the pickled particle_db (which is an OrderedDict)
+def load_particledb(filename):
+    with open(filename, "rb") as f:
+        particle_db = pickle.load(f)
+    return particle_db
+
+
+
+# 
+def read_qgsjetII_codes(filename, particle_db):
+    with open(filename) as f:
+        for line in f:
+            line = line.strip()
+            if len(line)==0 or line[0] == '#':
+                continue
+            line = line.split('#')[0]
+            print (line)
+            identifier, model_code, xsType = line.split()
+            try:
+                particle_db[identifier]["qgsjetII_code"] = int(model_code)
+                particle_db[identifier]["qgsjetII_xsType"] = int(xsType)
+            except KeyError as e:
+                raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
+
+
+            
+
+# generates the enum to access qgsjetII particles by readable names
+def generate_qgsjetII_enum(particle_db):
+    output = "enum class QgsjetIICode : int8_t {\n"
+    for identifier, pData in particle_db.items():
+        if pData.get('qgsjetII_code') != None:
+            output += "  {:s} = {:d},\n".format(identifier, pData['qgsjetII_code'])
+    output += "};\n"
+    return output
+
+
+
+# generates the look-up table to convert corsika codes to qgsjetII codes
+def generate_corsika2qgsjetII(particle_db):    
+    string = "std::array<QgsjetIICodeIntType, {:d}> constexpr corsika2qgsjetII = {{\n".format(len(particle_db))
+    for identifier, pData in particle_db.items():
+        modelCode = pData.get("qgsjetII_code", 0)
+        string += "  {:d}, // {:s}\n".format(modelCode, identifier if modelCode else identifier + " (not implemented in QGSJETII)")
+    string += "};\n"
+    return string
+    
+
+
+# generates the look-up table to convert corsika codes to qgsjetII codes
+def generate_corsika2qgsjetII_xsType(particle_db):    
+    string = "std::array<int, {:d}> constexpr corsika2qgsjetIIXStype = {{\n".format(len(particle_db))
+    for identifier, pData in particle_db.items():
+        modelCodeXS = pData.get("qgsjetII_xsType", -1)
+        string += "  {:d}, // {:s}\n".format(modelCodeXS, identifier if modelCodeXS else identifier + " (not implemented in QGSJETII)")
+    string += "};\n"
+    return string
+
+
+# generates the look-up table to convert qgsjetII codes to corsika codes    
+def generate_qgsjetII2corsika(particle_db) :
+    minID = 0
+    for identifier, pData in particle_db.items() :
+        if 'qgsjetII_code' in pData:
+            minID = min(minID, pData['qgsjetII_code'])
+
+    string = "QgsjetIICodeIntType constexpr minQgsjetII = {:d};\n\n".format(minID)
+
+    pDict = {}
+    for identifier, pData in particle_db.items() :
+        if 'qgsjetII_code' in pData:
+            model_code = pData['qgsjetII_code'] - minID
+            pDict[model_code] = identifier
+    
+    nPart = max(pDict.keys()) - min(pDict.keys()) + 1
+    string += "std::array<corsika::particles::Code, {:d}> constexpr qgsjetII2corsika = {{\n".format(nPart)
+    
+    for iPart in range(nPart) :
+        identifier = pDict.get(iPart, "Unknown")
+        qgsID = iPart + minID
+        string += "  corsika::particles::Code::{:s}, // {:d} \n".format(identifier, qgsID)
+    
+    string += "};\n"
+    return string
+
+if __name__ == "__main__":
+    if len(sys.argv) != 3:
+        print("usage: {:s} <particle_db.pkl> <qgsjetII_codes.dat>".format(sys.argv[0]), file=sys.stderr)
+        sys.exit(1)
+        
+    print("code_generator.py for QGSJETII")
+    
+    particle_db = load_particledb(sys.argv[1])
+    read_qgsjetII_codes(sys.argv[2], particle_db)
+    
+    with open("Generated.inc", "w") as f:
+        print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
+        print(generate_qgsjetII_enum(particle_db), file=f)
+        print(generate_corsika2qgsjetII(particle_db), file=f)
+        print(generate_qgsjetII2corsika(particle_db), file=f)
+        print(generate_corsika2qgsjetII_xsType(particle_db), file=f)

Processes/QGSJetII/qgsjet-II-04-codes.dat

+# input file for particle conversion to/from QGSJet
+# the format of this file is: "corsika-identifier" "qgsjet-id" "hadron-class for x-section"
+
+# class 0 (cannot interact)
+Electron    11 0
+Positron   -11 0
+
+# class 1
+Pi0          0 1
+PiPlus       1 1
+PiMinus     -1 1
+# in crmc: all particles from 100 to 999 ???
+Eta         10 1
+# Eta        -10 1 there is no anti-eta
+Rho0        19 1
+
+# class 2
+#Nucleus      40 2
+Neutron      3 2
+AntiNeutron -3 2 
+Proton       2 2
+AntiProton  -2 2 
+# in crmc: 1000 to 9999 ???
+Lambda0      6 2
+Lambda0Bar  -6 2
+LambdaCPlus  9 2
+LambdaCMinusBar -9 2
+
+# class 3
+K0Long      -5 3
+K0           5 3 
+K0Bar       -5 3 
+K0Short      5 3
+# ambiguity between the K0/b and K0s/l
+KPlus        4 3
+KMinus      -4 3
+
+# class 4
+D0           8 4
+D0Bar        8 4
+DPlus        7 4
+DMinus      -7 4
+#DS+/- (340)
+#etac (440)
+#j/psi (441)
+#h_1c0
+#psi'
+#Xi_0c0
+#Xi_1c0
+#Xi_2c0
+
+
+

Processes/QGSJetII/qgsjet-II-04.cc

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/process/qgsjetII/qgsjet-II-04.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <iostream>
+#include <random>
+
+datadir::datadir(const std::string& dir) {
+  if (dir.length() > 130) {
+    std::cerr << "QGSJetII error, will cut datadir \"" << dir
+              << "\" to 130 characters: " << std::endl;
+  }
+  int i = 0;
+  for (i = 0; i < std::min(130, int(dir.length())); ++i) data[i] = dir[i];
+  data[i + 0] = ' ';
+  data[i + 1] = '\0';
+}
+
+double qgran_(int&) {
+  static corsika::random::RNG& rng =
+      corsika::random::RNGManager::GetInstance().GetRandomStream("qgran");
+
+  std::uniform_real_distribution<double> dist;
+  return dist(rng);
+}
+
+void lzmaopenfile_(const char*, int) {}
+void lzmaclosefile_() {}
+void lzmafillarray_(const double&, const int&) {}

Processes/QGSJetII/qgsjet-II-04.h

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_qgsjetII_interface_h_
+#define _include_qgsjetII_interface_h_
+
+#include <string>
+
+//----------------------------------------------
+//  C++ interface for the QGSJetII event generator
+//----------------------------------------------
+// wrapper
+
+extern "C" {
+
+// data memory layout
+
+extern struct { int nsp; } qgarr12_;
+
+const int nptmax = 95000;
+const int iapmax = 208;
+
+extern struct {
+  double esp[nptmax][4];
+  int ich[nptmax];
+} qgarr14_;
+
+extern struct {
+  // c nsf - number of secondary fragments;
+  // c iaf(i) - mass of the i-th fragment
+  int nsf;
+  int iaf[iapmax];
+} qgarr13_;
+
+extern struct {
+  int nwt;
+  int nwp;
+} qgarr55_;
+
+/**
+   Small helper class to provide a data-directory name in the format qgsjetII expects
+ */
+class datadir {
+private:
+  datadir operator=(const std::string& dir);
+  datadir operator=(const datadir&);
+
+public:
+  datadir(const std::string& dir);
+  char data[132];
+};
+
+// functions
+void qgset_();
+void qgaini_(
+    const char* datdir); // Note: there is a length limiation 132 from fortran-qgsjet here
+
+/**
+   @function qgini_
+
+   additional initialization procedure per event
+
+   @parameter e0n  - interaction energy (per hadron/nucleon),
+   @parameter icp0 - hadron type (+-1 - pi+-, +-2 - p(p~), +-3 - n(n~), +-4 - K+-, +-5 -
+   K_l/s),
+   @parameter iap  - projectile mass number (1 - for a hadron),
+   @parameter iat  - target mass number
+*/
+void qgini_(const double& e0n, const int& icp0, const int& iap, const int& iat);
+
+/**
+   @function qgconf_
+
+   generate one event configuration
+*/
+void qgconf_();
+
+/**
+   @function qgsect_
+
+   hadron-nucleus (hadron-nucleus) particle production cross section
+
+   @parameter e0n lab. energy per projectile nucleon (hadron)
+   @parameter icz hadron class (1 - pion, 2 - nucleon, 3 - kaon)
+   @parameter iap projectile mass number (1=<iap<=iapmax),
+   @parameter iat target mass number     (1=<iat<=iapmax)
+ */
+double qgsect_(const double& e0n, const int& icz, const int& iap0, const int& iat0);
+
+/**
+   @function qgran
+
+   link to random number generation
+ */
+double qgran_(int&);
+
+/**
+   dummy function from CRMC
+ */
+void lzmaopenfile_(const char* name, int length);
+void lzmaclosefile_();
+void lzmafillarray_(const double& dum, const int& idum);
+}
+
+#endif

Processes/QGSJetII/testQGSJetII.cc

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/process/qgsjetII/Interaction.h>
+#include <corsika/process/qgsjetII/ParticleConversion.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process::qgsjetII;
+
+TEST_CASE("QgsjetII", "[processes]") {
+
+  SECTION("QgsjetII -> Corsika") {
+    REQUIRE(particles::PiPlus::GetCode() == process::qgsjetII::ConvertFromQgsjetII(
+                                                process::qgsjetII::QgsjetIICode::PiPlus));
+  }
+
+  SECTION("Corsika -> QgsjetII") {
+    REQUIRE(process::qgsjetII::ConvertToQgsjetII(particles::PiMinus::GetCode()) ==
+            process::qgsjetII::QgsjetIICode::PiMinus);
+    REQUIRE(process::qgsjetII::ConvertToQgsjetIIRaw(particles::Proton::GetCode()) == 2);
+  }
+
+  SECTION("canInteractInQgsjetII") {
+
+    REQUIRE(process::qgsjetII::CanInteract(particles::Proton::GetCode()));
+    REQUIRE(process::qgsjetII::CanInteract(particles::Code::KPlus));
+    REQUIRE(process::qgsjetII::CanInteract(particles::Nucleus::GetCode()));
+    // REQUIRE(process::qgsjetII::CanInteract(particles::Helium::GetCode()));
+
+    REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::EtaC::GetCode()));
+    REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::SigmaC0::GetCode()));
+  }
+
+  SECTION("cross-section type") {
+
+    REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) == 2);
+    REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) == 3);
+    REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) == 2);
+    REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) == 1);
+  }
+}
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/qgsjetII/qgsjet-II-04.h>
+
+using namespace corsika::units::si;
+using namespace corsika::units;
+
+TEST_CASE("QgsjetIIInterface", "[processes]") {
+
+  // setup environment, geometry
+  environment::Environment<environment::IMediumModel> env;
+  auto& universe = *(env.GetUniverse());
+
+  auto theMedium =
+      environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+          geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+          1_km * std::numeric_limits<double>::infinity());
+
+  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+  theMedium->SetModelProperties<MyHomogeneousModel>(
+      1_kg / (1_m * 1_m * 1_m),
+      environment::NuclearComposition(
+          std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+
+  auto const* nodePtr = theMedium.get();
+  universe.AddChild(std::move(theMedium));
+
+  const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+
+  random::RNGManager::GetInstance().RegisterRandomStream("qgran");
+
+  SECTION("InteractionInterface") {
+
+    setup::Stack stack;
+    const HEPEnergyType E0 = 100_GeV;
+    HEPMomentumType P0 =
+        sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+    auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+    geometry::Point pos(cs, 0_m, 0_m, 0_m);
+    auto particle =
+        stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+                                     corsika::stack::MomentumVector, geometry::Point,
+                                     units::si::TimeType, unsigned int, unsigned int>{
+            particles::Code::Nucleus, E0, plab, pos, 0_ns, 16, 8});
+    // corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+    //	  particles::Code::PiPlus, E0, plab, pos, 0_ns});
+
+    particle.SetNode(nodePtr);
+    corsika::stack::SecondaryView view(particle);
+    auto projectile = view.GetProjectile();
+
+    Interaction model;
+    model.Init();
+    [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
+    [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+  }
+}