From f5a0728b7836934654414740ba0a810b80d372bb Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Tue, 24 Sep 2024 00:22:19 +0000
Subject: [PATCH] Resolve "fix boost in pythia8 neutrino interface"
---
.../detail/modules/pythia8/NeutrinoInteraction.inl | 13 ++++++++-----
corsika/modules/pythia8/NeutrinoInteraction.hpp | 12 ++++++++++++
2 files changed, 20 insertions(+), 5 deletions(-)
diff --git a/corsika/detail/modules/pythia8/NeutrinoInteraction.inl b/corsika/detail/modules/pythia8/NeutrinoInteraction.inl
index 67f98ae79..3830e3917 100644
--- a/corsika/detail/modules/pythia8/NeutrinoInteraction.inl
+++ b/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));
diff --git a/corsika/modules/pythia8/NeutrinoInteraction.hpp b/corsika/modules/pythia8/NeutrinoInteraction.hpp
index 804c907ec..2a8946eb6 100644
--- a/corsika/modules/pythia8/NeutrinoInteraction.hpp
+++ b/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();
/**
--
GitLab