IAP GITLAB

Skip to content
Snippets Groups Projects

Neutrino interactions with pythia 8310

Merged Neutrino interactions with pythia 8310
neutrino-interactions-with-pythia-8245 into master
Merged Felix Riehn requested to merge neutrino-interactions-with-pythia-8245 into master
Viewing commit 0a57fa5d
Next
Show latest version
2 files
+ 276
− 0
Compare changes
  • Side-by-side
  • Inline
Files
2
corsika/detail/modules/pythia8/NeutrinoInteraction.inl 0 → 100644
+ 201
− 0
/*
* (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

data privacy - impressum