added primitive application to calculate interaction cross sections in hadronic interaction models

applications/CMakeLists.txt

@@ -13,3 +13,10 @@ endif()
 install (
   TARGETS c8_air_shower DESTINATION bin
 )
+
+add_executable (cross_section cross_section.cpp)
+target_link_libraries (cross_section CORSIKA8)
+
+install (
+  TARGETS cross_section DESTINATION bin
+)

applications/cross_section.cpp

+/*
+ * (c) Copyright 2019 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.
+ */
+
+#include <corsika/modules/Epos.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+#include <corsika/modules/Pythia8.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+//#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/random/RNGManager.hpp>
+
+// #include <corsika/media/Environment.hpp>
+// #include <corsika/media/HomogeneousMedium.hpp>
+// #include <corsika/media/MediumPropertyModel.hpp>
+// #include <corsika/media/NuclearComposition.hpp>
+// #include <corsika/media/ShowerAxis.hpp>
+// #include <corsika/media/UniformMagneticField.hpp>
+
+// #include <corsika/setup/SetupStack.hpp>
+
+#include <fstream>
+#include <string>
+#include <tuple>
+
+/*
+  NOTE, WARNING, ATTENTION
+
+  The Random.hpp implements the hook of epos to the C8 random
+  number generator. It has to occur excatly ONCE per linked
+  executable. If you include the header below in multiple "tests" and
+  link them togehter, it will fail.
+ */
+#include <corsika/modules/Random.hpp>
+
+using namespace corsika;
+// using namespace corsika::epos;
+
+template <typename TModel>
+void calculate_cross_sections(TModel & model, std::string const &model_name,
+                              bool lab_frame = true) {
+
+  std::stringstream fname;
+  fname << "cross-sections-" << model_name << ".txt";
+  std::ofstream out(fname.str());
+  out << "# cross section table for " << model_name << "\n";
+  out << "# i, elab, sqsnn, sig_pp, sig_pO \n";
+
+  CoordinateSystemPtr const &cs = get_root_CoordinateSystem();
+
+  float const amin = 2;
+  float const amax = 12;
+  int const nebins = 15;
+  float const da = (amax - amin) / (float)nebins;
+
+  for (int i = 0; i < nebins; ++i) {
+
+    corsika::units::si::HEPMomentumType const setP =
+        pow(10, da * i + amin) * 1_GeV;
+    
+    auto setE = calculate_total_energy(setP, Proton::mass);
+    corsika::units::si::HEPEnergyType const comEnn =
+        corsika::calculate_com_energy(setE, Proton::mass,
+                                      corsika::constants::nucleonMass);
+
+    //corsika::units::si::HEPMomentumType setP =
+    //    corsika::calculate_momentum(setE, Proton::mass);
+
+    FourMomentum p4protonProj(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
+    FourMomentum p4protonTarg(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
+    FourMomentum p4oxygenTarg(0_eV, {cs, 0_GeV, 0_GeV, 0_GeV});
+
+    if (lab_frame) {
+      // four momenta in lab frame
+      p4protonProj =
+          corsika::FourMomentum(setE, MomentumVector(cs, {0_eV, 0_eV, setP}));
+      p4protonTarg = corsika::FourMomentum(
+          Proton::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+      p4oxygenTarg = corsika::FourMomentum(
+          Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+    } else {
+      // four momenta in com frame
+      p4protonProj = corsika::FourMomentum(
+          corsika::calculate_total_energy(setP / 2, Proton::mass),
+          MomentumVector(cs, {0_eV, 0_eV, setP / 2}));
+      p4protonTarg = corsika::FourMomentum(
+          corsika::calculate_total_energy(setP / 2, Proton::mass),
+          MomentumVector(cs, {0_eV, 0_eV, -setP / 2}));
+      p4oxygenTarg = corsika::FourMomentum(
+          corsika::calculate_total_energy(setP / 2, Oxygen::mass),
+          MomentumVector(cs, {0_eV, 0_eV, -setP / 2}));
+    }
+
+    // p-p
+    CrossSectionType const xs_prod_pp = model->getCrossSection(
+        Code::Proton, Code::Proton, p4protonProj, p4protonTarg);
+    // p-O
+    CrossSectionType const xs_prod_pO = model->getCrossSection(
+        Code::Proton, Code::Oxygen, p4protonProj, p4oxygenTarg);
+
+    CORSIKA_LOG_INFO("Elab={:15.2f} GeV,Ecom={:15.2f} GeV, sig-p-p={:6.2f} mb, sig-p-O={:6.2f} mb",
+                      setE / 1_GeV, comEnn / 1_GeV, xs_prod_pp / 1_mb,
+                      xs_prod_pO / 1_mb);
+
+    out << i << " " << setE / 1_GeV << " " << comEnn / 1_GeV << " "
+        << xs_prod_pp / 1_mb << " " << xs_prod_pO / 1_mb << "\n";
+  }
+  out.close();
+}
+
+int main(int argc, char **argv) {
+
+  logging::set_level(logging::level::info);
+
+  if (argc != 2) {
+    std::cout
+        << "usage: check <interaction model> \n valid models are: sibyll, "
+           "epos, qgsII, pythia8"
+        << std::endl;
+    return 1;
+  }
+  std::string int_model_name = std::string(argv[1]);
+
+  if (int_model_name == "sibyll") {
+    RNGManager<>::getInstance().registerRandomStream("sibyll");
+    auto model = std::make_shared<corsika::sibyll::InteractionModel>(std::set<Code>{Code::Proton, Code::Oxygen});
+    calculate_cross_sections(model, int_model_name);    
+  } else if (int_model_name == "pythia8") {
+     RNGManager<>::getInstance().registerRandomStream("pythia");
+     auto model = std::make_shared<corsika::pythia8::InteractionModel>();
+     calculate_cross_sections(model, int_model_name);
+   
+   } else if (int_model_name == "epos") {
+     RNGManager<>::getInstance().registerRandomStream("epos");
+     auto model = std::make_shared<corsika::epos::InteractionModel>();
+     calculate_cross_sections(model, int_model_name);
+  }
+   else if (int_model_name == "qgsjet") {
+     RNGManager<>::getInstance().registerRandomStream("qgsjet");
+     auto model = std::make_shared<corsika::qgsjetII::InteractionModel>();
+     calculate_cross_sections(model, int_model_name);
+   }
+}