Neutrino interactions with pythia 8310

Never mind the branch name. This uses the latest pythia 8.310.

applications/c8_air_shower.cpp

@@ -147,7 +147,9 @@ int main(int argc, char** argv) {
                    ->needs(opt_Z)
                    ->check(CLI::Range(1, 58))
                    ->group("Primary");
-  app.add_option("-p,--pdg", "PDG code for primary.")
+  app.add_option("-p,--pdg",
+                 "PDG code for primary (p=2212, gamma=22, e-=11, nu_e=12, mu-=13, "
+                 "nu_mu=14, tau=15, nu_tau=16).")
       ->excludes(opt_A)
       ->excludes(opt_Z)
       ->group("Primary");
@@ -178,6 +180,14 @@ int main(int argc, char** argv) {
       ->default_val(0.3)
       ->check(CLI::Range(0.000001, 1.e13))
       ->group("Config");
+  bool track_neutrinos = false;
+  app.add_flag("--track-neutrinos", track_neutrinos, "switch on tracking of neutrinos")
+      ->group("Config");
+  app.add_option("--neutrino-interaction-type",
+                 "charged (CC) or neutral current (NC) or both")
+      ->default_val("both")
+      ->check(CLI::IsMember({"neutral", "NC", "charged", "CC", "both"}))
+      ->group("Misc.");
   app.add_option("--observation-level",
                  "Height above earth radius of the observation level (in m)")
       ->default_val(0.)
@@ -240,7 +250,6 @@ int main(int argc, char** argv) {
       ->default_val(0)
       ->check(CLI::Range(0, 20))
       ->group("Radio");
-
   // parse the command line options into the variables
   CLI11_PARSE(app, argc, argv);
 
@@ -397,13 +406,30 @@ int main(int argc, char** argv) {
 
   corsika::pythia8::Decay decayPythia;
 
+  // neutrino interactions with pythia (options are: NC, CC)
+  bool NC = false;
+  bool CC = false;
+  if (auto const nuIntStr = app["--neutrino-interaction-type"]->as<std::string>();
+      nuIntStr == "neutral" || nuIntStr == "NC") {
+    NC = true;
+    CC = false;
+  } else if (nuIntStr == "charged" || nuIntStr == "CC") {
+    NC = false;
+    CC = true;
+  } else if (nuIntStr == "both") {
+    NC = true;
+    CC = true;
+  }
+  corsika::pythia8::NeutrinoInteraction neutrinoPrimaryPythia(NC, CC);
+
   // hadronic photon interactions in resonance region
   corsika::sophia::InteractionModel sophia;
 
   HEPEnergyType const emcut = 1_GeV * app["--emcut"]->as<double>();
   HEPEnergyType const hadcut = 1_GeV * app["--hadcut"]->as<double>();
   HEPEnergyType const mucut = 1_GeV * app["--mucut"]->as<double>();
-  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut, true, dEdX);
+  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut,
+                                             !track_neutrinos, dEdX);
 
   // tell proposal that we are interested in all energy losses above the particle cut
   set_energy_production_threshold(Code::Electron, std::min({emcut, hadcut, mucut}));
@@ -541,9 +567,9 @@ int main(int argc, char** argv) {
   output.add("ZHS", zhs);
 
   // assemble the final process sequence with radio
-  auto sequence =
-      make_sequence(stackInspect, hadronSequence, decayPythia, emCascade, emContinuous,
-                    coreas, zhs, longprof, observationLevel, thinning, cut);
+  auto sequence = make_sequence(stackInspect, neutrinoPrimaryPythia, hadronSequence,
+                                decayPythia, emCascade, emContinuous, coreas, zhs,
+                                longprof, observationLevel, thinning, cut);
 
   /* === END: SETUP PROCESS LIST === */
 
@@ -554,6 +580,7 @@ int main(int argc, char** argv) {
   Cascade EAS(env, tracking, sequence, output, stack);
 
   // print our primary parameters all in one place
+  CORSIKA_LOG_INFO("Primary name: {}", beamCode);
   if (app["--pdg"]->count() > 0) {
     CORSIKA_LOG_INFO("Primary PDG ID:     {}", app["--pdg"]->as<int>());
   } else {
@@ -591,7+618,7 @@
     if (force_interaction) {
       CORSIKA_LOG_INFO("Fixing first interaction at injection point.");
       EAS.forceInteraction();
     }
 
     primaryWriter.recordPrimary(primaryProperties);