diff --git a/corsika/detail/modules/pythia8/NeutrinoInteraction.inl b/corsika/detail/modules/pythia8/NeutrinoInteraction.inl
index 9c256e72abdc889fb9cd13d32f14ae732137b0e4..183d360c27848da8d5bae7902041dc8742db048d 100644
--- a/corsika/detail/modules/pythia8/NeutrinoInteraction.inl
+++ b/corsika/detail/modules/pythia8/NeutrinoInteraction.inl
@@ -27,69 +27,72 @@ namespace corsika::pythia8 {
// Main Pythia object for managing the cascade evolution.
// Can also do decays, but no hard processes.
- pythiaMain_.readString("ProcessLevel:all = off");
+ // pythiaMain_.readString("ProcessLevel:all = off");
- // Reduce statistics printout to relevant ones.
+ // // 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");
+ // // 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");
+ // }
+ // if (!handle_nc_ && !handle_cc_) {
+ // CORSIKA_LOG_ERROR("wrong configuration. Need either CC or NC interactions!");
+ // }
+ // pythiaMain_.settings.parm("PhaseSpace:Q2Min", minQ2_ / 1_GeV);
+
+ // // 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
+ // // 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() {
@@ -102,11 +105,61 @@ namespace corsika::pythia8 {
FourMomentum const& projectileP4,
FourMomentum const& targetP4) {
- CORSIKA_LOG_DEBUG("Pythia::NeutrinoInteraction: doInteraction: {} ", projectileId);
+ CORSIKA_LOG_INFO("Pythia::NeutrinoInteraction: doInteraction: {} E={}GeV",
+ projectileId, projectileP4.getTimeLikeComponent() / 1_GeV);
+ // allow just one call to NeutrinoInteraction
if (!count_) {
- // References to the two event records. Clear main event record.
+ double eProton = Proton::mass / 1_GeV; // 0.93827;
+ double eElectron = projectileP4.getTimeLikeComponent() / 1_GeV; // 270.5;
+ double Q2min = minQ2_ / 1_GeV;
+
+ // Set up incoming beams, for frame with unequal beam energies.
+ pythiaMain_.readString("Beams:frameType = 2");
+ // BeamA = proton.
+ pythiaMain_.readString("Beams:idB = 2212");
+ pythiaMain_.settings.parm("Beams:eB", eProton);
+ // BeamB = electron.
+ int const idNow = static_cast<int>(get_PDG(projectileId));
+ pythiaMain_.readString("Beams:idA = 12");
+
+ // pythiaMain_.settings.parm("Beams:idA", idNow);
+ pythiaMain_.settings.parm("Beams:eA", eElectron);
+
+ // Set up DIS process within some phase space.
+ // Neutral current (with gamma/Z interference).
+ pythiaMain_.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
+ // Uncomment to allow charged current.
+ pythiaMain_.readString("WeakBosonExchange:ff2ff(t:W) = on");
+ // Phase-space cut: minimal Q2 of process.
+ pythiaMain_.settings.parm("PhaseSpace:Q2Min", Q2min);
+
+ // 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");
+
+ pythiaMain_.readString("Stat:showProcessLevel = off");
+ pythiaMain_.readString("Stat:showPartonLevel = off");
+
+ pythiaMain_.readString("Print:quiet = on");
+ pythiaMain_.readString("Check:epTolErr = 0.1");
+ pythiaMain_.readString("Check:epTolWarn = 0.0001");
+ pythiaMain_.readString("Check:mTolErr = 0.01");
+
+ pythiaMain_.init();
+
+ // References to the event record
Pythia8::Event& eventMain = pythiaMain_.event;
COMBoost const labFrameBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
@@ -155,7 +208,8 @@ namespace corsika::pythia8 {
// // 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
+ // // skip particles that have decayed / are initial particles in pythia's
+ // event
// // record
// if (!p8p.isFinal()) continue;
// try {
diff --git a/corsika/modules/pythia8/NeutrinoInteraction.hpp b/corsika/modules/pythia8/NeutrinoInteraction.hpp
index 8bb449b48fb943346bd32a013b92dafd76d4e7c8..bbc47dac1636e3fcc921de354273f46720f32638 100644
--- a/corsika/modules/pythia8/NeutrinoInteraction.hpp
+++ b/corsika/modules/pythia8/NeutrinoInteraction.hpp
@@ -68,6 +68,7 @@ namespace corsika::pythia8 {
bool const print_listing_ = false;
bool const handle_nc_;
bool const handle_cc_;
+ HEPEnergyType const minQ2_ = 25_GeV;
Pythia8::Pythia pythiaMain_;
};
} // namespace corsika::pythia8