Neutrino interactions with pythia 8310

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

corsika/detail/modules/pythia8/NeutrinoInteraction.inl

@@ -20,9 +20,14 @@ namespace corsika::pythia8 {
       , 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("configure Pythia for primary neutrino interactions. NC={}, CC={}",
+                     handle_nc_, handle_cc_);
+    CORSIKA_LOG_INFO("minimal Q2 in DIS: {} GeV2", minQ2_ / 1_GeV / 1_GeV);
+
     // 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.
@@ -105,27 +110,25 @@ namespace corsika::pythia8 {
                                           FourMomentum const& projectileP4,
                                           FourMomentum const& targetP4) {
 
-    CORSIKA_LOG_INFO("Pythia::NeutrinoInteraction: doInteraction: {} E={}GeV",
-                     projectileId, projectileP4.getTimeLikeComponent() / 1_GeV);
+    CORSIKA_LOG_INFO("Primary {} - p interaction with E_nu = {} GeV", projectileId,
+                     projectileP4.getTimeLikeComponent() / 1_GeV);
 
     // allow just one call to NeutrinoInteraction
     if (!count_) {
 
-      double eProton = Proton::mass / 1_GeV;                          // 0.93827;
-      double eElectron = projectileP4.getTimeLikeComponent() / 1_GeV; // 270.5;
-      double Q2min = minQ2_ / 1_GeV;
+      double const eProton = Proton::mass / 1_GeV;                          // 0.93827;
+      double const eElectron = projectileP4.getTimeLikeComponent() / 1_GeV; // 270.5;
+      double const Q2min = minQ2_ / 1_GeV / 1_GeV;
+      int const idNow = static_cast<int>(get_PDG(projectileId));
 
       // Set up incoming beams, for frame with unequal beam energies.
       // projectile is along -z
       pythiaMain_.readString("Beams:frameType = 2");
-      // BeamA = proton.
+      // BeamA = proton.(+z)
       pythiaMain_.readString("Beams:idA = 2212");
       pythiaMain_.settings.parm("Beams:eA", eProton);
-      // BeamB = electron.
-      int const idNow = static_cast<int>(get_PDG(projectileId));
-      pythiaMain_.readString("Beams:idB = 12");
-
-      // pythiaMain_.settings.parm("Beams:idA", idNow);
+      // BeamB = neutrino (-z direction)
+      pythiaMain_.settings.mode("Beams:idB", idNow);
       pythiaMain_.settings.parm("Beams:eB", eElectron);
 
       // Set up DIS process within some phase space.
@@ -175,7 +178,7 @@ namespace corsika::pythia8 {
 
       MomentumVector Plab_final{labFrameBoost.getOriginalCS()};
       auto Elab_final = HEPEnergyType::zero();
-      CORSIKA_LOG_INFO("pythia stack size {}", eventMain.size());
+      CORSIKA_LOG_INFO("particles generated in neutrino interaction:");
       for (int i = 0; i < eventMain.size(); ++i) {
         if (eventMain[i].isFinal()) {
           auto const& p8p = eventMain[i];
@@ -195,6 +198,8 @@ namespace corsika::pythia8 {
             auto pnew =
                 view.addSecondary(std::make_tuple(pyId, Ekin, pyPlab.normalized()));
 
+            CORSIKA_LOG_INFO("id = {}, E = {} GeV, p = {} GeV", pyId, Ekin / 1_GeV,
+                             pyPlab.getComponents() / 1_GeV);
             Plab_final += pnew.getMomentum();
             Elab_final += pnew.getEnergy();
           }
@@ -249,6 +254,8 @@ namespace corsika::pythia8 {
           "Elab_final= {}"
           ", Plab_final= {}",
           Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
+    } else {
+      CORSIKA_LOG_ERROR("Only one primary neutrino interaction allowed!");
     }
     count_++;
   }