CMakeLists.txt

@@ -87,7 +87,7 @@ include (corsikaDefines)
 # check if compiler is C++17 compliant
 #
 include (CheckCXXCompilerFlag)
-check_CXX_compiler_flag ("--std=c++17" COMPILER_SUPPORTS_CXX17)
+check_CXX_compiler_flag ("-std=c++17" COMPILER_SUPPORTS_CXX17)
 if (NOT COMPILER_SUPPORTS_CXX17)
  message (FATAL "| CORSIKA8 > The compiler ${CMAKE_CXX_COMPILER} has no C++17 support. Please use a different C++ compiler.")
 endif ()

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"]))