add neutrino interaction for pythia

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 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

corsika/modules/pythia8/NeutrinoInteraction.hpp

+/*
+ * (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