corsika/detail/framework/core/Cascade.inl

@@ -296,22 +296,6 @@ namespace corsika {
 
       CORSIKA_LOG_DEBUG("step limit reached (e.g. deflection). nothing further happens.");
 
-      // final sanity check, no actions
-      {
-        auto const* numericalNodeAfterStep =
-            environment_.getUniverse()->getContainingNode(particle.getPosition());
-        CORSIKA_LOG_TRACE(
-            "Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
-            fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
-        if (numericalNodeAfterStep != currentLogicalNode) {
-          CORSIKA_LOG_ERROR(
-              "expect to be in node currentLogicalNode={} but are in "
-              "numericalNodeAfterStep={}. Continue, but without guarantee.",
-              fmt::ptr(currentLogicalNode), fmt::ptr(numericalNodeAfterStep));
-        }
-      }
-      // we did not cross any volume boundary
-
       // step length limit
       return;
     }

corsika/detail/modules/pythia8/NeutrinoInteraction.inl

@@ -58,7 +58,10 @@ namespace corsika::pythia8 {
         corsika::RNGManager<>::getInstance().getRandomStream("pythia");
     Code const targetNucleonId = (nucleonChannelDist(rng) ? Code::Neutron : Code::Proton);
     int const idTarget_pythia = static_cast<int>(get_PDG(targetNucleonId));
-    CORSIKA_LOGGER_DEBUG(logger_, "selected {} target", 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;
@@ -115,10 +118,10 @@ namespace corsika::pythia8 {
     // References to the event record
     Pythia8::Event& eventMain = pythiaMain_.event;
     // define boost assuming target nucleon is at rest
-    COMBoost const labFrameBoost{projectileP4, get_mass(targetNucleonId)};
+    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 = labFrameBoost.getRotatedCS();
+    auto const& rotCS = boost.getRotatedCS();
 
     if (!pythiaMain_.next()) {
       throw std::runtime_error("Pythia neutrino collision failed ");
@@ -126,7 +129,7 @@ namespace corsika::pythia8 {
       CORSIKA_LOGGER_DEBUG(logger_, "pythia neutrino interaction done!");
     }
 
-    MomentumVector Plab_final{labFrameBoost.getOriginalCS()};
+    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) {
@@ -143,7 +146,7 @@ namespace corsika::pythia8 {
           // calculate 3-momentum and kinetic energy
           CORSIKA_LOGGER_DEBUG(logger_, "momentum in CoM: {} GeV",
                                pyPcom.getComponents() / 1_GeV);
-          auto const pyP4lab = labFrameBoost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPcom});
+          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));

corsika/modules/pythia8/NeutrinoInteraction.hpp

@@ -22,9 +22,21 @@ namespace corsika::pythia8 {
 
   using HEPEnergyTypeSqr = decltype(1_GeV * 1_GeV);
 
+  /**
+   * @brief Defines the interface to PYTHIA8. Configured for
+   * DIS neutrino - nucleon interactions.
+   *
+   */
+
   class NeutrinoInteraction : public InteractionProcess<NeutrinoInteraction> {
 
   public:
+    /**
+     * Constructs the interface for PYTHIA8 neutrino interactions
+     *
+     * @param handleNC - Switch on/off neutral current interactions
+     * @param handleCC - Switch on/off charged current interactions
+     */
     NeutrinoInteraction(bool const& handleNC = true, bool const& handleCC = true);
     ~NeutrinoInteraction();
     /**

python/examples/shower_profile.py

@@ -72,7 +72,7 @@ def draw_profiles(pr, pr_config):
     ax.set_title(pr_config["type"] + "\n" + title)
     unit_grammage_str = pr_config["units"]["grammage"]  # get units of simulation output
     ax.set_xlabel(f"slant depth, X ({unit_grammage_str})")
-    ax.set_ylabel("dN/dX")
+    ax.set_ylabel("N(X)")
     ax.legend()
     ax.set_yscale("log")
     ax.set_xlim(min(pr["X"]), max(pr["X"]))