Neutrino interactions with pythia 8310

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

applications/c8_air_shower.cpp

@@ -180,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.)
@@ -398,8 +406,21 @@ int main(int argc, char** argv) {
 
   corsika::pythia8::Decay decayPythia;
 
-  // neutrino interactions (options are: NC, CC)
-  corsika::pythia8::NeutrinoInteraction neutrinoPrimaryPythia(true, true);
+  // 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;
@@ -407,7 +428,8 @@ int main(int argc, char** argv) {
   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}));