From 0a57fa5d6e553d3770cce6918dc280fdce850bf1 Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Wed, 14 Feb 2024 18:56:17 +0100
Subject: [PATCH] add neutrino interaction for pythia
---
.../modules/pythia8/NeutrinoInteraction.inl | 201 ++++++++++++++++++
.../modules/pythia8/NeutrinoInteraction.hpp | 75 +++++++
2 files changed, 276 insertions(+)
create mode 100644 corsika/detail/modules/pythia8/NeutrinoInteraction.inl
create mode 100644 corsika/modules/pythia8/NeutrinoInteraction.hpp
diff --git a/corsika/detail/modules/pythia8/NeutrinoInteraction.inl b/corsika/detail/modules/pythia8/NeutrinoInteraction.inl
new file mode 100644
index 000000000..9c256e72a
--- /dev/null
+++ b/corsika/detail/modules/pythia8/NeutrinoInteraction.inl
@@ -0,0 +1,201 @@
+/*
+ * (c) Copyright 2024 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 <boost/filesystem/path.hpp>
+
+namespace corsika::pythia8 {
+
+ inline NeutrinoInteraction::NeutrinoInteraction(HEPEnergyType const energyLab,
+ bool const& handleNC,
+ bool const& handleCC,
+ bool const print_listing)
+ : print_listing_(print_listing)
+ , handle_nc_(handleNC)
+ , handle_cc_(handleCC)
+ , pythiaMain_{CORSIKA_Pythia8_XML_DIR, false} {
+ Pythia8::RndmEngine* rndm = new corsika::pythia8::Random();
+ pythiaMain_.setRndmEnginePtr(rndm);
+
+ // CORSIKA_LOG_INFO("Pythia8 MPI init file: {}", mpiInitFile.native());
+ // Main Pythia object for managing the cascade evolution.
+ // Can also do decays, but no hard processes.
+
+ pythiaMain_.readString("ProcessLevel:all = off");
+
+ // Reduce statistics printout to relevant ones.
+ // pythiaMain_.readString("Stat:showProcessLevel = off");
+ // pythiaMain_.readString("Stat:showPartonLevel = off");
+
+ // Set up for fixed-target collisions.
+ pythiaMain_.readString(
+ "Beams:frameType = 2"); // arbitrary frame, need to define full 4-momenta
+ pythiaMain_.settings.parm("Beams:eA", energyLab / 1_GeV);
+ pythiaMain_.readString("Beams:idA = 12"); // electron neutrino
+ pythiaMain_.settings.parm("Beams:eB", Proton::mass / 1_GeV);
+ pythiaMain_.readString("Beams:idB = 2212");
+
+ // // Set up incoming beams, for frame with unequal beam energies.
+ // pythia.readString("Beams:frameType = 2");
+ // // BeamA = proton.
+ // pythia.readString("Beams:idA = 2212");
+ // pythia.settings.parm("Beams:eA", eProton);
+ // // BeamB = electron.
+ // pythia.readString("Beams:idB = 11");
+ // pythia.settings.parm("Beams:eB", eElectron);
+
+ // switch on Z and W exchange. these won't affect the hadrons much but will allow
+ // for neutrino primaries
+ if (handle_nc_) {
+ CORSIKA_LOG_INFO("initializing pythia for neutral currents..");
+ pythiaMain_.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
+ }
+ if (handle_cc_) {
+ CORSIKA_LOG_INFO("initializing pythia for charged currents..");
+ pythiaMain_.readString("WeakBosonExchange:ff2ff(t:W) = on");
+ }
+ pythiaMain_.settings.parm("PhaseSpace:Q2Min", 25.);
+
+ // 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");
+
+ // no Decays to be done by pythiaMain_.
+ pythiaMain_.readString("HadronLevel:Decay = off");
+
+ // Reduce printout and relax energy-momentum conservation.
+ // pythiaMain_.readString("Print:quiet = on");
+ pythiaMain_.readString("Check:epTolErr = 0.1");
+ pythiaMain_.readString("Check:epTolWarn = 0.0001");
+ pythiaMain_.readString("Check:mTolErr = 0.01");
+
+ // Redure statistics printout to relevant ones.
+ pythiaMain_.readString("Stat:showProcessLevel = off");
+ pythiaMain_.readString("Stat:showPartonLevel = off");
+
+ // we can't test this block, LCOV_EXCL_START
+ if (!pythiaMain_.init())
+ throw std::runtime_error("Pythia::NeutrinoInteraction: Initialization failed!");
+ // LCOV_EXCL_STOP
+ }
+
+ inline NeutrinoInteraction::~NeutrinoInteraction() {
+ CORSIKA_LOG_INFO("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_LOG_DEBUG("Pythia::NeutrinoInteraction: doInteraction: {} ", projectileId);
+
+ if (!count_) {
+
+ // References to the two event records. Clear main event record.
+ Pythia8::Event& eventMain = pythiaMain_.event;
+
+ COMBoost const labFrameBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+ auto const proj4MomLab = labFrameBoost.toCoM(projectileP4);
+ auto const& rotCS = labFrameBoost.getRotatedCS();
+
+ if (!pythiaMain_.next()) {
+ throw std::runtime_error("Pythia neutrino collision failed ");
+ } else {
+ CORSIKA_LOG_INFO("pythia neutrino interaction done!");
+ }
+
+ MomentumVector Plab_final{labFrameBoost.getOriginalCS()};
+ auto Elab_final = HEPEnergyType::zero();
+ CORSIKA_LOG_INFO("pythia stack size {}", eventMain.size());
+ for (int i = 0; i < eventMain.size(); ++i) {
+ if (eventMain[i].isFinal()) {
+ auto const& p8p = eventMain[i];
+ try {
+ auto const volatile id = static_cast<PDGCode>(p8p.id());
+ auto const pyId = convert_from_PDG(id);
+
+ MomentumVector const pyPlab(
+ rotCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+ auto const pyP = labFrameBoost.fromCoM(FourVector{p8p.e() * 1_GeV, pyPlab});
+
+ 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, pyPlab.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
+ }
+ }
+
+ // // copy particles from pythia stack to corsika
+ // for (Pythia8::Particle const& p8p : eventMain) {
+ // // skip particles that have decayed / are initial particles in pythia's event
+ // // record
+ // if (!p8p.isFinal()) continue;
+ // try {
+ // auto const volatile id = static_cast<PDGCode>(p8p.id());
+ // auto const pyId = convert_from_PDG(id);
+
+ // MomentumVector const pyPlab(
+ // rotCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+ // auto const pyP = labFrameBoost.fromCoM(FourVector{p8p.e() * 1_GeV,
+ // pyPlab});
+
+ // 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,
+ // pyPlab.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
+ // }
+
+ // eventMain.clear();
+
+ CORSIKA_LOG_DEBUG(
+ "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
diff --git a/corsika/modules/pythia8/NeutrinoInteraction.hpp b/corsika/modules/pythia8/NeutrinoInteraction.hpp
new file mode 100644
index 000000000..8bb449b48
--- /dev/null
+++ b/corsika/modules/pythia8/NeutrinoInteraction.hpp
@@ -0,0 +1,75 @@
+/*
+ * (c) Copyright 2024 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 <tuple>
+#include <boost/filesystem/path.hpp>
+
+#include <corsika/framework/utility/CorsikaData.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/modules/pythia8/Pythia8.hpp>
+
+namespace corsika::pythia8 {
+
+ class NeutrinoInteraction : public InteractionProcess<NeutrinoInteraction> {
+
+ public:
+ NeutrinoInteraction(HEPEnergyType const, bool const& handleNC = true,
+ bool const& handleCC = true, bool const print_listing = false);
+ ~NeutrinoInteraction();
+ /**
+ * Returns inelastic (production) cross section.
+ *
+ * This cross section must correspond to the process described in doInteraction.
+ * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+ *
+ * @param projectile is the Code of the projectile
+ * @param target is the Code of the target
+ * @param projectileP4 is the four momentum of the projectile
+ * @param targetP4 is the four momentum of the target
+ *
+ * @return inelastic cross section
+ */
+ CrossSectionType getCrossSection(Code const projectile, Code const target,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ if (isValid(projectile, target, projectileP4, targetP4))
+ return cross_section_;
+ else
+ return CrossSectionType::zero();
+ };
+
+ bool isValid(Code const projectileId, [[maybe_unused]] Code const targetId,
+ [[maybe_unused]] FourMomentum const& projectileP4,
+ [[maybe_unused]] FourMomentum const& targetP4) const {
+ return is_neutrino(projectileId);
+ };
+
+ /**
+ * In this function PYTHIA is called to produce one event. The
+ * event is copied (and boosted) into the shower lab frame.
+ */
+ template <typename TView>
+ void doInteraction(TView& output, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+ private:
+ CrossSectionType const cross_section_ = 4_nb;
+ int count_ = 0;
+ bool const print_listing_ = false;
+ bool const handle_nc_;
+ bool const handle_cc_;
+ Pythia8::Pythia pythiaMain_;
+ };
+} // namespace corsika::pythia8
+
+#include <corsika/detail/modules/pythia8/NeutrinoInteraction.inl>
\ No newline at end of file
--
GitLab