corsika/detail/modules/proposal/ContinuousProcess.inl

+
 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <PROPOSAL/PROPOSAL.h>
 
 #include <corsika/media/IMediumModel.hpp>
 #include <corsika/modules/proposal/ContinuousProcess.hpp>
-#include <corsika/modules/proposal/Interaction.hpp>
+#include <corsika/modules/proposal/InteractionModel.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
 #include <corsika/framework/core/Logging.hpp>
 
-#include <corsika/setup/SetupEnvironment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <iostream>
-
 namespace corsika::proposal {
 
-  inline void ContinuousProcess::buildCalculator(Code code,
-                                                 NuclearComposition const& comp) {
+  template <typename TOutput>
+  inline void ContinuousProcess<TOutput>::buildCalculator(Code code,
+                                                          size_t const& component_hash) {
     // search crosssection builder for given particle
     auto p_cross = cross.find(code);
     if (p_cross == cross.end())
       throw std::runtime_error("PROPOSAL could not find corresponding builder");
+    if (code == Code::Photon) return; // no continuous builders needed for photons
 
     // interpolate the crosssection for given media and energy cut. These may
-    // take some minutes if you have to build the tables and cannot read the
+    // take some minutes if you have to build the tables and cannot read the tables
     // from disk
-    auto const emCut = get_energy_threshold(
-        code); //! energy thresholds globally defined for individual particles
-    auto c = p_cross->second(media.at(comp.getHash()), emCut);
+    auto c = p_cross->second(media.at(component_hash), proposal_energycutsettings[code]);
+
+    // choose multiple scattering model
+    static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::MoliereInterpol;
 
     // Build displacement integral and scattering object and interpolate them too and
     // saved in the calc map by a key build out of a hash of composed of the component and
     // particle code.
-    auto disp = PROPOSAL::make_displacement(c, true);
-    auto scatter =
-        PROPOSAL::make_scattering("highland", particle[code], media.at(comp.getHash()));
-    calc[std::make_pair(comp.getHash(), code)] =
-        std::make_tuple(std::move(disp), std::move(scatter));
+    auto calculator = Calculator{PROPOSAL::make_displacement(c, true),
+                                 PROPOSAL::make_multiple_scattering(
+                                     ms_type, particle[code], media.at(component_hash))};
+    calc[std::make_pair(component_hash, code)] = std::move(calculator);
   }
 
-  template <>
-  inline ContinuousProcess::ContinuousProcess(setup::Environment const& _env)
-      : ProposalProcessBase(_env) {}
+  template <typename TOutput>
+  template <typename TEnvironment, typename... TOutputArgs>
+  inline ContinuousProcess<TOutput>::ContinuousProcess(TEnvironment const& _env,
+                                                       TOutputArgs&&... args)
+      : ProposalProcessBase(_env)
+      , TOutput(args...) {
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildCalculator(particle_code, medium.first);
+      }
+    }
+  }
 
-  template <>
-  inline void ContinuousProcess::scatter(setup::Stack::particle_type& vP,
-                                         HEPEnergyType const& loss,
-                                         GrammageType const& grammage) {
+  template <typename TOutput>
+  template <typename TParticle>
+  inline void ContinuousProcess<TOutput>::scatter(Step<TParticle>& step,
+                                                  HEPEnergyType const& loss,
+                                                  GrammageType const& grammage) {
 
     // get or build corresponding calculators
-    auto c = getCalculator(vP, calc);
+    auto c = getCalculator(step.getParticlePre(), calc);
 
-    // Cast corsika vector to proposal vector
-    auto vP_dir = vP.getDirection();
-    auto d = vP_dir.getComponents();
-    auto direction = PROPOSAL::Vector3D(d.getX().magnitude(), d.getY().magnitude(),
-                                        d.getZ().magnitude());
+    auto initial_particle_dir = step.getDirectionPre();
 
-    auto E_f = vP.getEnergy() - loss;
+    auto E_i_total = step.getEkinPre() + step.getParticlePre().getMass();
+    auto E_f_total = E_i_total - loss;
 
-    // draw random numbers required for scattering process
+    // sample scattering_angle, which is the combination sqrt(theta_x^2 + theta_y^2),
+    // where theta_x and theta_y itself are independent angles drawn from the multiple
+    // scattering distribution
     std::uniform_real_distribution<double> distr(0., 1.);
-    auto rnd = array<double, 4>();
-    for (auto& it : rnd) it = distr(RNG_);
-
-    // calculate deflection based on particle energy, loss
-    auto [mean_dir, final_dir] = get<eSCATTERING>(c->second)->Scatter(
-        grammage / 1_g * square(1_cm), vP.getEnergy() / 1_MeV, E_f / 1_MeV, direction,
-        rnd);
-
-    // TODO: neglect mean direction deflection because Trajectory is a const ref
-    (void)mean_dir;
+    auto scattering_angle = (c->second).scatter->CalculateScatteringAngle2D(
+        grammage / 1_g * square(1_cm), E_i_total / 1_MeV, E_f_total / 1_MeV, distr(RNG_),
+        distr(RNG_));
+
+    auto const& root = initial_particle_dir.getCoordinateSystem();
+
+    // construct vector that is normal to initial direction.
+    DirectionVector normal_vec{root, {0, 0, 0}};
+    if (initial_particle_dir.getX(root) > 0.1) {
+      normal_vec = {
+          root, {-initial_particle_dir.getY(root), initial_particle_dir.getX(root), 0}};
+    } else {
+      // if x is small, use y and z to construct normal vector
+      normal_vec = {
+          root, {0, -initial_particle_dir.getZ(root), initial_particle_dir.getY(root)}};
+    }
+
+    // rotation of zenith by moliere_angle
+    CoordinateSystemPtr rotated1 =
+        make_rotation(root, normal_vec.getComponents(), scattering_angle);
+
+    // rotation of azimuth by random angle between 0 and 2*PI
+    std::uniform_real_distribution<double> distr_azimuth(0., 2 * M_PI);
+    CoordinateSystemPtr rotated2 = make_rotation(
+        rotated1, initial_particle_dir.getComponents(), distr_azimuth(RNG_));
+
+    DirectionVector scattered_particle_dir{root,
+                                           initial_particle_dir.getComponents(rotated2)};
 
     // update particle direction after continuous loss caused by multiple
     // scattering
-    auto vec = QuantityVector(final_dir.GetX() * E_f, final_dir.GetY() * E_f,
-                              final_dir.GetZ() * E_f);
-    vP.setMomentum(MomentumVector(vP_dir.getCoordinateSystem(), vec));
+    DirectionVector diff_dir_ = scattered_particle_dir - initial_particle_dir;
+    step.add_dU(diff_dir_);
   }
 
-  template <>
-  inline ProcessReturn ContinuousProcess::doContinuous(setup::Stack::particle_type& vP,
-                                                       setup::Trajectory const& vT,
-                                                       bool const) {
-
-    if (!canInteract(vP.getPID())) return ProcessReturn::Ok;
-    if (vT.getLength() == 0_m) return ProcessReturn::Ok;
-
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn ContinuousProcess<TOutput>::doContinuous(Step<TParticle>& step,
+                                                                bool const) {
+    if (!canInteract(step.getParticlePre().getPID())) return ProcessReturn::Ok;
+    if (step.getDisplacement().getSquaredNorm() == static_pow<2>(0_m))
+      return ProcessReturn::Ok;
+    if (step.getParticlePre().getPID() == Code::Photon)
+      return ProcessReturn::Ok; // no continuous energy losses, no scattering for photons
     // calculate passed grammage
-    auto dX =
-        vP.getNode()->getModelProperties().getIntegratedGrammage(vT, vT.getLength());
+    auto dX = step.getParticlePre().getNode()->getModelProperties().getIntegratedGrammage(
+        step.getStraightTrack());
 
     // get or build corresponding track integral calculator and solve the
     // integral
-    auto c = getCalculator(vP, calc);
-    auto final_energy = get<eDISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
-                            vP.getEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
-                        1_MeV;
-    auto dE = vP.getEnergy() - final_energy;
-    energy_lost_ += dE;
+    auto c = getCalculator(step.getParticlePre(), calc);
+    auto E_i_total = (step.getEkinPre() + step.getParticlePre().getMass());
+    auto E_f_total = (c->second).disp->UpperLimitTrackIntegral(
+                         E_i_total * (1 / 1_MeV), dX * ((1 / 1_g) * 1_cm * 1_cm)) *
+                     1_MeV;
+    auto dE = E_i_total - E_f_total;
 
     // if the particle has a charge take multiple scattering into account
-    if (vP.getChargeNumber() != 0) scatter(vP, dE, dX);
-    vP.setEnergy(final_energy);
-    vP.setMomentum(vP.getMomentum() * vP.getEnergy() / vP.getMomentum().getNorm());
+    if (step.getParticlePre().getChargeNumber() != 0) scatter(step, dE, dX);
+    step.add_dEkin(-dE); // on the stack, this is just kinetic energy, E-m
+
+    // also send to output
+    TOutput::write(step.getPositionPre(), step.getPositionPost(),
+                   step.getParticlePre().getPID(),
+                   step.getParticlePre().getWeight() * dE);
+
     return ProcessReturn::Ok;
   }
 
-  template <>
-  inline LengthType ContinuousProcess::getMaxStepLength(
-      setup::Stack::particle_type const& vP, setup::Trajectory const& vT) {
+  template <typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType ContinuousProcess<TOutput>::getMaxStepLength(
+      TParticle const& vP, TTrajectory const& track) {
     auto const code = vP.getPID();
     if (!canInteract(code)) return meter * std::numeric_limits<double>::infinity();
+    if (code == Code::Photon)
+      return meter *
+             std::numeric_limits<double>::infinity(); // no step limitation for photons
 
-    // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
-    // hyper parameter which must be adjusted.
+    // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be
+    // a hyper parameter which must be adjusted.
     //
     auto const energy = vP.getEnergy();
     auto const energy_lim = std::max(
         energy * 0.9, // either 10% relative loss max., or
-        get_energy_threshold(
-            code) // energy thresholds globally defined for individual particles
-            *
-            0.99 // need to go 1% below global e-cut to assure removal in ParticleCut. The
-                 // 1% does not matter since at cut-time the entire energy is removed.
+        (get_kinetic_energy_propagation_threshold(code) +
+         get_mass(code)) // energy thresholds globally defined for individual particles
+            * 0.9999     // need to go slightly below global e-cut to assure removal
+                         // in ParticleCut. This does not matter since at cut-time
+                         // the entire energy is removed.
     );
 
     // solving the track integral for giving energy lim
     auto c = getCalculator(vP, calc);
-    auto grammage = get<eDISPLACEMENT>(c->second)->SolveTrackIntegral(
-                        energy / 1_MeV, energy_lim / 1_MeV) *
-                    1_g / square(1_cm);
+    auto grammage =
+        (c->second).disp->SolveTrackIntegral(energy / 1_MeV, energy_lim / 1_MeV) * 1_g /
+        square(1_cm);
 
     // return it in distance aequivalent
-    auto dist = vP.getNode()->getModelProperties().getArclengthFromGrammage(vT, grammage);
+    auto dist =
+        vP.getNode()->getModelProperties().getArclengthFromGrammage(track, grammage);
     CORSIKA_LOG_TRACE("PROPOSAL::getMaxStepLength X={} g/cm2, l={} m ",
                       grammage / 1_g * square(1_cm), dist / 1_m);
+    if (dist < 0_m) {
+      CORSIKA_LOG_WARN(
+          "PROPOSAL::getMaxStepLength calculated a negative step length of l={} m. "
+          "Return 0_m instead.",
+          dist / 1_m);
+      return 0_m;
+    }
     return dist;
   }
 
-  inline void ContinuousProcess::showResults() const {
-    CORSIKA_LOG_DEBUG(
-        " ******************************\n"
-        " PROCESS::ContinuousProcess: \n"
-        " energy lost dE (GeV)      :  {}",
-        energy_lost_ / 1_GeV);
+  template <typename TOutput>
+  inline YAML::Node ContinuousProcess<TOutput>::getConfig() const {
+    return YAML::Node();
   }
 
-  inline void ContinuousProcess::reset() { energy_lost_ = 0_GeV; }
-
 } // namespace corsika::proposal

corsika/detail/modules/proposal/HadronicPhotonModel.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <tuple>
+#include <random>
+
+namespace corsika::proposal {
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  inline HadronicPhotonModel<THadronicLEModel, THadronicHEModel>::HadronicPhotonModel(
+      THadronicLEModel& _hadintLE, THadronicHEModel& _hadintHE,
+      HEPEnergyType const& _heenthresholdNN)
+      : leHadronicInteraction_(_hadintLE)
+      , heHadronicInteraction_(_hadintHE)
+      , heHadronicModelThresholdLabNN_(_heenthresholdNN) {
+    // check validity of threshold assuming photon-nucleon
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtS =
+        calculate_com_energy(_heenthresholdNN, Rho0::mass, Proton::mass);
+    if (!heHadronicInteraction_.isValid(Code::Rho0, Code::Proton, sqrtS)) {
+      CORSIKA_LOGGER_CRITICAL(
+          logger_,
+          "Invalid energy threshold for hadron interaction model. threshold_lab= {} GeV, "
+          "threshold_com={} GeV",
+          _heenthresholdNN / 1_GeV, sqrtS / 1_GeV);
+      throw std::runtime_error("Configuration error!");
+    }
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "Threshold for HE hadronic interactions in proposal set to Elab={} GeV",
+        _heenthresholdNN / 1_GeV);
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TStackView>
+  inline ProcessReturn
+  HadronicPhotonModel<THadronicLEModel, THadronicHEModel>::doHadronicPhotonInteraction(
+      TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
+      Code const& targetId) {
+
+    //  temporarily add to stack, will be removed after interaction in DoInteraction
+    typename TStackView::inner_stack_value_type photonStack;
+    Point const pDummy(labCS, {0_m, 0_m, 0_m});
+    TimeType const tDummy = 0_ns;
+    // target at rest
+    FourMomentum const targetP4(get_mass(targetId),
+                                MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+    auto hadronicPhoton = photonStack.addParticle(
+        std::make_tuple(Code::Photon, photonP4.getTimeLikeComponent(),
+                        photonP4.getSpaceLikeComponents().normalized(), pDummy, tDummy));
+    hadronicPhoton.setNode(view.getProjectile().getNode());
+    // create inelastic interaction of the hadronic photon
+    // create new StackView for the photon
+    TStackView photon_secondaries(hadronicPhoton);
+
+    // call inner hadronic event generator
+    CORSIKA_LOGGER_TRACE(logger_, "{} + {} interaction. Ekinlab = {} GeV", Code::Photon,
+                         targetId, photonP4.getTimeLikeComponent() / 1_GeV);
+    // check if had. model can handle configuration
+
+    if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
+      CORSIKA_LOGGER_TRACE(logger_, "HE photo-hadronic interaction!");
+      auto const sqrtSNN = (photonP4 + targetP4 / get_nucleus_A(targetId)).getNorm();
+      CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={} GeV", sqrtSNN / 1_GeV);
+      // when Sibyll is used for hadronic interactions Argon cannot be used as target
+      // nucleus. Since PROPOSAL has a non-zero cross section for Argon
+      // targets we have to check here if the model can handle Argon (see Issue #498)
+      if (!heHadronicInteraction_.isValid(Code::Rho0, targetId, sqrtSNN)) {
+        CORSIKA_LOGGER_WARN(
+            logger_,
+            "HE interaction model cannot handle configuration in photo-hadronic "
+            "interaction! projectile={}, target={} (A={}, Z={}), sqrt(S) per "
+            "nuc.={:8.2f} "
+            "GeV. Skipping secondary production!",
+            Code::Rho0, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
+            sqrtSNN / 1_GeV);
+        return ProcessReturn::Ok;
+      }
+      heHadronicInteraction_.doInteraction(photon_secondaries, Code::Rho0, targetId,
+                                           photonP4, targetP4);
+    } else {
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "LE photo-hadronic interaction! implemented via SOPHIA "
+                           "assuming a single nucleon as target");
+      // sample nucleon from nucleus A,Z
+      double const fProtons = get_nucleus_Z(targetId) / double(get_nucleus_A(targetId));
+      double const fNeutrons = 1. - fProtons;
+      std::discrete_distribution<int> nucleonChannelDist{fProtons, fNeutrons};
+      corsika::default_prng_type& rng =
+          corsika::RNGManager<>::getInstance().getRandomStream("proposal");
+      Code const nucleonId = (nucleonChannelDist(rng) ? Code::Neutron : Code::Proton);
+      // target passed to SOPHIA needs to be exactly on-shell!
+      FourMomentum const nucleonP4(get_mass(nucleonId),
+                                   MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+      CORSIKA_LOGGER_DEBUG(logger_,
+                           "selected {} as target nucleon (f_proton, f_neutron)={},{}",
+                           nucleonId, fProtons, fNeutrons);
+      auto const sqrtSNN = (photonP4 + nucleonP4).getNorm();
+      CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={} GeV", sqrtSNN / 1_GeV);
+
+      if (!leHadronicInteraction_.isValid(Code::Photon, nucleonId, sqrtSNN)) {
+        CORSIKA_LOGGER_WARN(
+            logger_,
+            "LE interaction model cannot handle configuration in photo-hadronic "
+            "interaction! projectile={}, target={} (A={}, Z={}), sqrt(S) per "
+            "nuc.={:8.2f} "
+            "GeV. Skipping secondary production!",
+            Code::Photon, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
+            sqrtSNN / 1_GeV);
+        return ProcessReturn::Ok;
+      }
+      leHadronicInteraction_.doInteraction(photon_secondaries, Code::Photon, nucleonId,
+                                           photonP4, nucleonP4);
+    }
+    for (const auto& pSec : photon_secondaries) {
+      auto const p3lab = pSec.getMomentum();
+      Code const pid = pSec.getPID();
+      HEPEnergyType const secEkin =
+          calculate_kinetic_energy(p3lab.getNorm(), get_mass(pid));
+      view.addSecondary(std::make_tuple(pid, secEkin, p3lab.normalized()));
+    }
+    return ProcessReturn::Ok;
+  }
+} // namespace corsika::proposal

corsika/detail/modules/proposal/Interaction.inl

-/*
- * (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/media/IMediumModel.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/modules/proposal/Interaction.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <corsika/setup/SetupEnvironment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <limits>
-#include <memory>
-#include <random>
-#include <tuple>
-
-namespace corsika::proposal {
-
-  template <>
-  inline Interaction::Interaction(setup::Environment const& _env)
-      : ProposalProcessBase(_env) {}
-
-  inline void Interaction::buildCalculator(Code code, NuclearComposition const& comp) {
-    // search crosssection builder for given particle
-    auto p_cross = cross.find(code);
-    if (p_cross == cross.end())
-      throw std::runtime_error("PROPOSAL could not find corresponding builder");
-
-    // interpolate the crosssection for given media and energy cut. These may
-    // take some minutes if you have to build the tables and cannot read the
-    // from disk
-    auto const emCut = get_energy_threshold(
-        code); //! energy thresholds globally defined for individual particles
-
-    auto c = p_cross->second(media.at(comp.getHash()), emCut);
-
-    // Look which interactions take place and build the corresponding
-    // interaction and secondarie builder. The interaction integral will
-    // interpolated too and saved in the calc map by a key build out of a hash
-    // of composed of the component and particle code.
-    auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
-    calc[std::make_pair(comp.getHash(), code)] = std::make_tuple(
-        PROPOSAL::make_secondaries(inter_types, particle[code], media.at(comp.getHash())),
-        PROPOSAL::make_interaction(c, true));
-  }
-
-  template <>
-  inline ProcessReturn Interaction::doInteraction(setup::StackView& view) {
-
-    auto const projectile = view.getProjectile();
-
-    if (canInteract(projectile.getPID())) {
-
-      // get or build corresponding calculators
-      auto c = getCalculator(projectile, calc);
-
-      // get the rates of the interaction types for every component.
-      std::uniform_real_distribution<double> distr(0., 1.);
-
-      // sample a interaction-type, loss and component
-      auto rates = get<eINTERACTION>(c->second)->Rates(projectile.getEnergy() / 1_MeV);
-      auto [type, comp_ptr, v] = get<eINTERACTION>(c->second)->SampleLoss(
-          projectile.getEnergy() / 1_MeV, rates, distr(RNG_));
-
-      // Read how much random numbers are required to calculate the secondaries.
-      // Calculate the secondaries and deploy them on the corsika stack.
-      auto rnd =
-          vector<double>(get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type));
-      for (auto& it : rnd) it = distr(RNG_);
-      Point const& place = projectile.getPosition();
-      CoordinateSystemPtr const& labCS = place.getCoordinateSystem();
-
-      auto point = PROPOSAL::Vector3D(place.getX(labCS) / 1_cm, place.getY(labCS) / 1_cm,
-                                      place.getZ(labCS) / 1_cm);
-      auto projectile_dir = projectile.getDirection();
-      auto d = projectile_dir.getComponents(labCS);
-      auto direction = PROPOSAL::Vector3D(d.getX().magnitude(), d.getY().magnitude(),
-                                          d.getZ().magnitude());
-      auto loss = make_tuple(static_cast<int>(type), point, direction,
-                             v * projectile.getEnergy() / 1_MeV, 0.);
-      auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries(
-          projectile.getEnergy() / 1_MeV, loss, *comp_ptr, rnd);
-      for (auto& s : sec) {
-        auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
-        auto vecProposal = get<PROPOSAL::Loss::DIRECTION>(s);
-        auto vec = QuantityVector(vecProposal.GetX() * E, vecProposal.GetY() * E,
-                                  vecProposal.GetZ() * E);
-        auto p = MomentumVector(labCS, vec);
-        auto sec_code =
-            convert_from_PDG(static_cast<PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
-        view.addSecondary(
-            make_tuple(sec_code, E, p, projectile.getPosition(), projectile.getTime()));
-      }
-    }
-    return ProcessReturn::Ok;
-  }
-
-  template <>
-  inline GrammageType Interaction::getInteractionLength(
-      setup::Stack::particle_type const& projectile) {
-
-    if (canInteract(projectile.getPID())) {
-      auto c = getCalculator(projectile, calc);
-      return get<eINTERACTION>(c->second)->MeanFreePath(projectile.getEnergy() / 1_MeV) *
-             1_g / (1_cm * 1_cm);
-    }
-    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-  }
-} // namespace corsika::proposal

corsika/detail/modules/proposal/InteractionModel.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <limits>
+#include <memory>
+#include <random>
+#include <tuple>
+#include <PROPOSAL/particle/Particle.h>
+#include <corsika/modules/proposal/InteractionModel.hpp>
+
+namespace corsika::proposal {
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TEnvironment>
+  inline InteractionModel<THadronicLEModel, THadronicHEModel>::InteractionModel(
+      TEnvironment const& _env, THadronicLEModel& _hadintLE, THadronicHEModel& _hadintHE,
+      HEPEnergyType const& _enthreshold)
+      : ProposalProcessBase(_env)
+      , HadronicPhotonModel<THadronicLEModel, THadronicHEModel>(_hadintLE, _hadintHE,
+                                                                _enthreshold) {
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildCalculator(particle_code, medium.first);
+      }
+    }
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  inline void InteractionModel<THadronicLEModel, THadronicHEModel>::buildCalculator(
+      Code code, size_t const& component_hash) {
+    // search crosssection builder for given particle
+    auto p_cross = cross.find(code);
+    if (p_cross == cross.end())
+      throw std::runtime_error("PROPOSAL could not find corresponding builder");
+
+    // interpolate the crosssection for given media and energy cut. These may
+    // take some minutes if you have to build the tables and cannot read the tables
+    // from disk
+    auto c = p_cross->second(media.at(component_hash), proposal_energycutsettings[code]);
+
+    // Look which interactions take place and build the corresponding
+    // interaction and secondary builder. The interaction integral will
+    // interpolated too and saved in the calc map by a key build out of a hash
+    // of composed of the component and particle code.
+    auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
+    calc_[std::make_pair(component_hash, code)] = std::make_tuple(
+        PROPOSAL::make_secondaries(inter_types, particle[code], media.at(component_hash)),
+        PROPOSAL::make_interaction(c, true, true),
+        std::make_unique<LPM_calculator>(media.at(component_hash), code, inter_types));
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TStackView>
+  inline ProcessReturn
+  InteractionModel<THadronicLEModel, THadronicHEModel>::doInteraction(
+      TStackView& view, Code const projectileId,
+      [[maybe_unused]] FourMomentum const& projectileP4) {
+
+    auto const projectile = view.getProjectile();
+
+    if (canInteract(projectileId)) {
+
+      // get or build corresponding calculators
+      auto c = getCalculator(projectile, calc_);
+
+      // get the rates of the interaction types for every component.
+      std::uniform_real_distribution<double> distr(0., 1.);
+
+      // sample a interaction-type, loss and component
+      auto rates =
+          std::get<eINTERACTION>(c->second)->Rates(projectile.getEnergy() / 1_MeV);
+      auto [type, target_hash, v] = std::get<eINTERACTION>(c->second)->SampleLoss(
+          projectile.getEnergy() / 1_MeV, rates, distr(RNG_));
+
+      // TODO: This should become obsolete as soon #482 is fixed
+      if (type == PROPOSAL::InteractionType::Undefined) {
+        CORSIKA_LOG_WARN(
+            "PROPOSAL: No particle interaction possible. "
+            "Put initial particle back on stack.");
+        view.addSecondary(std::make_tuple(projectileId,
+                                          projectile.getEnergy() - get_mass(projectileId),
+                                          projectile.getDirection()));
+        return ProcessReturn::Ok;
+      }
+
+      // Read how much random numbers are required to calculate the secondaries.
+      // Calculate the secondaries and deploy them on the corsika stack.
+      auto rnd = std::vector<double>(
+          std::get<eSECONDARIES>(c->second)->RequiredRandomNumbers(type));
+      for (auto& it : rnd) it = distr(RNG_);
+      Point const& place = projectile.getPosition();
+      CoordinateSystemPtr const& labCS = place.getCoordinateSystem();
+
+      auto point = PROPOSAL::Cartesian3D(
+          place.getX(labCS) / 1_cm, place.getY(labCS) / 1_cm, place.getZ(labCS) / 1_cm);
+      auto projectile_dir = projectile.getDirection();
+      auto d = projectile_dir.getComponents(labCS);
+      auto direction = PROPOSAL::Cartesian3D(d.getX().magnitude(), d.getY().magnitude(),
+                                             d.getZ().magnitude());
+
+      auto loss = PROPOSAL::StochasticLoss(
+          static_cast<int>(type), v * projectile.getEnergy() / 1_MeV, point, direction,
+          projectile.getTime() / 1_s, 0., projectile.getEnergy() / 1_MeV);
+      PROPOSAL::Component target;
+      if (type != PROPOSAL::InteractionType::Ioniz)
+        target = PROPOSAL::Component::GetComponentForHash(target_hash);
+
+      auto sec =
+          std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
+
+      // Check for LPM effect suppression!
+
+      if (std::get<eLPM_SUPPRESSION>(c->second)) {
+        auto lpm_suppression = CheckForLPM(
+            *std::get<eLPM_SUPPRESSION>(c->second), projectile.getEnergy(), type, sec,
+            projectile.getNode()->getModelProperties().getMassDensity(place), target, v);
+        if (lpm_suppression) {
+          // Discard interaction - put initial particle back on stack
+          CORSIKA_LOG_DEBUG("LPM suppression detected!");
+          view.addSecondary(std::make_tuple(
+              projectileId, projectile.getEnergy() - get_mass(projectileId),
+              projectile.getDirection()));
+          return ProcessReturn::Ok;
+        }
+      }
+
+      for (auto& s : sec) {
+        auto E = s.energy * 1_MeV;
+        auto vecProposal = s.direction;
+        auto dir = DirectionVector(
+            labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
+
+        if (static_cast<PROPOSAL::ParticleType>(s.type) ==
+            PROPOSAL::ParticleType::Hadron) {
+          FourMomentum const photonP4(E, E * dir);
+          // for PROPOSAL media target.GetAtomicNum() returns the atomic number in units
+          // of atomic mass, so Nitrogen is 14.0067. When using media in CORSIKA the same
+          // field is filled with the number of nucleons (ie. 14 for Nitrogen) To be sure
+          // we explicitly cast to int
+          auto const A = int(target.GetAtomicNum());
+          auto const Z = int(target.GetNucCharge());
+          Code const targetId = get_nucleus_code(A, Z);
+          CORSIKA_LOGGER_DEBUG(
+              logger_,
+              "photo-hadronic interaction of projectile={} with target={}! Energy={} GeV",
+              projectileId, targetId, E / 1_GeV);
+          this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
+        } else {
+          auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
+          // use mass provided by PROPOSAL to ensure correct conversion to kinetic energy
+          auto massProposal =
+              PROPOSAL::ParticleDef::GetParticleDefForType(s.type).mass * 1_MeV;
+          view.addSecondary(std::make_tuple(sec_code, E - massProposal, dir));
+        }
+      }
+    }
+    return ProcessReturn::Ok;
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  template <typename TParticle>
+  inline CrossSectionType
+  InteractionModel<THadronicLEModel, THadronicHEModel>::getCrossSection(
+      TParticle const& projectile, Code const projectileId,
+      FourMomentum const& projectileP4) {
+
+    // ==============================================
+    // this block better diappears. RU 26.10.2021
+    //
+    // determine the volume where the particle is (last) known to be
+    auto const* currentLogicalNode = projectile.getNode();
+    NuclearComposition const& composition =
+        currentLogicalNode->getModelProperties().getNuclearComposition();
+    auto const meanMass = composition.getAverageMassNumber() * constants::u;
+    // ==============================================
+
+    if (canInteract(projectileId)) {
+      auto c = getCalculator(projectile, calc_);
+      return meanMass / (std::get<eINTERACTION>(c->second)->MeanFreePath(
+                             projectileP4.getTimeLikeComponent() / 1_MeV) *
+                         1_g / (1_cm * 1_cm));
+    }
+
+    return CrossSectionType::zero();
+  }
+
+  template <typename THadronicLEModel, typename THadronicHEModel>
+  bool InteractionModel<THadronicLEModel, THadronicHEModel>::CheckForLPM(
+      const LPM_calculator& lpm_calculator, const HEPEnergyType projectile_energy,
+      const PROPOSAL::InteractionType type,
+      const std::vector<PROPOSAL::ParticleState>& sec, const MassDensityType mass_density,
+      const PROPOSAL::Component& target, const double v) {
+
+    double suppression_factor;
+    double current_density = mass_density * ((1_cm * 1_cm * 1_cm) / 1_g);
+
+    if (type == PROPOSAL::InteractionType::Photopair && lpm_calculator.photo_pair_lpm_) {
+      if (sec.size() != 2)
+        throw std::runtime_error(
+            "Invalid number of secondaries for Photopair in CheckForLPM");
+      auto x =
+          sec[0].energy / (sec[0].energy + sec[1].energy); // particle energy asymmetry
+      double density_correction = current_density / lpm_calculator.mass_density_baseline_;
+      suppression_factor = lpm_calculator.photo_pair_lpm_->suppression_factor(
+          projectile_energy / 1_MeV, x, target, density_correction);
+    } else if (type == PROPOSAL::InteractionType::Brems && lpm_calculator.brems_lpm_) {
+      double density_correction = current_density / lpm_calculator.mass_density_baseline_;
+      suppression_factor = lpm_calculator.brems_lpm_->suppression_factor(
+          projectile_energy / 1_MeV, v, target, density_correction);
+    } else if (type == PROPOSAL::InteractionType::Epair && lpm_calculator.epair_lpm_) {
+      if (sec.size() != 3)
+        throw std::runtime_error(
+            "Invalid number of secondaries for EPair in CheckForLPM");
+      double rho = (sec[1].energy - sec[2].energy) /
+                   (sec[1].energy + sec[2].energy); // todo: check
+      // it would be better if these variables are calculated within PROPOSAL
+      double beta = (v * v) / (2 * (1 - v));
+      double xi = (lpm_calculator.particle_mass_ * lpm_calculator.particle_mass_) /
+                  (PROPOSAL::ME * PROPOSAL::ME) * v * v / 4 * (1 - rho * rho) / (1 - v);
+      double density_correction = current_density / lpm_calculator.mass_density_baseline_;
+      suppression_factor = lpm_calculator.epair_lpm_->suppression_factor(
+          projectile_energy / 1_MeV, v, rho * rho, beta, xi, density_correction);
+    } else {
+      return false;
+    }
+
+    // rejection sampling
+    std::uniform_real_distribution<double> distr(0., 1.);
+    auto rnd_num = distr(RNG_);
+
+    CORSIKA_LOGGER_TRACE(logger_,
+                         "Checking for LPM suppression! energy: {}, v: {}, type: {}, "
+                         "target: {}, mass_density: {}, mass_density_baseline: {}, "
+                         "suppression_factor: {}, rnd: {}, result: {}, pmass: {} ",
+                         projectile_energy, v, type, target.GetName(),
+                         mass_density / (1_g / (1_cm * 1_cm * 1_cm)),
+                         lpm_calculator.mass_density_baseline_, suppression_factor,
+                         rnd_num, (rnd_num > suppression_factor),
+                         lpm_calculator.particle_mass_);
+
+    if (rnd_num > suppression_factor)
+      return true;
+    else
+      return false;
+  }
+} // namespace corsika::proposal

corsika/detail/modules/proposal/ProposalProcessBase.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/media/IMediumModel.hpp>
@@ -11,10 +10,7 @@
 #include <corsika/modules/proposal/ProposalProcessBase.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
-
-#include <corsika/setup/SetupEnvironment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
 
 #include <cstdlib>
 #include <iostream>
@@ -26,22 +22,41 @@
 namespace corsika::proposal {
 
   inline bool ProposalProcessBase::canInteract(Code pcode) const {
-    if (std::find(begin(tracked), end(tracked), pcode) != end(tracked)) return true;
-    return false;
+    return std::find(begin(tracked), end(tracked), pcode) != end(tracked);
+  }
+
+  inline HEPEnergyType ProposalProcessBase::getOptimizedEmCut(Code code) const {
+    // get energy above which energy losses need to be considered
+    auto const production_threshold = get_energy_production_threshold(code);
+
+    HEPEnergyType lowest_table_value = 0_GeV;
+
+    // find tables for EnergyCuts closest (but still smaller than) production_threshold
+    for (auto const& table_energy : energycut_table_values) {
+      if (table_energy <= production_threshold && table_energy > lowest_table_value) {
+        lowest_table_value = table_energy;
+      }
+    }
+
+    if (lowest_table_value == 0_GeV) {
+      // no appropriate table available
+      return production_threshold;
+    }
+
+    return lowest_table_value;
   }
 
-  inline ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env)
-      : RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
+  template <typename TEnvironment>
+  inline ProposalProcessBase::ProposalProcessBase(TEnvironment const& _env) {
     _env.getUniverse()->walk([&](auto& vtn) {
       if (vtn.hasModelProperties()) {
         const auto& prop = vtn.getModelProperties();
-        const auto& medium = mediumData(
-            prop.getMedium(Point(_env.getCoordinateSystem(), 0_cm, 0_cm, 0_cm)));
+        const auto& medium = mediumData(prop.getMedium());
 
-        auto comp_vec = std::vector<PROPOSAL::Components::Component>();
+        auto comp_vec = std::vector<PROPOSAL::Component>();
         const auto& comp = prop.getNuclearComposition();
         auto frac_iter = comp.getFractions().cbegin();
-        for (auto& pcode : comp.getComponents()) {
+        for (auto const pcode : comp.getComponents()) {
           comp_vec.emplace_back(std::string(get_name(pcode)), get_nucleus_Z(pcode),
                                 get_nucleus_A(pcode), *frac_iter);
           ++frac_iter;
@@ -54,21 +69,38 @@ namespace corsika::proposal {
       }
     });
 
-    PROPOSAL::InterpolationDef::order_of_interpolation = 2;
-    PROPOSAL::InterpolationDef::nodes_cross_section = 100;
-    PROPOSAL::InterpolationDef::nodes_propagate = 1000;
-
     //! If corsika data exist store interpolation tables to the corresponding
     //! path, otherwise interpolation tables would only stored in main memory if
     //! no explicit intrpolation def is specified.
-    if (auto data_path = std::getenv("CORSIKA_DATA")) {
-      PROPOSAL::InterpolationDef::path_to_tables = std::string(data_path) + "/PROPOSAL";
-    } else {
-      throw std::runtime_error(
-          "It is not recommended to run PROPOSAL without its tables in "
-          "$CORSIKA_DATA/PROPOSAL. This would be extremely slow. Please provide the "
-          "table directory. ");
+    PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
+
+    //! Initialize EnergyCutSettings
+    for (auto const particle_code : tracked) {
+      if (particle_code == Code::Photon) {
+        // no EnergyCut for photon, only-stochastic propagation
+        continue;
+      }
+      // use optimized emcut for which tables should be available
+      auto optimized_ecut = getOptimizedEmCut(particle_code);
+      CORSIKA_LOG_INFO("PROPOSAL: Use tables with EnergyCut of {} for particle {}.",
+                       optimized_ecut, particle_code);
+      proposal_energycutsettings[particle_code] = optimized_ecut;
     }
+
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto const& medium : media) {
+      for (auto const particle_code : tracked) {
+        buildTables(medium.second, particle_code,
+                    proposal_energycutsettings[particle_code]);
+      }
+    }
+  }
+
+  void ProposalProcessBase::buildTables(PROPOSAL::Medium medium, Code particle_code,
+                                        HEPEnergyType emcut) {
+    CORSIKA_LOG_DEBUG("PROPOSAL: Initialize tables for particle {} in {}", particle_code,
+                      medium.GetName());
+    cross[particle_code](medium, emcut);
   }
 
   inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const

corsika/detail/modules/pythia8/Decay.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/modules/pythia8/Pythia8.hpp>
@@ -12,9 +11,6 @@
 
 #include <corsika/framework/utility/COMBoost.hpp>
 
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
 namespace corsika::pythia8 {
 
   inline Decay::Decay(bool const print_listing)
@@ -36,8 +32,7 @@ namespace corsika::pythia8 {
     // run this only once during construction
 
     // link random number generator in pythia to CORSIKA8
-    Pythia8::RndmEngine* rndm = new corsika::pythia8::Random();
-    Pythia8::Pythia::setRndmEnginePtr(rndm);
+    Pythia8::Pythia::setRndmEnginePtr(std::make_shared<corsika::pythia8::Random>());
 
     Pythia8::Pythia::readString("Next:numberShowInfo = 0");
     Pythia8::Pythia::readString("Next:numberShowProcess = 0");
@@ -62,8 +57,10 @@ namespace corsika::pythia8 {
 
     //    Pythia8::Pythia::particleData.readString("59:m0 = 101.00");
 
+    // LCOV_EXCL_START, we don't validate pythia8 internals
     if (!Pythia8::Pythia::init())
       throw std::runtime_error("Pythia::Decay: Initialization failed!");
+    // LCOV_EXCL_STOP
   }
 
   inline bool Decay::canHandleDecay(Code const vParticleCode) {
@@ -74,7 +71,7 @@ namespace corsika::pythia8 {
         vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
         vParticleCode == Code::Electron || vParticleCode == Code::Positron)
       return false;
-    else if (canDecay(vParticleCode)) // non-zero for particles known to sibyll
+    else if (canDecay(vParticleCode)) // check pythia8 internal
       return true;
     else
       return false;
@@ -91,7 +88,7 @@ namespace corsika::pythia8 {
 
   inline void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
     handleAllDecays_ = false;
-    for (auto p : vParticleList) setHandleDecay(p);
+    for (auto const p : vParticleList) setHandleDecay(p);
   }
 
   inline bool Decay::isDecayHandled(Code const vParticleCode) {
@@ -102,7 +99,7 @@ namespace corsika::pythia8 {
   }
 
   inline void Decay::setStable(std::vector<Code> const& particleList) {
-    for (auto p : particleList) Decay::setStable(p);
+    for (auto const p : particleList) Decay::setStable(p);
   }
 
   inline void Decay::setUnstable(Code const pCode) {
@@ -137,7 +134,7 @@ namespace corsika::pythia8 {
     if (handleAllDecays_)
       CORSIKA_LOG_INFO(" all particles known to Pythia are handled by Pythia::Decay!");
     else
-      for (auto& pCode : handledDecays_)
+      for (auto const pCode : handledDecays_)
         CORSIKA_LOG_INFO("Decay of {} is handled by Pythia!", pCode);
   }
 
@@ -171,11 +168,8 @@ namespace corsika::pythia8 {
 
     auto projectile = view.getProjectile();
 
-    auto const& decayPoint = projectile.getPosition();
-    auto const t0 = projectile.getTime();
-
     auto const& labMomentum = projectile.getMomentum();
-    CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
+    [[maybe_unused]] CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
 
     // define target kinematics in lab frame
     // define boost to and from CoM frame
@@ -204,38 +198,52 @@ namespace corsika::pythia8 {
     double const m = en;
 
     // add particle to pythia stack
-    event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
+    event.append(pdgCode, 1, 0, 0, // PID, status, col, acol
+                 px, py, pz, en, m);
 
+    // LCOV_EXCL_START, we don't validate pythia8 internals
     if (!Pythia8::Pythia::next())
       throw std::runtime_error("Pythia::Decay: decay failed!");
-    else
-      CORSIKA_LOG_DEBUG("Pythia::Decay: particles after decay: {} ", event.size());
+    // LCOV_EXCL_STOP
 
+    CORSIKA_LOG_DEBUG("Pythia::Decay: particles after decay: {} ", event.size());
+
+    // LCOV_EXCL_START, we don't validate pythia8 internals
     if (print_listing_) {
       // list final state
       event.list();
     }
+    // LCOV_EXCL_STOP
 
     // loop over final state
-    for (int i = 0; i < event.size(); ++i)
+    for (int i = 0; i < event.size(); ++i) {
       if (event[i].isFinal()) {
-        auto const pyId = convert_from_PDG(static_cast<PDGCode>(event[i].id()));
-        HEPEnergyType const Erest = event[i].e() * 1_GeV;
-        MomentumVector const pRest(
-            rotatedCS,
-            {event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
-        FourVector const fourMomRest{Erest, pRest};
-        auto const fourMomLab = boost.fromCoM(fourMomRest);
-
-        CORSIKA_LOG_TRACE(
-            "particle: id={} momentum={} energy={} ", pyId,
-            fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
-            fourMomLab.getTimeLikeComponent());
-
-        view.addSecondary(std::make_tuple(pyId, fourMomLab.getTimeLikeComponent(),
-                                          fourMomLab.getSpaceLikeComponents(), decayPoint,
-                                          t0));
+        try {
+          auto const id = static_cast<PDGCode>(event[i].id());
+          auto const pyId = convert_from_PDG(id);
+          HEPEnergyType const Erest = event[i].e() * 1_GeV;
+          MomentumVector const pRest(
+              rotatedCS,
+              {event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
+          FourVector const fourMomRest{Erest, pRest};
+          auto const fourMomLab = boost.fromCoM(fourMomRest);
+          auto const p3 = fourMomLab.getSpaceLikeComponents();
+
+          HEPEnergyType const mass = get_mass(pyId);
+          HEPEnergyType const Ekin = sqrt(p3.getSquaredNorm() + mass * mass) - mass;
+
+          CORSIKA_LOG_TRACE(
+              "particle: id={} momentum={} energy={} ", pyId,
+              fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
+              fourMomLab.getTimeLikeComponent());
+
+          view.addSecondary(std::make_tuple(pyId, Ekin, p3.normalized()));
+        } catch (std::out_of_range const& ex) {
+          CORSIKA_LOG_CRITICAL("Pythia ID {} unknown in C8", event[i].id());
+          throw ex;
+        }
       }
+    }
 
     // set particle stable
     Decay::setStable(particleId);

corsika/detail/modules/pythia8/Interaction.inl

-/*
- * (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/pythia8/Interaction.hpp>
-
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-
-#include <tuple>
-
-namespace corsika::pythia8 {
-
-  inline Interaction::~Interaction() {
-    CORSIKA_LOG_INFO("Pythia::Interaction n= {}", count_);
-  }
-
-  inline Interaction::Interaction(bool const print_listing)
-      : Pythia8::Pythia(CORSIKA_Pythia8_XML_DIR)
-      , print_listing_(print_listing) {
-
-    CORSIKA_LOG_INFO("Configuring Pythia8 from: {}", CORSIKA_Pythia8_XML_DIR);
-
-    // initialize Pythia
-
-    // reduce output from pythia if set to "on"
-    Pythia8::Pythia::readString("Print:quiet = on");
-    // check if data in particle data file is minimally consistent. Very verbose! set to
-    // "off"! we do not change the basic file provided by pythia.
-    Pythia8::Pythia::readString("Check:particleData = off");
-    Pythia8::Pythia::readString("Check:event = on");             // default: on
-    Pythia8::Pythia::readString("Check:levelParticleData = 12"); // 1 is default
-    /** \TODO: proper process initialization for MinBias needed, see
-        also Issue https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/369 **/
-    Pythia8::Pythia::readString("HardQCD:all = on");
-    Pythia8::Pythia::readString("ProcessLevel:resonanceDecays = off");
-
-    if (!Pythia8::Pythia::init())
-      throw std::runtime_error("Pythia::Interaction: Initialization failed!");
-
-    // any decays in pythia? if yes need to define which particles
-    if (internalDecays_) {
-      // define which particles are passed to corsika, i.e. which particles make it into
-      // history even very shortlived particles like charm or pi0 are of interest here
-      std::vector<Code> const HadronsWeWantTrackedByCorsika = {
-          Code::PiPlus, Code::PiMinus, Code::Pi0,        Code::KMinus,     Code::KPlus,
-          Code::K0Long, Code::K0Short, Code::SigmaPlus,  Code::SigmaMinus, Code::Lambda0,
-          Code::Xi0,    Code::XiMinus, Code::OmegaMinus, Code::DPlus,      Code::DMinus,
-          Code::D0,     Code::D0Bar};
-
-      Interaction::setStable(HadronsWeWantTrackedByCorsika);
-    }
-
-    // basic initialization of cross section routines
-    sigma_.init(&(Pythia8::Pythia::info), Pythia8::Pythia::settings,
-                &(Pythia8::Pythia::particleData), &(Pythia8::Pythia::rndm));
-  }
-
-  inline void Interaction::setStable(std::vector<Code> const& particleList) {
-    for (auto p : particleList) Interaction::setStable(p);
-  }
-
-  inline void Interaction::setUnstable(Code const pCode) {
-    CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} unstable..", pCode);
-    Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
-  }
-
-  inline void Interaction::setStable(Code const pCode) {
-    CORSIKA_LOG_DEBUG("Pythia::Interaction: setting {} stable..", pCode);
-    Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
-  }
-
-  inline void Interaction::configureLabFrameCollision(Code const BeamId,
-                                                      Code const TargetId,
-                                                      HEPEnergyType const BeamEnergy) {
-    // Pythia configuration of the current event
-    // very clumsy. I am sure this can be done better..
-
-    // set beam
-    // beam id for pythia
-    auto const pdgBeam = static_cast<int>(get_PDG(BeamId));
-    std::stringstream stBeam;
-    stBeam << "Beams:idA = " << pdgBeam;
-    Pythia8::Pythia::readString(stBeam.str());
-    // set target
-    auto pdgTarget = static_cast<int>(get_PDG(TargetId));
-    // replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
-    if (TargetId == Code::Hydrogen) pdgTarget = static_cast<int>(get_PDG(Code::Proton));
-    std::stringstream stTarget;
-    stTarget << "Beams:idB = " << pdgTarget;
-    Pythia8::Pythia::readString(stTarget.str());
-    // set frame to lab. frame
-    Pythia8::Pythia::readString("Beams:frameType = 2");
-    // set beam energy
-    double const Elab = BeamEnergy / 1_GeV;
-    std::stringstream stEnergy;
-    stEnergy << "Beams:eA = " << Elab;
-    Pythia8::Pythia::readString(stEnergy.str());
-    // target at rest
-    Pythia8::Pythia::readString("Beams:eB = 0.");
-    // initialize this config
-
-    if (!Pythia8::Pythia::init())
-      throw std::runtime_error("Pythia::Interaction: Initialization failed!");
-  }
-
-  inline bool Interaction::canInteract(Code const pCode) {
-    return pCode == Code::Proton || pCode == Code::Neutron || pCode == Code::AntiProton ||
-           pCode == Code::AntiNeutron || pCode == Code::PiMinus || pCode == Code::PiPlus;
-  }
-
-  inline std::tuple<CrossSectionType, CrossSectionType> Interaction::getCrossSection(
-      Code const BeamId, Code const TargetId, HEPEnergyType const CoMenergy) {
-    // interaction possible in pythia?
-    if (TargetId == Code::Proton || TargetId == Code::Hydrogen) {
-      if (canInteract(BeamId) && isValidCoMEnergy(CoMenergy)) {
-        // input particle PDG
-        auto const pdgCodeBeam = static_cast<int>(get_PDG(BeamId));
-        auto const pdgCodeTarget = static_cast<int>(get_PDG(TargetId));
-        double const ecm = CoMenergy / 1_GeV;
-
-        // calculate cross section
-        sigma_.calc(pdgCodeBeam, pdgCodeTarget, ecm);
-        if (sigma_.hasSigmaTot()) {
-          double const sigEla = sigma_.sigmaEl();
-          double const sigProd = sigma_.sigmaTot() - sigEla;
-
-          return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
-
-        } else
-          throw std::runtime_error("pythia cross section init failed");
-
-      } else {
-        return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
-                               std::numeric_limits<double>::infinity() * 1_mb);
-      }
-    } else {
-      throw std::runtime_error("invalid target for pythia");
-    }
-  }
-
-  inline GrammageType Interaction::getInteractionLength(
-      corsika::setup::Stack::particle_type const& particle) {
-
-    // coordinate system, get global frame of reference
-    MomentumVector const& pMomentum = particle.getMomentum();
-    CoordinateSystemPtr const& labCS = pMomentum.getCoordinateSystem();
-
-    Code const corsikaBeamId = particle.getPID();
-
-    // beam particles for pythia : 1, 2, 3 for p, pi, k
-    // read from cross section code table
-    bool const kInteraction = canInteract(corsikaBeamId);
-
-    // FOR NOW: assume target is at rest
-    MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV});
-
-    // total momentum and energy
-    HEPEnergyType Elab = particle.getEnergy() + constants::nucleonMass;
-    MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
-    pTotLab += pMomentum;
-    pTotLab += pTarget;
-    auto const pTotLabNorm = pTotLab.getNorm();
-    // calculate cm. energy
-    HEPEnergyType const ECoM = sqrt(
-        (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
-
-    CORSIKA_LOG_DEBUG(
-        "Interaction: LambdaInt: \n"
-        " input energy: {} GeV"
-        " beam can interact: {}"
-        " beam pid: {}",
-        particle.getEnergy() / 1_GeV, kInteraction, particle.getPID());
-
-    // TODO: move limits into variables
-    if (kInteraction && Elab >= 8.5_GeV && isValidCoMEnergy(ECoM)) {
-
-      // 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();
-      auto const mediumComposition =
-          currentNode->getModelProperties().getNuclearComposition();
-      // determine average interaction length
-
-      auto const weightedProdCrossSection =
-          mediumComposition.getWeightedSum([=](auto vTargetID) {
-            return std::get<0>(this->getCrossSection(corsikaBeamId, vTargetID, ECoM));
-          });
-
-      CORSIKA_LOG_DEBUG(
-          "Interaction: IntLength: weighted CrossSection (mb): {} "
-          "Interaction: IntLength: average mass number: {} ",
-          weightedProdCrossSection / 1_mb, mediumComposition.getAverageMassNumber());
-
-      // 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);
-  }
-
-  template <class TView>
-  inline void Interaction::doInteraction(TView& view) {
-
-    auto projectile = view.getProjectile();
-
-    const auto corsikaBeamId = projectile.getPID();
-    CORSIKA_LOG_DEBUG(
-        "Pythia::Interaction: "
-        "DoInteraction: {} interaction? ",
-        corsikaBeamId, corsika::pythia8::Interaction::canInteract(corsikaBeamId));
-
-    if (is_nucleus(corsikaBeamId)) {
-      // nuclei handled by different process, this should not happen
-      throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
-    }
-
-    if (corsika::pythia8::Interaction::canInteract(corsikaBeamId)) {
-
-      // define projectile
-      HEPEnergyType const eProjectileLab = projectile.getEnergy();
-      auto const pProjectileLab = projectile.getMomentum();
-      CoordinateSystemPtr const& labCS = pProjectileLab.getCoordinateSystem();
-
-      // position and time of interaction, not used in Sibyll
-      Point pOrig = projectile.getPosition();
-      TimeType tOrig = projectile.getTime();
-
-      // define target
-      // FOR NOW: target is always at rest
-      auto const eTargetLab = 0_GeV + constants::nucleonMass;
-      auto const pTargetLab = MomentumVector(labCS, 0_GeV, 0_GeV, 0_GeV);
-      FourVector const PtargLab(eTargetLab, pTargetLab);
-
-      CORSIKA_LOG_DEBUG(
-          "Interaction: ebeam lab: {} GeV"
-          "Interaction: pbeam lab: {} GeV",
-          eProjectileLab / 1_GeV, pProjectileLab.getComponents() / 1_GeV);
-
-      CORSIKA_LOG_DEBUG(
-          "Interaction: etarget lab: {} GeV"
-          "Interaction: ptarget lab: {} GeV ",
-          eTargetLab / 1_GeV, pTargetLab.getComponents() / 1_GeV);
-
-      FourVector const PprojLab(eProjectileLab, pProjectileLab);
-
-      // define target kinematics in lab frame
-      // define boost to and from CoM frame
-      // CoM frame definition in Pythia projectile: +z
-      COMBoost const boost(PprojLab, constants::nucleonMass);
-
-      // just for show:
-      // boost projecticle
-      auto const PprojCoM = boost.toCoM(PprojLab);
-
-      // boost target
-      auto const PtargCoM = boost.toCoM(PtargLab);
-
-      CORSIKA_LOG_DEBUG(
-          "Interaction: ebeam CoM: {} GeV"
-          "Interaction: pbeam CoM: {} GeV",
-          PprojCoM.getTimeLikeComponent() / 1_GeV,
-          PprojCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
-
-      CORSIKA_LOG_DEBUG(
-          "Interaction: etarget CoM: {} GeV"
-          "Interaction: ptarget CoM: {} GeV",
-          PtargCoM.getTimeLikeComponent() / 1_GeV,
-          PtargCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
-
-      CORSIKA_LOG_DEBUG("Interaction: position of interaction: ", pOrig.getCoordinates());
-      CORSIKA_LOG_DEBUG("Interaction: time: {}", tOrig);
-
-      HEPEnergyType Etot = eProjectileLab + eTargetLab;
-      MomentumVector Ptot = projectile.getMomentum();
-      // invariant mass, i.e. cm. energy
-      HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.getSquaredNorm());
-
-      // sample target mass number
-      auto const* currentNode = projectile.getNode();
-      auto const& mediumComposition =
-          currentNode->getModelProperties().getNuclearComposition();
-      // get cross sections for target materials
-      /*
-        Here we read the cross section from the interaction model again,
-        should be passed from getInteractionLength if possible
-       */
-      //#warning reading interaction cross section again, should not be necessary
-      auto const& compVec = mediumComposition.getComponents();
-      std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
-
-      for (size_t i = 0; i < compVec.size(); ++i) {
-        auto const targetId = compVec[i];
-        auto const [sigProd, sigEla] = getCrossSection(corsikaBeamId, targetId, Ecm);
-        [[maybe_unused]] auto const& dummy_sigEla = sigEla;
-        cross_section_of_components[i] = sigProd;
-      }
-
-      auto const corsikaTargetId =
-          mediumComposition.sampleTarget(cross_section_of_components, RNG_);
-      CORSIKA_LOG_DEBUG("Interaction: target selected: {}", corsikaTargetId);
-
-      if (corsikaTargetId != Code::Hydrogen && corsikaTargetId != Code::Neutron &&
-          corsikaTargetId != Code::Proton)
-        throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
-
-      CORSIKA_LOG_DEBUG(
-          "Interaction: "
-          " DoInteraction: E(GeV): {}"
-          " Ecm(GeV): {}",
-          eProjectileLab / 1_GeV, Ecm / 1_GeV);
-
-      if (eProjectileLab < 8.5_GeV || !isValidCoMEnergy(Ecm)) {
-        CORSIKA_LOG_DEBUG(
-            "Interaction: "
-            " DoInteraction: should have dropped particle.. "
-            "THIS IS AN ERROR");
-        throw std::runtime_error("energy too low for PYTHIA");
-
-      } else {
-        count_++;
-
-        configureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
-
-        // create event in pytia
-        if (!Pythia8::Pythia::next())
-          throw std::runtime_error("Pythia::DoInteraction: failed!");
-
-        // link to pythia stack
-        Pythia8::Event& event = Pythia8::Pythia::event;
-
-        if (print_listing_) {
-          // print final state
-          event.list();
-        }
-
-        MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
-        HEPEnergyType Elab_final = 0_GeV;
-        for (int i = 0; i < event.size(); ++i) {
-          Pythia8::Particle& p8p = event[i];
-          // skip particles that have decayed in pythia
-          if (!p8p.isFinal()) continue;
-
-          auto const pyId = convert_from_PDG(static_cast<PDGCode>(p8p.id()));
-
-          MomentumVector const pyPlab(
-              labCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
-          HEPEnergyType const pyEn = p8p.e() * 1_GeV;
-
-          // add to corsika stack
-          auto pnew =
-              projectile.addSecondary(std::make_tuple(pyId, pyEn, pyPlab, pOrig, tOrig));
-
-          Plab_final += pnew.getMomentum();
-          Elab_final += pnew.getEnergy();
-        }
-        CORSIKA_LOG_DEBUG(
-            "conservation (all GeV): "
-            "Elab_final= {}"
-            ", Plab_final= {}",
-            Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
-      }
-    }
-  }
-
-} // namespace corsika::pythia8

corsika/detail/modules/pythia8/InteractionModel.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <utility>
+#include <stdexcept>
+
+#include <boost/filesystem/path.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+#include <fmt/core.h>
+#include <cnpy.hpp>
+
+#include <corsika/modules/pythia8/InteractionModel.hpp>
+
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/utility/CrossSectionTable.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+namespace corsika::pythia8 {
+
+  inline InteractionModel::~InteractionModel() {}
+
+  inline InteractionModel::InteractionModel(std::set<Code> const& stableParticles,
+                                            boost::filesystem::path const& dataPath,
+                                            bool const printListing)
+      : crossSectionPPElastic_{loadPPTable(dataPath, "sig_el")}
+      , crossSectionPPInelastic_{loadPPTable(dataPath, "sig_inel")}
+      , print_listing_(printListing) {
+    auto rndm = std::make_shared<corsika::pythia8::Random>();
+    pythia_.setRndmEnginePtr(rndm);
+
+    CORSIKA_LOG_INFO("Pythia8 reuse files: {}", dataPath.native());
+
+    // projectile and target init to p Nitrogen to initialize pythia with angantyr
+    pythia_.readString("Beams:allowIDASwitch = on");
+    pythia_.settings.mode("Beams:idA",
+                          static_cast<PDGCodeIntType>(get_PDG(Code::Proton)));
+    pythia_.settings.mode("Beams:idB",
+                          static_cast<PDGCodeIntType>(get_PDG(Code::Oxygen)));
+
+    pythia_.readString("Beams:allowVariableEnergy = on");
+    pythia_.readString("Beams:frameType = 1"); // CoM
+    // initialize a maximum energy event
+    pythia_.settings.parm("Beams:eCM", 400_TeV / 1_GeV);
+
+    // needed when regular Pythia (not Angantyr) takes over for pp
+    pythia_.readString("SoftQCD:all = on");
+
+    /* reuse file settings
+    reuseInit = 2: initialization data is reuse file, but if the file is not found,
+    initialization fails reuseInit = 3: same, but if the file is not found, it will be
+    generated and saved after normal initialization */
+    pythia_.readString("MultipartonInteractions:reuseInit = 2");
+    pythia_.readFile(dataPath.native() + "/pythia8313.cmnd");
+    // switch on decays for all hadrons except the ones defined as tracked by C8
+    if (!stableParticles.empty()) {
+      pythia_.readString("HadronLevel:Decay = on");
+      for (auto pCode : stableParticles)
+        pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
+    } else {
+      // all hadrons stable
+      pythia_.readString("HadronLevel:Decay = on");
+    }
+    // Reduce printout and relax energy-momentum conservation.
+    pythia_.readString("Print:quiet = on");
+    pythia_.readString("Check:epTolErr = 0.1");
+    pythia_.readString("Check:epTolWarn = 0.0001");
+    pythia_.readString("Check:mTolErr = 0.01");
+
+    // Reduce statistics printout to relevant ones.
+    pythia_.readString("Stat:showProcessLevel = off");
+    pythia_.readString("Stat:showPartonLevel = off");
+
+    // initialize
+    // we can't test this block, LCOV_EXCL_START
+    if (!pythia_.init())
+      throw std::runtime_error("Pythia::InteractionModel: Initialization failed!");
+    // LCOV_EXCL_STOP
+
+    // create/load tables of cross sections
+    for (Code const projectile : validProjectiles_) {
+      for (Code const target : validTargets_) {
+        if (xs_map_.find(projectile) == xs_map_.end()) {
+          // projectile not mapped, load table directly
+          auto const tablePath =
+              dataPath / "pythia8_xs_npz" /
+              fmt::format("xs_{}_{}.npz",
+                          static_cast<PDGCodeIntType>(get_PDG(projectile)),
+                          static_cast<PDGCodeIntType>(get_PDG(target)));
+          auto const tables = cnpy::npz_load(tablePath.native());
+          // NOTE: tables are calculated for plab. In C8 we we assume elab. starting at
+          // plab=100GeV the difference is at most (1+m**2/p**2)=1.000088035
+          auto const energies = tables.at("plab").as_vec<float>();
+          auto const total_xs = tables.at("sig_tot").as_vec<float>();
+
+          if (auto const e_size = energies.size(), xs_size = total_xs.size();
+              xs_size != e_size) {
+            throw std::runtime_error{
+                fmt::format("Pythia8 table corrupt (plab size = {}; sig_tot size = {})",
+                            e_size, xs_size)};
+          }
+
+          auto xs_it = boost::make_transform_iterator(total_xs.cbegin(), millibarn_mult);
+          auto e_it_beg = boost::make_transform_iterator(energies.cbegin(), GeV_mult);
+          decltype(e_it_beg) e_it_end =
+              boost::make_transform_iterator(energies.cend(), GeV_mult);
+          crossSectionTables_.insert(
+              {std::pair{projectile, target},
+               CrossSectionTable<InterpolationTransforms::Log>{
+                   std::move(e_it_beg), std::move(e_it_end), std::move(xs_it)}});
+        }
+      }
+    }
+  }
+
+  inline CrossSectionTable<InterpolationTransforms::Log> InteractionModel::loadPPTable(
+      boost::filesystem::path const& dataPath, char const* key) {
+    auto const ppTablePath =
+        dataPath / "pythia8_xs_npz" /
+        fmt::format("xs_{}_{}.npz", static_cast<PDGCodeIntType>(get_PDG(Code::Proton)),
+                    static_cast<PDGCodeIntType>(get_PDG(Code::Proton)));
+    auto const tables = cnpy::npz_load(ppTablePath.native());
+    // NOTE: tables are calculated for plab. In C8 we we assume elab. starting at
+    // plab=100GeV the difference is at most (1+m**2/p**2)=1.000088035
+    auto const energies = tables.at("plab").as_vec<float>();
+    auto const xs = tables.at(key).as_vec<float>();
+
+    if (auto e_size = energies.size(), xs_size = xs.size(); xs_size != e_size) {
+      throw std::runtime_error{fmt::format(
+          "Pythia8 table corrupt (plab size = {}; xs size = {})", e_size, xs_size)};
+    }
+
+    auto xs_it = boost::make_transform_iterator(xs.cbegin(), millibarn_mult);
+    auto e_it_beg = boost::make_transform_iterator(energies.cbegin(), GeV_mult);
+    decltype(e_it_beg) e_it_end =
+        boost::make_transform_iterator(energies.cend(), GeV_mult);
+    return CrossSectionTable<InterpolationTransforms::Log>{
+        std::move(e_it_beg), std::move(e_it_end), std::move(xs_it)};
+  }
+
+  inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+                                        HEPEnergyType const sqrtSNN) const {
+
+    CORSIKA_LOG_DEBUG("pythia isValid: {} + {} at sqrtSNN = {} GeV", projectileId,
+                      targetId, sqrtSNN / 1_GeV);
+    if (is_nucleus(targetId))
+      CORSIKA_LOG_DEBUG("nucleus = {}-{}", get_nucleus_A(targetId),
+                        get_nucleus_Z(targetId));
+    if (is_nucleus(projectileId)) // not yet possible with Pythia
+      return false;
+
+    if (!canInteract(projectileId)) return false;
+
+    auto const mass_target =
+        get_mass(targetId) / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.);
+    HEPEnergyType const labE =
+        calculate_lab_energy(static_pow<2>(sqrtSNN), get_mass(projectileId), mass_target);
+    if (labE < eKinMinLab_) return false;
+
+    bool const validProjectile =
+        std::find(validProjectiles_.begin(), validProjectiles_.end(), projectileId) !=
+        validProjectiles_.end();
+    bool const validTarget = std::find(validTargets_.begin(), validTargets_.end(),
+                                       targetId) != validTargets_.end();
+    return validProjectile && validTarget;
+  }
+
+  inline bool InteractionModel::canInteract(Code const pCode) const {
+    return is_hadron(pCode) && !is_nucleus(pCode);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
+                                           FourMomentum const& projectileP4,
+                                           FourMomentum const& targetP4) const {
+    // elastic cross sections only available for proton
+    if (projectileId != Code::Proton || targetId != Code::Proton) {
+      return std::tuple{CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+
+    // NOTE: here we calculate SNN (per nucleon energy) AND S (per particle energy)
+    // because isValid expects sqrtSNN as input but the tables need Elab
+    auto const Aprojectile =
+        (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1.);
+    auto const Atarget = (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.);
+    auto const SNN = (projectileP4 / Aprojectile + targetP4 / Atarget).getNormSqr();
+    HEPEnergyType const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN))
+      return std::tuple{CrossSectionType::zero(), CrossSectionType::zero()};
+
+    // read cross section from table
+    // define system (the tables were created for lab. energies) since we cannot assume
+    // the input is in the lab. frame we calculate the lab. energy from SNN =
+    // (projectileP4 / Aprojectile + targetP4 / Atarget)**2 and Elab = S/(2. * mTarget /
+    // Atarget) where A* is the number of nucleons in a nucleus. A* is 1 if projectile or
+    // target are simple hadrons.
+    HEPEnergyType const Elab = calculate_lab_energy(
+        SNN, get_mass(projectileId) / Aprojectile, get_mass(targetId) / Atarget);
+
+    CrossSectionType const xsEl = crossSectionPPElastic_.interpolate(Elab);
+    CrossSectionType const xsInel = crossSectionPPInelastic_.interpolate(Elab);
+    return std::pair{xsInel, xsEl};
+  }
+
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+
+    auto const Aprojectile =
+        (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1.);
+    auto const Atarget = (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1.);
+
+    auto const SNN = (projectileP4 / Aprojectile + targetP4 / Atarget).getNormSqr();
+    HEPEnergyType const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN)) return CrossSectionType::zero();
+
+    // read cross section from table
+    // define system (the tables were created for lab. energies) since we cannot assume
+    // the input is in the lab. frame we calculate the lab. energy from SNN =
+    // (projectileP4 / Aprojectile + targetP4 / Atarget)**2 and Elab = S/(2. * mTarget /
+    // Atarget) where A* is the number of nucleons in a nucleus. A* is 1 if projectile or
+    // target are simple hadrons.
+    auto const it = xs_map_.find(projectileId);
+    auto const mapped_projectile = (it == xs_map_.end()) ? projectileId : it->second;
+    auto const& table = crossSectionTables_.at(std::pair{mapped_projectile, targetId});
+    HEPEnergyType const Elab = calculate_lab_energy(
+        SNN, get_mass(projectileId) / Aprojectile, get_mass(targetId) / Atarget);
+    CORSIKA_LOG_DEBUG("pythia getCrossSection: {}+{} at Elab= {} GeV, sqrtSNN = {} GeV",
+                      projectileId, get_name(targetId), Elab / 1_GeV, sqrtSNN / 1_GeV);
+    return table.interpolate(Elab);
+  }
+
+  template <class TView>
+  inline void InteractionModel::doInteraction(TView& view, Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    CORSIKA_LOG_DEBUG(
+        "Pythia::InteractionModel: "
+        "doInteraction: {} interaction? {} ",
+        projectileId, corsika::pythia8::InteractionModel::canInteract(projectileId));
+
+    // define system nucleon-nucleon
+    auto const projectileP4NN =
+        projectileP4 / (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1);
+    auto const targetP4NN =
+        targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+    auto const SNN = (projectileP4NN + targetP4NN).getNormSqr();
+    auto const sqrtSNN = sqrt(SNN);
+    HEPEnergyType const eProjectileLab = calculate_lab_energy(
+        SNN, get_mass(projectileId), get_mass(targetId) / get_nucleus_A(targetId));
+
+    CORSIKA_LOG_DEBUG("InteractionModel: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
+    CORSIKA_LOG_DEBUG(
+        "InteractionModel: "
+        " doInteraction: E(GeV): {}"
+        " Ecm_NN(GeV): {}",
+        eProjectileLab / 1_GeV, sqrtSNN / 1_GeV);
+
+    if (!isValid(projectileId, targetId, sqrtSNN)) {
+      throw std::runtime_error("invalid target,projectile,energy combination.");
+    }
+
+    COMBoost const boost{projectileP4NN, targetP4NN};
+    auto const& rotCS = boost.getRotatedCS();
+
+    // variables to talk to Pythia
+    double const eCM_pythia = sqrtSNN / 1_GeV; // CoM energy hadron-Nucleon
+    double const idA_pythia = static_cast<double>(get_PDG(projectileId));
+    double const idB_pythia = static_cast<double>(get_PDG(targetId));
+    CORSIKA_LOG_DEBUG("re-configuring pythia idA={}, idB={}, ecm={} GeV", idA_pythia,
+                      idB_pythia, eCM_pythia);
+    // configure event
+    pythia_.setBeamIDs(idA_pythia, idB_pythia);
+    pythia_.setKinematics(eCM_pythia);
+
+    // generate event (one chance only!)
+    if (!pythia_.next()) {
+      throw std::runtime_error("Pythia::InteractionModel: interaction failed!");
+    }
+
+    // LCOV_EXCL_START, we don't validate pythia8 internals
+    if (print_listing_) {
+      // print final state
+      pythia_.event.list();
+    }
+    // LCOV_EXCL_STOP
+
+    MomentumVector Plab_final{boost.getOriginalCS()};
+    auto Elab_final = HEPEnergyType::zero();
+
+    for (Pythia8::Particle const& p8p : pythia_.event) {
+      // skip particles that have decayed and initial particles
+      if (!p8p.isFinal()) continue;
+      try {
+        auto const volatile id = static_cast<PDGCode>(p8p.id());
+        auto const pyId = convert_from_PDG(id);
+
+        // skip nuclear remnants
+        if (is_nucleus(pyId)) continue;
+
+        MomentumVector const pyPcom(
+            rotCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+        auto const pyP = boost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPcom});
+
+        HEPEnergyType const mass = get_mass(pyId);
+        HEPEnergyType const Ekin =
+            sqrt(pyP.getSpaceLikeComponents().getSquaredNorm() + mass * mass) - mass;
+
+        // add to corsika stack
+        auto pnew = view.addSecondary(
+            std::make_tuple(pyId, Ekin, pyP.getSpaceLikeComponents().normalized()));
+
+        Plab_final += pnew.getMomentum();
+        Elab_final += pnew.getEnergy();
+      }
+      // irreproducible in tests, LCOV_EXCL_START
+      catch (std::out_of_range const& ex) {
+        CORSIKA_LOG_CRITICAL("Pythia ID {} unknown in C8", p8p.id());
+        throw ex;
+      }
+      // LCOV_EXCL_STOP
+    }
+
+    pythia_.event.clear();
+
+    CORSIKA_LOG_DEBUG(
+        "conservation (all GeV): "
+        "Elab_final= {}"
+        ", Plab_final= {}",
+        Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
+  }
+
+} // namespace corsika::pythia8

corsika/detail/modules/pythia8/NeutrinoInteraction.inl

+/*
+ * (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <boost/filesystem/path.hpp>
+#include "corsika/framework/core/EnergyMomentumOperations.hpp"
+
+namespace corsika::pythia8 {
+
+  inline NeutrinoInteraction::NeutrinoInteraction(std::set<Code> const& list,
+                                                  bool const& handleNC,
+                                                  bool const& handleCC)
+      : stable_particles_(list)
+      , handle_nc_(handleNC)
+      , handle_cc_(handleCC)
+      , pythiaMain_{CORSIKA_Pythia8_XML_DIR, false} {
+
+    pythiaMain_.setRndmEnginePtr(std::make_shared<corsika::pythia8::Random>());
+  }
+
+  inline NeutrinoInteraction::~NeutrinoInteraction() {
+    CORSIKA_LOGGER_DEBUG(logger_, "Pythia::NeutrinoInteraction n= {}", count_);
+  }
+
+  template <class TView>
+  void NeutrinoInteraction::doInteraction(TView& view, Code const projectileId,
+                                          Code const targetId,
+                                          FourMomentum const& projectileP4,
+                                          FourMomentum const& targetP4) {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "Primary {} - {} interaction with E_nu = {} GeV",
+                         projectileId, targetId,
+                         projectileP4.getTimeLikeComponent() / 1_GeV);
+    CORSIKA_LOGGER_DEBUG(
+        logger_, "configure Pythia for primary neutrino interactions. NC={}, CC={}",
+        handle_nc_, handle_cc_);
+    if (!handle_nc_ && !handle_cc_) {
+      CORSIKA_LOGGER_ERROR(
+          logger_,
+          "no neutrino interaction channel configured! Select either NC, CC or both!");
+      throw std::runtime_error("Configuration error!");
+    }
+    CORSIKA_LOGGER_DEBUG(logger_, "minimal Q2 in DIS: {} GeV2", minQ2_ / 1_GeV / 1_GeV);
+
+    if (!isValid(projectileId, targetId, projectileP4, targetP4)) {
+      CORSIKA_LOGGER_ERROR(logger_, "wrong projectile, target or energy configuration!");
+      throw std::runtime_error("Configuration error!");
+    }
+    // sample nucleon from nucleus A,Z
+    double const fProtons = get_nucleus_Z(targetId) / double(get_nucleus_A(targetId));
+    double const fNeutrons = 1. - fProtons;
+    std::discrete_distribution<int> nucleonChannelDist{fProtons, fNeutrons};
+    corsika::default_prng_type& rng =
+        corsika::RNGManager<>::getInstance().getRandomStream("pythia");
+    Code const targetNucleonId = (nucleonChannelDist(rng) ? Code::Neutron : Code::Proton);
+    int const idTarget_pythia = static_cast<int>(get_PDG(targetNucleonId));
+    auto const targetP4NN =
+        targetP4 / (is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+    CORSIKA_LOGGER_DEBUG(logger_, "selected {} target. P4NN={}", targetNucleonId,
+                         targetP4NN.getSpaceLikeComponents());
+
+    // set projectile
+    double const Q2min = minQ2_ / 1_GeV / 1_GeV;
+    int const idProjectile_pythia = static_cast<int>(get_PDG(projectileId));
+
+    // calculate CoM energy of the neutrino nucleon interaction
+    double const ecm_pythia =
+        calculate_com_energy(projectileP4.getTimeLikeComponent(), get_mass(projectileId),
+                             get_mass(targetNucleonId)) /
+        1_GeV;
+
+    CORSIKA_LOGGER_DEBUG(logger_, "center-of-mass energy: {} GeV", ecm_pythia);
+    // Set up CoM interaction
+    pythiaMain_.readString("Beams:frameType = 1");
+    // center-of-mass energy
+    pythiaMain_.settings.parm("Beams:eCM", ecm_pythia);
+    // BeamA (+z in CoM)
+    pythiaMain_.settings.mode("Beams:idA", idProjectile_pythia);
+    // BeamB (-z direction in CoM)
+    pythiaMain_.settings.mode("Beams:idB", idTarget_pythia);
+
+    // Set up DIS process within some phase space.
+    // Neutral current (with gamma/Z interference).
+    if (handle_nc_) pythiaMain_.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
+    // Charged current.
+    if (handle_cc_) pythiaMain_.readString("WeakBosonExchange:ff2ff(t:W) = on");
+    // Phase-space cut: minimal Q2 of process.
+    pythiaMain_.settings.parm("PhaseSpace:Q2Min", Q2min);
+
+    // Set dipole recoil on. Necessary for DIS + shower.
+    pythiaMain_.readString("SpaceShower:dipoleRecoil = on");
+
+    // Allow emissions up to the kinematical limit,
+    // since rate known to match well to matrix elements everywhere.
+    pythiaMain_.readString("SpaceShower:pTmaxMatch = 2");
+
+    // QED radiation off lepton not handled yet by the new procedure.
+    pythiaMain_.readString("PDF:lepton = off");
+    pythiaMain_.readString("TimeShower:QEDshowerByL = off");
+
+    // switch on decays for all hadrons except the ones defined as tracked by C8
+    if (!stable_particles_.empty()) {
+      pythiaMain_.readString("HadronLevel:Decay = on");
+      for (auto pCode : stable_particles_)
+        pythiaMain_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
+    } else {
+      // all hadrons stable
+      pythiaMain_.readString("HadronLevel:Decay = off");
+    }
+
+    pythiaMain_.readString("Stat:showProcessLevel = off");
+    pythiaMain_.readString("Stat:showPartonLevel = off");
+
+    pythiaMain_.readString("Print:quiet = on");
+    pythiaMain_.readString("Check:epTolErr = 0.1");
+    pythiaMain_.readString("Check:epTolWarn = 0.0001");
+    pythiaMain_.readString("Check:mTolErr = 0.01");
+
+    pythiaMain_.init();
+
+    // References to the event record
+    Pythia8::Event& eventMain = pythiaMain_.event;
+    // define boost assuming target nucleon is at rest
+    COMBoost const boost{projectileP4, targetP4NN};
+    // the boost is along the total momentum axis and does not include a rotation
+    // get rotated frame where momentum of projectile in the CoM is along +z
+    auto const& rotCS = boost.getRotatedCS();
+
+    if (!pythiaMain_.next()) {
+      throw std::runtime_error("Pythia neutrino collision failed ");
+    } else {
+      CORSIKA_LOGGER_DEBUG(logger_, "pythia neutrino interaction done!");
+    }
+
+    MomentumVector Plab_final{boost.getOriginalCS()};
+    auto Elab_final = HEPEnergyType::zero();
+    CORSIKA_LOGGER_DEBUG(logger_, "particles generated in neutrino interaction:");
+    for (int i = 0; i < eventMain.size(); ++i) {
+      auto const& p8p = eventMain[i];
+      if (p8p.isFinal()) {
+        try {
+          // get particle ids from pythia stack and convert to corsika id
+          auto const volatile id = static_cast<PDGCode>(p8p.id());
+          auto const pyId = convert_from_PDG(id);
+          // get pythia momentum in CoM (rotated cs!)
+          MomentumVector const pyPcom(
+              rotCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+          // boost pythia 4Momentum in CoM to lab. frame. This includes the rotation.
+          // calculate 3-momentum and kinetic energy
+          CORSIKA_LOGGER_DEBUG(logger_, "momentum in CoM: {} GeV",
+                               pyPcom.getComponents() / 1_GeV);
+          auto const pyP4lab = boost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPcom});
+          auto const pyPlab = pyP4lab.getSpaceLikeComponents();
+          HEPEnergyType const pyEkinlab =
+              calculate_kinetic_energy(pyPlab.getNorm(), get_mass(pyId));
+
+          // add to corsika stack
+          auto pnew =
+              view.addSecondary(std::make_tuple(pyId, pyEkinlab, pyPlab.normalized()));
+
+          CORSIKA_LOGGER_DEBUG(logger_, "id = {}, E = {} GeV, p = {} GeV", pyId,
+                               pyEkinlab / 1_GeV, pyPlab.getComponents() / 1_GeV);
+          Plab_final += pnew.getMomentum();
+          Elab_final += pnew.getEnergy();
+        }
+        // irreproducible in tests, LCOV_EXCL_START
+        catch (std::out_of_range const& ex) {
+          CORSIKA_LOGGER_CRITICAL(logger_, "Pythia ID {} unknown in C8", p8p.id());
+          throw ex;
+        }
+        // LCOV_EXCL_STOP
+      }
+    }
+
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "conservation (all GeV): "
+                         "Elab_final= {}"
+                         ", Plab_final= {}",
+                         Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
+
+    count_++;
+  }
+
+} // namespace corsika::pythia8
\ No newline at end of file

corsika/detail/modules/pythia8/Random.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once

corsika/detail/modules/qgsjetII/Interaction.inl

-/*
- * (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/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-
-#include <sstream>
-#include <string>
-#include <tuple>
-
-#include <qgsjet-II-04.hpp>
-
-namespace corsika::qgsjetII {
-
-  inline Interaction::Interaction(const std::string& dataPath)
-      : data_path_(dataPath) {
-    if (dataPath == "") {
-      if (std::getenv("CORSIKA_DATA")) {
-        data_path_ = std::string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
-        CORSIKA_LOG_DEBUG("Searching for QGSJetII data tables in {}", data_path_);
-      }
-    }
-
-    // initialize QgsjetII
-    static bool initialized = false;
-    if (!initialized) {
-      qgset_();
-      datadir DIR(data_path_);
-      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 > 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 > 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& vP) const {
-
-    // coordinate system, get global frame of reference
-    CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
-
-    const Code corsikaBeamId = vP.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 Elab = vP.getEnergy();
-
-    CORSIKA_LOG_DEBUG(
-        "Interaction: LambdaInt: \n"
-        " input energy: {} GeV"
-        " beam can interact: {}"
-        " beam pid: {}",
-        vP.getEnergy() / 1_GeV, kInteraction, vP.getPID());
-
-    if (kInteraction) {
-
-      int Abeam = 0;
-      if (is_nucleus(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();
-
-      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();
-
-    const auto corsikaBeamId = projectile.getPID();
-    CORSIKA_LOG_DEBUG(
-        "ProcessQgsjetII: "
-        "DoInteraction: {} interaction possible? {}",
-        corsikaBeamId, corsika::qgsjetII::canInteract(corsikaBeamId));
-
-    if (corsika::qgsjetII::canInteract(corsikaBeamId)) {
-
-      CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
-
-      // position and time of interaction, not used in QgsjetII
-      Point pOrig = projectile.getPosition();
-      TimeType 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 = projectile.getNuclearA();
-
-      HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
-
-      CORSIKA_LOG_DEBUG(
-          "Interaction: ebeam lab: {} GeV"
-          "Interaction: pbeam lab: {} GeV",
-          projectileEnergyLab / 1_GeV, projectileMomentumLab.getComponents() / 1_GeV);
-      CORSIKA_LOG_DEBUG(
-          "Interaction: etarget lab: {} GeV"
-          "Interaction: ptarget lab: {} GeV",
-          targetEnergyLab / 1_GeV, targetMomentumLab.getComponents() / 1_GeV);
-      CORSIKA_LOG_DEBUG("Interaction: position of interaction: {}",
-                        pOrig.getCoordinates());
-      CORSIKA_LOG_DEBUG("Interaction: 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_);
-      CORSIKA_LOG_DEBUG("Interaction: target selected: {}", targetCode);
-
-      int targetMassNumber = 1;     // proton
-      if (is_nucleus(targetCode)) { // nucleus
-        targetMassNumber = get_nucleus_A(targetCode);
-        if (targetMassNumber > 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("Interaction: target qgsjetII code/A: {}", targetMassNumber);
-
-      int projectileMassNumber = 1; // "1" means "hadron"
-      QgsjetIIHadronType qgsjet_hadron_type =
-          qgsjetII::getQgsjetIIHadronType(corsikaBeamId);
-      if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
-        projectileMassNumber = projectile.getNuclearA();
-        if (projectileMassNumber > 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 (corsikaBeamId != Code::Nucleus) {
-        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
-      }
-
-      CORSIKA_LOG_DEBUG("Interaction: DoInteraction: E(GeV): {} GeV",
-                        projectileEnergyLab / 1_GeV);
-      count_++;
-      int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
-      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) {
-        Code idFragm = Code::Nucleus;
-        int A = fragm.getFragmentSize();
-        int Z = 0;
-        switch (A) {
-          case 1: { // proton/neutron
-            std::uniform_real_distribution<double> select;
-            if (select(rng_) > 0.5) {
-              idFragm = Code::Proton;
-              Z = 1;
-            } else {
-              idFragm = Code::Neutron;
-              Z = 0;
-            }
-
-            const HEPMassType nucleonMass = get_mass(idFragm);
-
-            auto momentum = Vector(
-                zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
-                                0.0_GeV, 0.0_GeV,
-                                sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
-                                     (projectileEnergyLabPerNucleon - nucleonMass))});
-
-            auto const energy = sqrt(momentum.getSquaredNorm() + square(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, 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 == Code::Nucleus) { // thus: not p or n
-          const HEPMassType nucleusMass = Proton::mass * Z + Neutron::mass * (A - Z);
-          auto momentum = Vector(
-              zAxisFrame, QuantityVector<hepmomentum_d>{
-                              0.0_GeV, 0.0_GeV,
-                              sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
-                                   (projectileEnergyLabPerNucleon * A - nucleusMass))});
-
-          auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleusMass));
-          momentum.rebase(originalCS); // transform back into standard lab frame
-          CORSIKA_LOG_DEBUG(
-              "secondary fragment> id={}"
-              " p={}"
-              " A={}"
-              " Z= {}",
-              idFragm, momentum.getComponents(), A, Z);
-
-          auto pnew = view.addSecondary(
-              std::make_tuple(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
-        CORSIKA_LOG_DEBUG(
-            "secondary fragment> id= {}"
-            " p= {}",
-            corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
-            momentum.getComponents());
-        auto pnew = view.addSecondary(
-            std::make_tuple(corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), energy,
-                            momentum, pOrig, tOrig));
-        Plab_final += pnew.getMomentum();
-        Elab_final += pnew.getEnergy();
-      }
-      CORSIKA_LOG_DEBUG(
-          "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
-          "Elab_final="
-          ", Plab_final={}"
-          ", 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

corsika/detail/modules/qgsjetII/InteractionModel.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/qgsjetII/InteractionModel.hpp>
+#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+#include <corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp>
+#include <corsika/modules/qgsjetII/QGSJetIIStack.hpp>
+
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.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) {
+    // initialize QgsjetII
+    corsika::connect_random_stream(rng_, ::qgsjetII::set_rng_function);
+
+    static bool initialized = false;
+    if (!initialized) {
+      CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
+      qgset_();
+      datadir DIR(dataPath.string() + "/");
+      qgaini_(DIR.data);
+      initialized = true;
+    }
+  }
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_);
+  }
+
+  inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
+                                        HEPEnergyType const sqrtS) const {
+
+    if (sqrtS < sqrtSmin_) { return false; }
+    if (is_nucleus(targetId)) {
+      size_t iTarget = get_nucleus_A(targetId);
+      if (iTarget > int(maxMassNumber_) || iTarget <= 0) { return false; }
+    } else if (targetId != Proton::code) {
+      return false;
+    }
+
+    if (is_nucleus(projectileId)) {
+      size_t iProjectile = get_nucleus_A(projectileId);
+      if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) { return false; }
+    } else if (!is_hadron(projectileId)) {
+      return false;
+    }
+    return true;
+  }
+
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+
+    if (!corsika::qgsjetII::canInteract(projectileId)) {
+      return CrossSectionType::zero();
+    }
+
+    auto const AfactorProjectile =
+        is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1;
+    auto const AfactorTarget = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
+
+    // define projectile, in lab frame
+    auto const S = (projectileP4 + targetP4).getNormSqr();
+    auto const SNN =
+        (projectileP4 / AfactorProjectile + targetP4 / AfactorTarget).getNormSqr();
+    auto const sqrtSNN = sqrt(SNN);
+    if (!isValid(projectileId, targetId, sqrtSNN)) { return CrossSectionType::zero(); }
+
+    auto const projMass = get_mass(projectileId);
+    auto const targetMass = get_mass(targetId);
+
+    // lab-frame energy per projectile nucleon as required by qgsect()
+    HEPEnergyType const ElabN =
+        calculate_lab_energy(S, projMass, targetMass) / AfactorProjectile;
+
+    int const iBeam = static_cast<QgsjetIIXSClassIntType>(
+        corsika::qgsjetII::getQgsjetIIXSCode(projectileId));
+
+    CORSIKA_LOG_DEBUG(
+        "QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
+        " iProjectile= {} iTarget= {}",
+        ElabN / 1_GeV, iBeam, AfactorProjectile, AfactorTarget);
+    double const ElabNGeV{ElabN * (1 / 1_GeV)};
+    double const sigProd = qgsect_(ElabNGeV, iBeam, AfactorProjectile, AfactorTarget);
+    CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
+    return sigProd * 1_mb;
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code projCode, Code targetCode,
+                                           FourMomentum const& proj4mom,
+                                           FourMomentum const& target4mom) const {
+    return {getCrossSection(projCode, targetCode, proj4mom, target4mom),
+            CrossSectionType::zero()};
+  }
+
+  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));
+
+    // define projectile, in lab frame
+    auto const AfactorProjectile =
+        is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1;
+    auto const AfactorTarget = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
+
+    auto const S = (projectileP4 + targetP4).getNormSqr();
+    auto const SNN =
+        (projectileP4 / AfactorProjectile + targetP4 / AfactorTarget).getNormSqr();
+    auto const sqrtSNN = sqrt(SNN);
+
+    if (!corsika::qgsjetII::canInteract(projectileId) ||
+        !isValid(projectileId, targetId, sqrtSNN)) {
+      throw std::runtime_error(fmt::format(
+          "invalid target [{}]/ projectile [{}] /energy [{} GeV] combination.",
+          get_name(targetId, full_name{}), get_name(projectileId, full_name{}),
+          sqrtSNN / 1_GeV));
+    }
+
+    auto const projMass = get_mass(projectileId);
+    auto const targetMass = get_mass(targetId);
+
+    // lab-frame energy per projectile nucleon
+    HEPEnergyType const Elab = calculate_lab_energy(S, projMass, targetMass);
+    auto const ElabN = Elab / AfactorProjectile;
+
+    CORSIKA_LOG_DEBUG("ebeam lab: {} GeV per projectile nucleon", ElabN / 1_GeV);
+
+    int const targetMassNumber = AfactorTarget;
+    CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetId, targetMassNumber);
+
+    // select QGSJetII internal projectile type
+    QgsjetIIHadronType qgsjet_hadron_type = qgsjetII::getQgsjetIIHadronType(projectileId);
+    if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
+      qgsjet_hadron_type = bernoulli_(rng_) ? QgsjetIIHadronType::ProtonType
+                                            : QgsjetIIHadronType::NeutronType;
+    } 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 const qgsjet_hadron_type_int =
+        static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
+    CORSIKA_LOG_DEBUG(
+        "qgsjet_hadron_type_int={} projectileMassNumber={} targetMassNumber={}",
+        qgsjet_hadron_type_int, AfactorProjectile, AfactorTarget);
+    qgini_(ElabN / 1_GeV, qgsjet_hadron_type_int, AfactorProjectile, AfactorTarget);
+    qgconf_();
+
+    CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
+
+    // bookkeeping
+    MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType Elab_final = 0_GeV;
+
+    // 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 / AfactorProjectile, targetP4 / AfactorTarget);
+
+    auto const& originalCS = boost.getOriginalCS();
+    auto const& csPrime =
+        boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
+
+    HEPMomentumType const pLabMag = calculate_momentum(Elab, get_mass(projectileId));
+    MomentumVector const pLab{csPrime, {0_eV, 0_eV, pLabMag}};
+
+    // internal QGSJetII system: hadron-nucleon lab. frame!
+    COMBoost const boostInternal({Elab / AfactorProjectile, pLab / AfactorProjectile},
+                                 targetMass / AfactorTarget); // felix
+
+    // fragments
+    QGSJetIIFragmentsStack qfs;
+    for (auto& fragm : qfs) {
+      int const A = fragm.getFragmentSize();
+      if (A == 1) { // nucleon
+        Code const idFragm = bernoulli_(rng_) ? Code::Proton : Code::Neutron;
+
+        HEPMassType const nucleonMass = get_mass(idFragm);
+        // no pT, fragments just go forward
+        MomentumVector const momentum{
+            csPrime, {0_eV, 0_eV, calculate_momentum(ElabN, 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{ElabN, momentum});
+        auto const P4output = boost.fromCoM(P4com);
+        auto p3output = P4output.getSpaceLikeComponents();
+        p3output.rebase(originalCS); // transform back into standard lab frame
+
+        HEPEnergyType const Ekin =
+            calculate_kinetic_energy(p3output.getNorm(), nucleonMass);
+
+        CORSIKA_LOG_DEBUG(
+            "secondary fragment> id= {}"
+            " p={}",
+            idFragm, p3output.getComponents());
+
+        auto pnew =
+            view.addSecondary(std::make_tuple(idFragm, Ekin, p3output.normalized()));
+        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);
+          }
+        }
+
+        Code const idFragm = get_nucleus_code(A, Z);
+        HEPEnergyType const mass = get_mass(idFragm);
+        // no pT, frgments just go forward
+        MomentumVector momentum{csPrime,
+                                {0.0_GeV, 0.0_GeV, calculate_momentum(ElabN * A, mass)}};
+
+        // 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{ElabN * A, momentum});
+        auto const P4output = boost.fromCoM(P4com);
+        auto p3output = P4output.getSpaceLikeComponents();
+        p3output.rebase(originalCS); // transform back into standard lab frame
+
+        HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), mass);
+
+        CORSIKA_LOG_DEBUG(
+            "secondary fragment> id={}"
+            " p={}"
+            " A={}"
+            " Z={}",
+            idFragm, p3output.getComponents(), A, Z);
+
+        auto pnew =
+            view.addSecondary(std::make_tuple(idFragm, Ekin, p3output.normalized()));
+        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
+
+      Code const pid = corsika::qgsjetII::convertFromQgsjetII(psec.getPID());
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), mass);
+
+      CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}", pid, p3output.getComponents());
+      auto pnew = view.addSecondary(std::make_tuple(pid, Ekin, p3output.normalized()));
+      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

corsika/detail/modules/qgsjetII/ParticleConversion.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/framework/core/ParticleProperties.hpp>

corsika/detail/modules/qgsjetII/QGSJetIIStack.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once

corsika/detail/modules/qgsjetII/qgsjet-II-04.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -26,14 +25,6 @@ inline datadir::datadir(const std::string& dir) {
   data[i + 1] = '\0';
 }
 
-inline double qgran_(int&) {
-  static corsika::default_prng_type& rng =
-      corsika::RNGManager::getInstance().GetRandomStream("qgran");
-
-  std::uniform_real_distribution<double> dist;
-  return dist(rng);
-}
-
 inline void lzmaopenfile_(const char*, int) {}
 inline void lzmaclosefile_() {}
 inline void lzmafillarray_(const double&, const int&) {}

corsika/detail/modules/radio/CoREAS.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/radio/CoREAS.hpp>
+
+namespace corsika {
+
+  template <typename TRadioDetector, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn CoREAS<TRadioDetector, TPropagator>::simulate(
+      Step<Particle> const& step) {
+
+    // get the global simulation time for that track.
+    auto const startTime_{
+        step.getTimePre()}; // time at the start point of the track. I
+                            // should use something similar to fCoreHitTime (?)
+    auto const endTime_{step.getTimePost()}; // time at end point of track.
+
+    if (startTime_ == endTime_) {
+      return ProcessReturn::Ok;
+    } else {
+
+      // get start and end position of the track
+      Point const startPoint_{step.getPositionPre()};
+      Point const endPoint_{step.getPositionPost()};
+      // get the coordinate system of the startpoint and hence the track
+      auto const cs_{startPoint_.getCoordinateSystem()};
+
+      auto const currDirection{(endPoint_ - startPoint_).normalized()};
+      // calculate the track length
+      auto const tracklength_{(endPoint_ - startPoint_).getNorm()};
+
+      // beta is velocity / speed of light. Start & end should be the same in endpoints!
+      auto const corrBetaValue{(endPoint_ - startPoint_).getNorm() /
+                               (constants::c * (endTime_ - startTime_))};
+      auto const beta_{currDirection * corrBetaValue};
+
+      // get particle charge
+      auto const charge_{get_charge(step.getParticlePre().getPID())};
+
+      // get thinning weight
+      auto const thinningWeight{step.getParticlePre().getWeight()};
+
+      // constants for electric field vector calculation
+      auto const constants_{charge_ * emConstant_ * thinningWeight};
+
+      // set threshold for application of ZHS-like approximation.
+      const double approxThreshold_{1.0e-3};
+
+      // loop over each observer in the observer collection (detector)
+      for (auto& observer : observers_.getObservers()) {
+
+        // get the SignalPathCollection (path1) from the start "endpoint" to the observer.
+        auto paths1{this->propagator_.propagate(step.getParticlePre(), startPoint_,
+                                                observer.getLocation())};
+
+        // get the SignalPathCollection (path2) from the end "endpoint" to the observer.
+        auto paths2{this->propagator_.propagate(step.getParticlePre(), endPoint_,
+                                                observer.getLocation())};
+
+        // LCOV_EXCL_START
+        // This should never happen unless someone implements a bad propagator
+        // good to catch this, hard to test without writing a custom and broken propagator
+        if (paths1.size() != paths2.size()) {
+          CORSIKA_LOG_CRITICAL(
+              "Signal Paths do not have the same size in CoREAS. Path starts: {} and "
+              "path ends: {}",
+              paths1.size(), paths2.size());
+        }
+        // LCOV_EXCL_STOP
+
+        // loop over both paths at once and directly compare 'start' and 'end'
+        // attributes
+        for (size_t i = 0; (i < paths1.size() && i < paths2.size()); i++) {
+
+          // calculate preDoppler factor
+          double preDoppler_{1.0 - paths1[i].refractive_index_source_ *
+                                       beta_.dot(paths1[i].emit_)};
+
+          // check if preDoppler has become zero in case of refractive index of unity
+          // because of numerical limitations here you might need std::fabs(preDoppler)
+          // in the if statement - same with post & mid
+          // LCOV_EXCL_START
+          if (preDoppler_ == 0) {
+            CORSIKA_LOG_WARN("preDoppler factor numerically zero in COREAS");
+            // redo calculation with higher precision
+            auto const& beta_components_{beta_.getComponents(cs_)};
+            auto const& emit_components_{paths1[i].emit_.getComponents(cs_)};
+            long double const indexL_{paths1[i].refractive_index_source_};
+            long double const betaX_{static_cast<double>(beta_components_.getX())};
+            long double const betaY_{static_cast<double>(beta_components_.getY())};
+            long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+            long double const startX_{static_cast<double>(emit_components_.getX())};
+            long double const startY_{static_cast<double>(emit_components_.getY())};
+            long double const startZ_{static_cast<double>(emit_components_.getZ())};
+            long double const doppler =
+                1.0l - indexL_ * (betaX_ * startX_ + betaY_ * startY_ + betaZ_ * startZ_);
+            preDoppler_ = doppler;
+          }
+          // LCOV_EXCL_STOP
+
+          // calculate postDoppler factor
+          double postDoppler_{1.0 - paths2[i].refractive_index_source_ *
+                                        beta_.dot(paths2[i].emit_)};
+
+          // check if postDoppler has become zero in case of refractive index of unity
+          // because of numerical limitations
+          // LCOV_EXCL_START
+          if (postDoppler_ == 0) {
+            CORSIKA_LOG_WARN("postDoppler factor numerically zero in CoREAS");
+            // redo calculation with higher precision
+            auto const& beta_components_{beta_.getComponents(cs_)};
+            auto const& emit_components_{paths2[i].emit_.getComponents(cs_)};
+            long double const indexL_{paths2[i].refractive_index_source_};
+            long double const betaX_{static_cast<double>(beta_components_.getX())};
+            long double const betaY_{static_cast<double>(beta_components_.getY())};
+            long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+            long double const endX_{static_cast<double>(emit_components_.getX())};
+            long double const endY_{static_cast<double>(emit_components_.getY())};
+            long double const endZ_{static_cast<double>(emit_components_.getZ())};
+            long double const doppler =
+                1.0l - indexL_ * (betaX_ * endX_ + betaY_ * endY_ + betaZ_ * endZ_);
+            postDoppler_ = doppler;
+          }
+          // LCOV_EXCL_STOP
+
+          // calculate receive time for startpoint (aka time delay)
+          auto startPointReceiveTime_{
+              startTime_ +
+              paths1[i].propagation_time_}; // TODO: time 0 is when the imaginary
+                                            // primary hits the ground
+
+          // calculate receive time for endpoint
+          auto endPointReceiveTime_{endTime_ + paths2[i].propagation_time_};
+
+          // get unit vector for startpoint at observer location
+          auto ReceiveVectorStart_{paths1[i].receive_};
+
+          // get unit vector for endpoint at observer location
+          auto ReceiveVectorEnd_{paths2[i].receive_};
+
+          // perform ZHS-like calculation close to Cherenkov angle and for refractive
+          // index at observer location greater than 1
+          if ((paths1[i].refractive_index_destination_ > 1) &&
+              ((std::fabs(preDoppler_) < approxThreshold_) ||
+               (std::fabs(postDoppler_) < approxThreshold_))) {
+            CORSIKA_LOG_DEBUG("Used ZHS-like approximation in CoREAS - radio");
+
+            // clear the existing paths for this particle and track, since we don't need
+            // them anymore
+            paths1.clear();
+            paths2.clear();
+
+            auto const halfVector_{(startPoint_ - endPoint_) * 0.5};
+            auto const midPoint_{endPoint_ + halfVector_};
+
+            // get global simulation time for the middle point of that track.
+            TimeType const midTime_{(startTime_ + endTime_) * 0.5};
+
+            // get the SignalPathCollection (path3) from the middle "endpoint" to the
+            // observer.
+            auto paths3{this->propagator_.propagate(step.getParticlePre(), midPoint_,
+                                                    observer.getLocation())};
+
+            // now loop over the paths for endpoint that we got above
+            for (auto const& path : paths3) {
+
+              auto const midPointReceiveTime_{midTime_ + path.propagation_time_};
+              double midDoppler_{1.0 -
+                                 path.refractive_index_source_ * beta_.dot(path.emit_)};
+
+              // check if midDoppler has become zero because of numerical limitations
+              // LCOV_EXCL_START
+              if (midDoppler_ == 0) {
+                CORSIKA_LOG_WARN("midDoppler factor numerically zero in COREAS");
+                // redo calculation with higher precision
+                auto const& beta_components_{beta_.getComponents(cs_)};
+                auto const& emit_components_{path.emit_.getComponents(cs_)};
+                long double const indexL_{path.refractive_index_source_};
+                long double const betaX_{static_cast<double>(beta_components_.getX())};
+                long double const betaY_{static_cast<double>(beta_components_.getY())};
+                long double const betaZ_{static_cast<double>(beta_components_.getZ())};
+                long double const midX_{static_cast<double>(emit_components_.getX())};
+                long double const midY_{static_cast<double>(emit_components_.getY())};
+                long double const midZ_{static_cast<double>(emit_components_.getZ())};
+                long double const doppler =
+                    1.0l - indexL_ * (betaX_ * midX_ + betaY_ * midY_ + betaZ_ * midZ_);
+                midDoppler_ = doppler;
+              }
+              // LCOV_EXCL_STOP
+
+              // change the values of the receive unit vectors of start and end
+              ReceiveVectorStart_ = path.receive_;
+              ReceiveVectorEnd_ = path.receive_;
+
+              // CoREAS calculation -> get ElectricFieldVector for "midPoint"
+              ElectricFieldVector EVmid_ = (path.emit_.cross(path.emit_.cross(beta_))) /
+                                           midDoppler_ / path.R_distance_ * constants_ *
+                                           observer.getSampleRate();
+
+              ElectricFieldVector EV1_{EVmid_};
+              ElectricFieldVector EV2_{EVmid_ * (-1.0)};
+
+              TimeType deltaT_{tracklength_ / (constants::c * corrBetaValue) *
+                               std::fabs(midDoppler_)}; // TODO: Caution with this!
+
+              if (startPointReceiveTime_ <
+                  endPointReceiveTime_) // EVstart_ arrives earlier
+              {
+                startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+                endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+              } else // EVend_ arrives earlier
+              {
+                startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
+                endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
+              }
+
+              TimeType const gridResolution_{1 / observer.getSampleRate()};
+              deltaT_ = endPointReceiveTime_ - startPointReceiveTime_;
+
+              // redistribute contributions over time scale defined by the observation
+              // time resolution
+              if (abs(deltaT_) < (gridResolution_)) {
+
+                EV1_ *= std::fabs((deltaT_ / gridResolution_));
+                EV2_ *= std::fabs((deltaT_ / gridResolution_));
+
+                // ToDO: be careful with times in C8!!! where is the zero (time). Is it
+                // close-by?
+                long const startBin = static_cast<long>(
+                    std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+                long const endBin = static_cast<long>(
+                    std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+                double const startBinFraction =
+                    (startPointReceiveTime_ / gridResolution_) -
+                    std::floor(startPointReceiveTime_ / gridResolution_);
+                double const endBinFraction =
+                    (endPointReceiveTime_ / gridResolution_) -
+                    std::floor(endPointReceiveTime_ / gridResolution_);
+
+                // only do timing modification if contributions would land in same bin
+                if (startBin == endBin) {
+
+                  // if startE arrives before endE
+                  if ((deltaT_) >= 0_s) {
+                    if ((startBinFraction >= 0.5) &&
+                        (endBinFraction >= 0.5)) // both points left of bin center
+                    {
+                      startPointReceiveTime_ -=
+                          gridResolution_; // shift EV1_ to previous gridpoint
+                    } else if ((startBinFraction < 0.5) &&
+                               (endBinFraction < 0.5)) // both points right of bin center
+                    {
+                      endPointReceiveTime_ +=
+                          gridResolution_; // shift EV2_ to next gridpoint
+                    } else                 // points on both sides of bin center
+                    {
+                      double const leftDist = 1.0 - startBinFraction;
+                      double const rightDist = endBinFraction;
+                      // check if asymmetry to right or left
+                      if (rightDist >= leftDist) {
+                        endPointReceiveTime_ +=
+                            gridResolution_; // shift EV2_ to next gridpoint
+                      } else {
+                        startPointReceiveTime_ -=
+                            gridResolution_; // shift EV1_ to previous gridpoint
+                      }
+                    }
+                  } else // if endE arrives before startE
+                  {
+                    if ((startBinFraction >= 0.5) &&
+                        (endBinFraction >= 0.5)) // both points left of bin center
+                    {
+                      endPointReceiveTime_ -=
+                          gridResolution_; // shift EV2_ to previous gridpoint
+                    } else if ((startBinFraction < 0.5) &&
+                               (endBinFraction < 0.5)) // both points right of bin center
+                    {
+                      startPointReceiveTime_ +=
+                          gridResolution_; // shift EV1_ to next gridpoint
+                    } else                 // points on both sides of bin center
+                    {
+                      double const leftDist = 1.0 - endBinFraction;
+                      double const rightDist = startBinFraction;
+                      // check if asymmetry to right or left
+                      if (rightDist >= leftDist) {
+                        startPointReceiveTime_ +=
+                            gridResolution_; // shift EV1_ to next gridpoint
+                      } else {
+                        endPointReceiveTime_ -=
+                            gridResolution_; // shift EV2_ to previous gridpoint
+                      }
+                    }
+                  } // End of else statement
+                }   // End of if for startbin == endbin
+              }     // End of if deltaT < gridresolution
+
+              // TODO: Be very careful with this. Maybe the EVs should be fed after the
+              // for loop of paths3
+              observer.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+              observer.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
+            } // End of looping over paths3
+
+          } // end of ZHS-like approximation
+          else {
+
+            // calculate electric field vector for startpoint
+            ElectricFieldVector EV1_ =
+                (paths1[i].emit_.cross(paths1[i].emit_.cross(beta_))) / preDoppler_ /
+                paths1[i].R_distance_ * constants_ * observer.getSampleRate();
+
+            // calculate electric field vector for endpoint
+            ElectricFieldVector EV2_ =
+                (paths2[i].emit_.cross(paths2[i].emit_.cross(beta_))) / postDoppler_ /
+                paths2[i].R_distance_ * constants_ * (-1.0) * observer.getSampleRate();
+
+            if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) {
+
+              TimeType const gridResolution_{1 / observer.getSampleRate()};
+              TimeType deltaT_{endPointReceiveTime_ - startPointReceiveTime_};
+
+              if (abs(deltaT_) < (gridResolution_)) {
+
+                EV1_ *= std::fabs(deltaT_ / gridResolution_); // Todo: rename EV1 and 2
+                EV2_ *= std::fabs(deltaT_ / gridResolution_);
+
+                long const startBin = static_cast<long>(
+                    std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+                long const endBin = static_cast<long>(
+                    std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+                double const startBinFraction =
+                    (startPointReceiveTime_ / gridResolution_) -
+                    std::floor(startPointReceiveTime_ / gridResolution_);
+                double const endBinFraction =
+                    (endPointReceiveTime_ / gridResolution_) -
+                    std::floor(endPointReceiveTime_ / gridResolution_);
+
+                // only do timing modification if contributions would land in same bin
+                if (startBin == endBin) {
+
+                  if ((startBinFraction >= 0.5) &&
+                      (endBinFraction >= 0.5)) // both points left of bin center
+                  {
+                    startPointReceiveTime_ -=
+                        gridResolution_; // shift EV1_ to previous gridpoint
+                  } else if ((startBinFraction < 0.5) &&
+                             (endBinFraction < 0.5)) // both points right of bin center
+                  {
+                    endPointReceiveTime_ +=
+                        gridResolution_; // shift EV2_ to next gridpoint
+                  } else                 // points on both sides of bin center
+                  {
+                    double const leftDist = 1.0 - startBinFraction;
+                    double const rightDist = endBinFraction;
+                    // check if asymmetry to right or left
+                    if (rightDist >= leftDist) {
+                      endPointReceiveTime_ +=
+                          gridResolution_; // shift EV2_ to next gridpoint
+                    } else {
+                      startPointReceiveTime_ -=
+                          gridResolution_; // shift EV1_ to previous gridpoint
+                    }
+                  }
+
+                } // End of if for startbin == endbin
+              }   // End of if deltaT < gridresolution
+            }     // End of if that checks small doppler factors
+            observer.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_);
+            observer.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_);
+          } // End of else that does not perform ZHS-like approximation
+
+        } // End of loop over both paths to get signal info
+      }   // End of looping over observer
+
+      return ProcessReturn::Ok;
+    }
+  } // End of simulate method
+
+} // namespace corsika

corsika/detail/modules/radio/LimitedRadioProcess.inl

+/*
+ * (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+namespace corsika {
+
+  template <class TUnderlyingProcess>
+  inline LimitedRadioProcess<TUnderlyingProcess>::LimitedRadioProcess(
+      TUnderlyingProcess&& process,
+      std::function<LimitedRadioProcess::functor_signature> runThis)
+      : process_{std::forward<TUnderlyingProcess>(process)}
+      , runThis_{std::move(runThis)} {}
+
+  template <class TUnderlyingProcess>
+  template <typename Particle>
+  inline ProcessReturn LimitedRadioProcess<TUnderlyingProcess>::doContinuous(
+      const Step<Particle>& step, const bool arg) {
+
+    if ((step.getParticlePre().getPID() != Code::Electron) &&
+        (step.getParticlePre().getPID() != Code::Positron)) {
+      CORSIKA_LOG_DEBUG("Not e+/-, will not run RadioProcess");
+      return ProcessReturn::Ok;
+    }
+
+    // Check for running the wrapped RadioProcess function
+    if (runThis_(step.getPositionPre())) {
+      CORSIKA_LOG_TRACE("Will run RadioProcess for particle at {}",
+                        step.getPositionPre());
+      return process_.doContinuous(step, arg);
+    }
+    CORSIKA_LOG_DEBUG("Particle at {} is not selected by LimitedRadioProcess",
+                      step.getPositionPre());
+    return ProcessReturn::Ok;
+  }
+
+  template <class TUnderlyingProcess>
+  template <typename Particle, typename Track>
+  inline LengthType LimitedRadioProcess<TUnderlyingProcess>::getMaxStepLength(
+      Particle const& vParticle, Track const& vTrack) const {
+    return process_.getMaxStepLength(vParticle, vTrack);
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/modules/radio/RadioProcess.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/RadioProcess.hpp>
+
+namespace corsika {
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline TRadioImpl&
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::implementation() {
+    return static_cast<TRadioImpl&>(*this);
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline TRadioImpl const&
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::implementation() const {
+    return static_cast<TRadioImpl const&>(*this);
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::RadioProcess(
+      TObserverCollection& observers, TPropagator& propagator)
+      : observers_(observers)
+      , propagator_(propagator) {}
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn RadioProcess<TObserverCollection, TRadioImpl,
+                                    TPropagator>::doContinuous(const Step<Particle>& step,
+                                                               const bool) {
+
+    // return immediately if radio process does not have any observers
+    if (observers_.size() == 0) return ProcessReturn::Ok;
+
+    // we want the following particles:
+    // Code::Electron & Code::Positron
+
+    // we wrap Simulate() in doContinuous as the plan is to add particle level
+    // filtering or thinning for calculation of the radio emission. This is
+    // important for controlling the runtime of radio (by ignoring particles
+    // that aren't going to contribute i.e. heavy hadrons)
+    // if (valid(step)) {
+    auto const particleID_{step.getParticlePre().getPID()};
+    if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
+      CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
+      return this->implementation().simulate(step);
+    } else {
+      CORSIKA_LOG_DEBUG("Particle {} is irrelevant for radio", particleID_);
+      return ProcessReturn::Ok;
+    }
+    //}
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  template <typename Particle, typename Track>
+  inline LengthType
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::getMaxStepLength(
+      [[maybe_unused]] const Particle& vParticle,
+      [[maybe_unused]] const Track& vTrack) const {
+    return meter * std::numeric_limits<double>::infinity();
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline void RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::startOfLibrary(
+      const boost::filesystem::path& directory) {
+
+    // setup the streamer
+    output_.initStreamer((directory / ("observers.parquet")).string());
+
+    // enable compression with the default level
+    output_.enableCompression();
+
+    // LCOV_EXCL_START
+    // build the schema
+    output_.addField("Time", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ex", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ey", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ez", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+    // LCOV_EXCL_STOP
+    // and build the streamer
+    output_.buildStreamer();
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline void RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::endOfShower(
+      const unsigned int) {
+
+    // loop over every observer and instruct them to
+    // flush data to disk, and then reset the observer
+    // before the next event
+    for (auto& observer : observers_.getObservers()) {
+
+      auto const sampleRate = observer.getSampleRate() * 1_s;
+      auto const radioImplementation =
+          static_cast<std::string>(this->implementation().algorithm);
+
+      // get the axis labels for this observer and write the first row.
+      axistype axis = observer.implementation().getAxis();
+
+      // get the copy of the waveform data for this event
+      std::vector<double> const& dataX = observer.implementation().getWaveformX();
+      std::vector<double> const& dataY = observer.implementation().getWaveformY();
+      std::vector<double> const& dataZ = observer.implementation().getWaveformZ();
+
+      // check for the axis name
+      std::string label = "Unknown";
+      if (observer.getDomainLabel() == "Time") {
+        label = "Time";
+      }
+      // LCOV_EXCL_START
+      else if (observer.getDomainLabel() == "Frequency") {
+        label = "Frequency";
+      }
+      // LCOV_EXCL_STOP
+      if (radioImplementation == "ZHS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          auto time = (axis.at(i + 1) + axis.at(i)) / 2.;
+          auto Ex = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
+          auto Ey = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
+          auto Ez = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
+
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(time) << static_cast<double>(Ex)
+              << static_cast<double>(Ey) << static_cast<double>(Ez) << parquet::EndRow;
+        }
+      } else if (radioImplementation == "CoREAS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(axis[i])
+              << static_cast<double>(dataX[i]) << static_cast<double>(dataY[i])
+              << static_cast<double>(dataZ[i]) << parquet::EndRow;
+        }
+      }
+
+      observer.reset();
+    }
+    output_.closeStreamer();
+
+    // increment our event counter
+    showerId_++;
+  }
+
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  inline YAML::Node
+  RadioProcess<TObserverCollection, TRadioImpl, TPropagator>::getConfig() const {
+
+    // top-level YAML node
+    YAML::Node config;
+
+    // fill in some basics
+    config["type"] = "RadioProcess";
+    config["algorithm"] = this->implementation().algorithm;
+    config["units"]["time"] = "ns";
+    config["units"]["frequency"] = "GHz";
+    config["units"]["electric field"] = "V/m";
+    config["units"]["distance"] = "m";
+
+    for (auto& observer : observers_.getObservers()) {
+      // get the name/location of this observer
+      auto name = observer.getName();
+      auto location = observer.getLocation().getCoordinates();
+
+      // get the observers config
+      config["observers"][name] = observer.getConfig();
+
+      // write the location of this observer
+      config["observers"][name]["location"].push_back(location.getX() / 1_m);
+      config["observers"][name]["location"].push_back(location.getY() / 1_m);
+      config["observers"][name]["location"].push_back(location.getZ() / 1_m);
+    }
+
+    return config;
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/ZHS.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/radio/ZHS.hpp>
+
+namespace corsika {
+
+  template <typename TRadioDetector, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn ZHS<TRadioDetector, TPropagator>::simulate(
+      Step<Particle> const& step) {
+    auto const startTime{step.getTimePre()};
+    auto const endTime{step.getTimePost()};
+
+    // LCOV_EXCL_START
+    if (startTime == endTime) {
+      return ProcessReturn::Ok;
+      // LCOV_EXCL_STOP
+    } else {
+
+      auto const startPoint{step.getPositionPre()};
+      auto const endPoint{step.getPositionPost()};
+      LengthType const trackLength{(startPoint - endPoint).getNorm()};
+
+      auto const betaModule{(endPoint - startPoint).getNorm() /
+                            (constants::c * (endTime - startTime))};
+      auto const beta{(endPoint - startPoint).normalized() * betaModule};
+
+      auto const charge{get_charge(step.getParticlePre().getPID())};
+
+      // // get "mid" position of the track geometrically
+
+      auto const halfVector{(startPoint - endPoint) / 2};
+      auto const midPoint{endPoint + halfVector};
+
+      // get thinning weight
+      auto const thinningWeight{step.getParticlePre().getWeight()};
+
+      auto const constants{charge * emConstant_ * thinningWeight};
+
+      // we loop over each observer in the collection
+      for (auto& observer : observers_.getObservers()) {
+        auto midPaths{this->propagator_.propagate(step.getParticlePre(), midPoint,
+                                                  observer.getLocation())};
+        // Loop over midPaths, first check Fraunhoffer limit
+        for (size_t i{0}; i < midPaths.size(); i++) {
+          double const uTimesK{beta.dot(midPaths[i].emit_) / betaModule};
+          double const sinTheta2{1. - uTimesK * uTimesK};
+          LengthType const lambda{constants::c / observer.getSampleRate()};
+          double const fraunhLimit{sinTheta2 * trackLength * trackLength /
+                                   midPaths[i].R_distance_ / lambda * 2 * M_PI};
+          // Checks if we are in fraunhoffer domain
+          if (fraunhLimit > 1.0) {
+            /// code for dividing track and calculating field.
+            double const nSubTracks{sqrt(fraunhLimit) + 1};
+            auto const step_{(endPoint - startPoint) / nSubTracks};
+            TimeType const timeStep{(endTime - startTime) / nSubTracks};
+            // energy should be divided up when it is possible to get the energy at end of
+            // track!!!!
+            auto point1{startPoint};
+            TimeType time1{startTime};
+            for (int j{0}; j < nSubTracks; j++) {
+              auto const point2{point1 + step_};
+              TimeType const time2{time1 + timeStep};
+              auto const newHalfVector{(point1 - point2) / 2.};
+              auto const newMidPoint{point2 + newHalfVector};
+              auto const newMidPaths{this->propagator_.propagate(
+                  step.getParticlePre(), newMidPoint, observer.getLocation())};
+              // A function for calculating the field should be made since it is repeated
+              // later
+              for (size_t k{0}; k < newMidPaths.size(); k++) {
+
+                double const n_source{newMidPaths[k].refractive_index_source_};
+
+                double const betaTimesK{beta.dot(newMidPaths[k].emit_)};
+                TimeType const midTime{(time1 + time2) / 2.};
+                TimeType detectionTime1{time1 + newMidPaths[k].propagation_time_ -
+                                        n_source * betaTimesK * (time1 - midTime)};
+                TimeType detectionTime2{time2 + newMidPaths[k].propagation_time_ -
+                                        n_source * betaTimesK * (time2 - midTime)};
+
+                // make detectionTime1_ be the smallest time =>
+                // changes step function order so sign is changed to account for it
+                double sign = 1.;
+                if (detectionTime1 > detectionTime2) {
+                  detectionTime1 = time2 + newMidPaths[k].propagation_time_ -
+                                   n_source * betaTimesK * (time2 - midTime);
+                  detectionTime2 = time1 + newMidPaths[k].propagation_time_ -
+                                   n_source * betaTimesK * (time1 - midTime);
+                  sign = -1.;
+                } // end if statement for time structure
+
+                double const startBin{
+                    std::floor((detectionTime1 - observer.getStartTime()) *
+                                   observer.getSampleRate() +
+                               0.5)};
+                double const endBin{
+                    std::floor((detectionTime2 - observer.getStartTime()) *
+                                   observer.getSampleRate() +
+                               0.5)};
+
+                auto const betaPerp{
+                    newMidPaths[k].emit_.cross(beta.cross(newMidPaths[k].emit_))};
+                double const denominator{1. - n_source * betaTimesK};
+
+                if (startBin == endBin) {
+                  // track contained in bin
+                  // if not in Cerenkov angle then
+                  if (std::fabs(denominator) > 1.e-15) {
+                    double const f{
+                        std::fabs((detectionTime2 * observer.getSampleRate() -
+                                   detectionTime1 * observer.getSampleRate()))};
+                    VectorPotential const Vp = betaPerp * sign * constants * f /
+                                               denominator / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime2, betaPerp, Vp);
+                  } else { // If emission in Cerenkov angle => approximation
+                    double const f{time2 * observer.getSampleRate() -
+                                   time1 * observer.getSampleRate()};
+                    VectorPotential const Vp =
+                        betaPerp * sign * constants * f / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime2, betaPerp, Vp);
+                  } // end if Cerenkov angle approx
+                } else {
+                  /*Track is contained in more than one bin*/
+                  int const numberOfBins{static_cast<int>(endBin - startBin)};
+                  // first contribution/ plus 1 bin minus 0.5 from new observer ruonding
+                  double f{std::fabs(startBin + 0.5 -
+                                     (detectionTime1 - observer.getStartTime()) *
+                                         observer.getSampleRate())};
+                  VectorPotential Vp = betaPerp * sign * constants * f / denominator /
+                                       newMidPaths[k].R_distance_;
+                  observer.receive(detectionTime1, betaPerp, Vp);
+                  // intermidiate contributions
+                  for (int it{1}; it < numberOfBins; ++it) {
+                    Vp = betaPerp * constants / denominator / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime1 + static_cast<double>(it) /
+                                                          observer.getSampleRate(),
+                                     betaPerp, Vp);
+                  } // end loop over bins in which potential vector is not zero
+                  // final contribution// f +0.5 from new observer rounding
+                  f = std::fabs((detectionTime2 - observer.getStartTime()) *
+                                    observer.getSampleRate() +
+                                0.5 - endBin);
+                  Vp = betaPerp * sign * constants * f / denominator /
+                       newMidPaths[k].R_distance_;
+                  observer.receive(detectionTime2, betaPerp, Vp);
+                } // end if statement for track in multiple bins
+
+              } // end of loop over newMidPaths
+              // update points for next sub track
+              point1 = point1 + step_;
+              time1 = time1 + timeStep;
+            }
+
+          } else // Calculate vector potential of whole track
+          {
+            double const n_source{midPaths[i].refractive_index_source_};
+
+            double const betaTimesK{beta.dot(midPaths[i].emit_)};
+            TimeType const midTime{(startTime + endTime) / 2};
+            TimeType detectionTime1{startTime + midPaths[i].propagation_time_ -
+                                    n_source * betaTimesK * (startTime - midTime)};
+            TimeType detectionTime2{endTime + midPaths[i].propagation_time_ -
+                                    n_source * betaTimesK * (endTime - midTime)};
+
+            // make detectionTime1_ be the smallest time =>
+            // changes step function order so sign is changed to account for it
+            double sign = 1.;
+            if (detectionTime1 > detectionTime2) {
+              detectionTime1 = endTime + midPaths[i].propagation_time_ -
+                               n_source * betaTimesK * (endTime - midTime);
+              detectionTime2 = startTime + midPaths[i].propagation_time_ -
+                               n_source * betaTimesK * (startTime - midTime);
+              sign = -1.;
+            } // end if statement for time structure
+
+            double const startBin{std::floor((detectionTime1 - observer.getStartTime()) *
+                                                 observer.getSampleRate() +
+                                             0.5)};
+            double const endBin{std::floor((detectionTime2 - observer.getStartTime()) *
+                                               observer.getSampleRate() +
+                                           0.5)};
+
+            auto const betaPerp{midPaths[i].emit_.cross(beta.cross(midPaths[i].emit_))};
+            double const denominator{1. -
+                                     midPaths[i].refractive_index_source_ * betaTimesK};
+
+            if (startBin == endBin) {
+              // track contained in bin
+              // if not in Cerenkov angle then
+              if (std::fabs(denominator) > 1.e-15) {
+                double const f{std::fabs((detectionTime2 * observer.getSampleRate() -
+                                          detectionTime1 * observer.getSampleRate()))};
+
+                VectorPotential const Vp = betaPerp * sign * constants * f / denominator /
+                                           midPaths[i].R_distance_;
+                observer.receive(detectionTime2, betaPerp, Vp);
+              } else { // If emission in Cerenkov angle => approximation
+                double const f{endTime * observer.getSampleRate() -
+                               startTime * observer.getSampleRate()};
+                VectorPotential const Vp =
+                    betaPerp * sign * constants * f / midPaths[i].R_distance_;
+                observer.receive(detectionTime2, betaPerp, Vp);
+              } // end if Cerenkov angle approx
+            } else {
+              /*Track is contained in more than one bin*/
+              int const numberOfBins{static_cast<int>(endBin - startBin)};
+              // TODO: should we check for Cerenkov angle?
+              // first contribution
+              double f{std::fabs(startBin + 0.5 -
+                                 (detectionTime1 - observer.getStartTime()) *
+                                     observer.getSampleRate())};
+              VectorPotential Vp =
+                  betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+              observer.receive(detectionTime1, betaPerp, Vp);
+              // intermediate contributions
+              for (int it{1}; it < numberOfBins; ++it) {
+                Vp = betaPerp * sign * constants / denominator / midPaths[i].R_distance_;
+                observer.receive(
+                    detectionTime1 + static_cast<double>(it) / observer.getSampleRate(),
+                    betaPerp, Vp);
+              } // end loop over bins in which potential vector is not zero
+              // final contribution
+              f = std::fabs((detectionTime2 - observer.getStartTime()) *
+                                observer.getSampleRate() +
+                            0.5 - endBin);
+              Vp =
+                  betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+              observer.receive(detectionTime2, betaPerp, Vp);
+            } // end if statement for track in multiple bins
+
+          } // finish if statement of track in fraunhoffer or not
+
+        } // end loop over mid paths
+
+      } // END: loop over observers
+      return ProcessReturn::Ok;
+    }
+  } // end simulate
+
+} // namespace corsika

corsika/detail/modules/radio/detectors/ObserverCollection.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <corsika/modules/radio/detectors/ObserverCollection.hpp>
+
+namespace corsika {
+
+  template <typename TObserverImpl>
+  inline void ObserverCollection<TObserverImpl>::addObserver(
+      TObserverImpl const& observer) {
+    observers_.push_back(observer);
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl& ObserverCollection<TObserverImpl>::at(std::size_t const i) {
+    return observers_.at(i);
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl const& ObserverCollection<TObserverImpl>::at(
+      std::size_t const i) const {
+    return observers_.at(i);
+  }
+
+  template <typename TObserverImpl>
+  inline int ObserverCollection<TObserverImpl>::size() const {
+    return observers_.size();
+  }
+
+  template <typename TObserverImpl>
+  inline std::vector<TObserverImpl>& ObserverCollection<TObserverImpl>::getObservers() {
+    return observers_;
+  }
+
+  template <typename TObserverImpl>
+  inline void ObserverCollection<TObserverImpl>::reset() {
+    std::for_each(observers_.begin(), observers_.end(),
+                  std::mem_fn(&TObserverImpl::reset));
+  }
+
+} // namespace corsika