documentation/rtfd-requirements.txt

-sphinx_rtd_theme == 0.5.1
-sphinx == 3.5.3
-breathe == 4.27.0
-recommonmark == 0.7.1
+sphinx_rtd_theme >= 0.5.1
+sphinx >= 3.5.3
+recommonmark >= 0.7.1
+breathe >= 4.27.0
\ No newline at end of file

documentation/tracking.rst

+Tracking
+========
+
+Tracking is the process of moving a kinetic object from one point A to another point B. 
+Physical forces and laws should be applied explicitly or numerically in this process. 
+
+There are two main aspects of tracking:
+
+1. The actual proposal of a new point B, thus, a movement. 
+2. The description of the movement within a 3D geometry with boundaries and different volumina.
+
+Besides this, there furthermore is also the determination of the effects of the
+movement (e.g. magnetic fields, refractivity), or the accumulation of local 
+quantities (e.g. grammage). 
+For these purposes typically numerical integration routines are needed. 
+
+

examples/CMakeLists.txt

@@ -9,50 +9,98 @@ find_package (corsika CONFIG REQUIRED)
 # this is only for CORSIKA_REGISTER_EXAMPLE
 include ("${CMAKE_CURRENT_SOURCE_DIR}/corsikaExamples.cmake")
 
-add_executable (helix_example helix_example.cpp)
-target_link_libraries (helix_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (helix_example)
+# Put binaries in bin directory
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
 
-add_executable (geometry_example geometry_example.cpp)
-target_link_libraries (geometry_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (geometry_example)
+###################
+## cascade_examples
+###################
 
-add_executable (stack_example stack_example.cpp)
-target_link_libraries (stack_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (stack_example)
+add_executable (em_shower cascade_examples/em_shower.cpp)
+target_link_libraries (em_shower PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS 100 8472)
 
-add_executable (cascade_example cascade_example.cpp)
-target_link_libraries (cascade_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (cascade_example)
+if (WITH_FLUKA)
+  add_executable (mars cascade_examples/mars.cpp)
+  target_link_libraries (mars PRIVATE CORSIKA8::CORSIKA8)
+  message("FLUKA found, will make 'mars' example")
+else()
+  message("FLUKA not found, will not make 'mars' example")
+endif()
 
-add_executable (boundary_example boundary_example.cpp)
-target_link_libraries (boundary_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (boundary_example)
+add_executable (mc_conex cascade_examples/mc_conex.cpp)
+target_link_libraries (mc_conex  PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (mc_conex RUN_OPTIONS 4 2 10000.)
 
-add_executable (cascade_proton_example cascade_proton_example.cpp)
-target_link_libraries (cascade_proton_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (cascade_proton_example)
+add_executable (radio_em_shower cascade_examples/radio_em_shower.cpp)
+target_link_libraries (radio_em_shower PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (radio_em_shower RUN_OPTIONS 10 1121673489)
 
-add_executable (vertical_EAS vertical_EAS.cpp)
-target_link_libraries (vertical_EAS CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000. 1)
+if (WITH_FLUKA)
+  add_executable (water cascade_examples/water.cpp)
+  target_link_libraries (water PRIVATE CORSIKA8::CORSIKA8)
+  message("FLUKA found, will make 'water' example")
+else()
+  message("FLUKA not found, will not make 'water' example")
+endif()
 
-add_executable (stopping_power stopping_power.cpp)
-target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (stopping_power)
 
-add_executable (staticsequence_example staticsequence_example.cpp)
-target_link_libraries (staticsequence_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (staticsequence_example)
+#####################
+## framework_examples
+#####################
+
+add_executable (boundary_crossing framework_examples/boundary_crossing.cpp)
+target_link_libraries (boundary_crossing PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (boundary_crossing)
+
+add_executable (environment framework_examples/environment.cpp)
+target_link_libraries (environment PRIVATE CORSIKA8::CORSIKA8)
+
+add_executable (geometry framework_examples/geometry.cpp)
+target_link_libraries (geometry PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (geometry)
+
+add_executable (helix_trajectory framework_examples/helix_trajectory.cpp)
+target_link_libraries (helix_trajectory PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (helix_trajectory)
+
+add_executable (known_particles framework_examples/known_particles.cpp)
+target_link_libraries (known_particles PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (known_particles)
+
+add_executable (stack framework_examples/stack.cpp)
+target_link_libraries (stack PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (stack)
+
+add_executable (static_sequence framework_examples/static_sequence.cpp)
+target_link_libraries (static_sequence PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (static_sequence)
+
+
+###################
+## physics_examples
+###################
+
+add_executable (cross_section physics_examples/cross_section.cpp)
+target_link_libraries (cross_section PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (cross_section RUN_OPTIONS sibyll)
+
+add_executable (had_interactions physics_examples/had_interactions.cpp)
+target_link_libraries (had_interactions PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (had_interactions RUN_OPTIONS sibyll)
+
+add_executable (stopping_power physics_examples/stopping_power.cpp)
+target_link_libraries (stopping_power PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (stopping_power)
 
-add_executable (particle_list_example particle_list_example.cpp)
-target_link_libraries (particle_list_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (particle_list_example)
+add_executable (synchrotron_clover_leaf physics_examples/synchrotron_clover_leaf.cpp)
+target_link_libraries (synchrotron_clover_leaf PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_clover_leaf)
 
-add_executable (em_shower em_shower.cpp)
-target_link_libraries (em_shower CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS "100.")
+add_executable (synchrotron_test_C8tracking physics_examples/synchrotron_test_C8tracking.cpp)
+target_link_libraries (synchrotron_test_C8tracking PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_test_C8tracking)
 
-add_executable (hybrid_MC hybrid_MC.cpp)
-target_link_libraries (hybrid_MC  CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
+add_executable (synchrotron_test_manual_tracking physics_examples/synchrotron_test_manual_tracking.cpp)
+target_link_libraries (synchrotron_test_manual_tracking PRIVATE CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_test_manual_tracking)

examples/README.md

+# CORSIKA 8 Examples
+
+This directory contains several example scripts that can help you get started to simulate your own showers and/or to begin developing your own interfaces within the C8 framework.
+The sub-directories contain the several types of examples:
+
+ * **cascade_examples**: scripts to set up running full cascade simulations for various use cases. This is a good starting point if you want to simulate showers for your own experiments.
+ * **framework_examples**: show how to use and interact with the internal framework of C8
+ * **physics_examples**: demonstrate various tests of the particle interactions and/or cascades
+ 
\ No newline at end of file

examples/cascade_example.cpp

-/*
- * (c) Copyright 2018 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/framework/core/Cascade.hpp>
-#include <corsika/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Logging.hpp>
-
-#include <corsika/output/OutputManager.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/UniformMagneticField.hpp>
-
-#include <corsika/setup/SetupEnvironment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <corsika/modules/BetheBlochPDG.hpp>
-#include <corsika/modules/StackInspector.hpp>
-#include <corsika/modules/Sibyll.hpp>
-#include <corsika/modules/ParticleCut.hpp>
-#include <corsika/modules/TrackWriter.hpp>
-
-/*
-  NOTE, WARNING, ATTENTION
-
-  The .../Random.hpp implement the hooks of external modules to the C8 random
-  number generator. It has to occur excatly ONCE per linked
-  executable. If you include the header below multiple times and
-  link this togehter, it will fail.
- */
-#include <corsika/modules/sibyll/Random.hpp>
-#include <corsika/modules/urqmd/Random.hpp>
-
-#include <iostream>
-#include <limits>
-
-using namespace corsika;
-using namespace std;
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-
-  logging::set_level(logging::level::info);
-
-  std::cout << "cascade_example" << std::endl;
-
-  const LengthType height_atmosphere = 112.8_km;
-
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  RNGManager::getInstance().registerRandomStream("cascade");
-
-  // setup environment, geometry
-  setup::Environment env;
-  auto& universe = *(env.getUniverse());
-
-  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
-
-  auto world =
-      setup::Environment::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
-
-  using MyHomogeneousModel = MediumPropertyModel<
-      UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
-
-  // fraction of oxygen
-  float const fox = 0.20946;
-  auto const props = world->setModelProperties<MyHomogeneousModel>(
-      Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 0_T),
-      1_kg / (1_m * 1_m * 1_m),
-      NuclearComposition(std::vector<Code>{Code::Nitrogen, Code::Oxygen},
-                         std::vector<float>{1.f - fox, fox}));
-
-  auto innerMedium =
-      setup::Environment::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
-
-  innerMedium->setModelProperties(props);
-  world->addChild(std::move(innerMedium));
-  universe.addChild(std::move(world));
-
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.clear();
-  const Code beamCode = Code::Nucleus;
-  const int nuclA = 4;
-  const int nuclZ = int(nuclA / 2.15 + 0.7);
-  const HEPMassType mass = get_nucleus_mass(nuclA, nuclZ);
-  const HEPEnergyType E0 = nuclA * 1_TeV;
-  double theta = 0.;
-  double phi = 0.;
-
-  Point const injectionPos(
-      rootCS, 0_m, 0_m,
-      height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
-
-  OutputManager output("cascade_outputs");
-
-  ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
-
-  {
-    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-      return sqrt((Elab - m) * (Elab + m));
-    };
-    HEPMomentumType P0 = elab2plab(E0, mass);
-    auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
-      return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-                             -ptot * cos(theta));
-    };
-    auto const [px, py, pz] =
-        momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
-    auto plab = MomentumVector(rootCS, {px, py, pz});
-    cout << "input particle: " << beamCode << endl;
-    cout << "input angles: theta=" << theta << " phi=" << phi << endl;
-    cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
-    stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, nuclA, nuclZ));
-  }
-
-  // setup processes, decays and interactions
-  setup::Tracking tracking;
-  StackInspector<setup::Stack> stackInspect(1, true, E0);
-
-  RNGManager::getInstance().registerRandomStream("sibyll");
-  RNGManager::getInstance().registerRandomStream("pythia");
-  corsika::sibyll::Interaction sibyll;
-  corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  corsika::sibyll::Decay decay;
-  // cascade with only HE model ==> HE cut
-  ParticleCut cut(80_GeV, true, true);
-
-  TrackWriter trackWriter;
-  output.add("tracks", trackWriter); // register TrackWriter
-
-  BetheBlochPDG eLoss{showerAxis};
-
-  // assemble all processes into an ordered process list
-  auto sequence =
-      make_sequence(stackInspect, sibyll, sibyllNuc, decay, eLoss, cut, trackWriter);
-
-  // define air shower object, run simulation
-  Cascade EAS(env, tracking, sequence, output, stack);
-
-  EAS.run();
-
-  eLoss.printProfile(); // print longitudinal profile
-
-  cut.showResults();
-  const HEPEnergyType Efinal =
-      cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy();
-  cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-  cout << "total dEdX energy (GeV): " << eLoss.getTotal() / 1_GeV << endl
-       << "relative difference (%): " << eLoss.getTotal() / E0 * 100 << endl;
-  cut.reset();
-
-  output.endOfLibrary();
-}

examples/em_shower.cpp → examples/cascade_examples/em_shower.cpp

 /*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2022 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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/process/InteractionCounter.hpp>
 #include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/core/Logging.hpp>
 
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
 #include <corsika/output/OutputManager.hpp>
 
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
 #include <corsika/media/Environment.hpp>
 #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/media/ShowerAxis.hpp>
 #include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/ShowerAxis.hpp>
 #include <corsika/media/UniformMagneticField.hpp>
 
 #include <corsika/modules/LongitudinalProfile.hpp>
 #include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
 #include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/Sophia.hpp>
 #include <corsika/modules/TrackWriter.hpp>
-#include <corsika/modules/PROPOSAL.hpp>
 
 #include <corsika/setup/SetupStack.hpp>
 #include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
 
 #include <iomanip>
 #include <iostream>
@@ -42,88 +49,81 @@
 #include <string>
 #include <typeinfo>
 
-/*
-  NOTE, WARNING, ATTENTION
-
-  The .../Random.hpppp implement the hooks of external modules to the C8 random
-  number generator. It has to occur excatly ONCE per linked
-  executable. If you include the header below multiple times and
-  link this togehter, it will fail.
- */
-#include <corsika/modules/sibyll/Random.hpp>
-#include <corsika/modules/urqmd/Random.hpp>
-
 using namespace corsika;
 using namespace std;
 
-void registerRandomStreams() {
-  RNGManager::getInstance().registerRandomStream("cascade");
-  RNGManager::getInstance().registerRandomStream("proposal");
-  RNGManager::getInstance().seedAll();
+void registerRandomStreams(int seed) {
+  RNGManager<>::getInstance().registerRandomStream("cascade");
+  RNGManager<>::getInstance().registerRandomStream("proposal");
+  RNGManager<>::getInstance().registerRandomStream("sibyll");
+  RNGManager<>::getInstance().registerRandomStream("sophia");
+  if (seed == 0) {
+    std::random_device rd;
+    seed = rd();
+    CORSIKA_LOG_INFO("random seed (auto) {} ", seed);
+  } else {
+    CORSIKA_LOG_INFO("random seed {} ", seed);
+  }
+  RNGManager<>::getInstance().setSeed(seed);
 }
 
+using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using EnvType = Environment<EnvironmentInterface>;
 template <typename T>
 using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
 
 int main(int argc, char** argv) {
 
-  logging::set_level(logging::level::info);
+  logging::set_level(logging::level::warn);
 
-  if (argc != 2) {
-    std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
+  if (!(argc == 2 || argc == 3)) {
+    std::cerr << "usage: em_shower <energy/GeV> [seed]" << std::endl
+              << "seed = 0 for randomized seed" << std::endl;
     return 1;
   }
   feenableexcept(FE_INVALID);
+  int seed = 0;
+
+  if (argc >= 3) { seed = std::stoi(std::string(argv[2])); }
   // initialize random number sequence(s)
-  registerRandomStreams();
+  registerRandomStreams(seed);
 
   // setup environment, geometry
-  using EnvType = setup::Environment;
   EnvType env;
   CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
   Point const center{rootCS, 0_m, 0_m, 0_m};
-  auto builder = make_layered_spherical_atmosphere_builder<
-      setup::EnvironmentInterface, MyExtraEnv>::create(center,
-                                                       constants::EarthRadius::Mean,
-                                                       Medium::AirDry1Atm,
-                                                       Vector{rootCS, 0_T, 50_uT, 0_T});
-  builder.setNuclearComposition(
-      {{Code::Nitrogen, Code::Oxygen},
-       {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
-
-  builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
-  builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
-  builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
-  builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
-  builder.addLinearLayer(1e9_cm, 112.8_km);
-  builder.assemble(env);
 
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.clear();
+  // build a Linsley US Standard atmosphere into `env`
+  MagneticFieldVector bField{rootCS, 50_uT, 0_T, 0_T};
+  create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+      env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, bField);
+
+  std::unordered_map<Code, HEPEnergyType> energy_resolution = {
+      {Code::Electron, 5_MeV},
+      {Code::Positron, 5_MeV},
+      {Code::Photon, 5_MeV},
+  };
+  for (auto const& [pcode, energy] : energy_resolution)
+    set_energy_production_threshold(pcode, energy);
+
   const Code beamCode = Code::Electron;
   auto const mass = get_mass(beamCode);
   const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
   double theta = 0.;
   auto const thetaRad = theta / 180. * M_PI;
 
-  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-    return sqrt((Elab - m) * (Elab + m));
-  };
-  HEPMomentumType P0 = elab2plab(E0, mass);
+  HEPMomentumType P0 = calculate_momentum(E0, mass);
   auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
     return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
   };
 
   auto const [px, py, pz] = momentumComponents(thetaRad, P0);
   auto plab = MomentumVector(rootCS, {px, py, pz});
-  cout << "input particle: " << beamCode << endl;
-  cout << "input angles: theta=" << theta << endl;
-  cout << "input momentum: " << plab.getComponents() / 1_GeV
-       << ", norm = " << plab.getNorm() << endl;
 
-  auto const observationHeight = 1.4_km + builder.getEarthRadius();
-  auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+  auto const observationHeight = 0.0_km + constants::EarthRadius::Mean;
+  auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
   auto const t = -observationHeight * cos(thetaRad) +
                  sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
                       static_pow<2>(injectionHeight));
@@ -131,62 +131,86 @@ int main(int argc, char** argv) {
   Point const injectionPos =
       showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
 
-  std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
+                              false, 1000};
+  auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
 
-  stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+  CORSIKA_LOG_INFO("Primary particle:   {}", beamCode);
+  CORSIKA_LOG_INFO("Zenith angle:       {} (rad)", theta);
+  CORSIKA_LOG_INFO("Momentum:           {} (GeV)", plab.getComponents() / 1_GeV);
+  CORSIKA_LOG_INFO("Propagation dir:    {}", plab.getNorm());
+  CORSIKA_LOG_INFO("Injection point:    {}", injectionPos.getCoordinates());
+  CORSIKA_LOG_INFO("shower axis length: {} ",
+                   (showerCore - injectionPos).getNorm() * 1.02);
 
-  std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
-            << std::endl;
+  // setup processes, decays and interactions
+  EnergyLossWriter energyloss{showerAxis, dX};
+  ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV,
+                                                   100_GeV, true, energyloss);
+
+  corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
+                                      corsika::setup::C7trackedParticles);
+  corsika::sophia::InteractionModel sophia;
+  HEPEnergyType heThresholdNN = 80_GeV;
+  corsika::proposal::Interaction emCascade(
+      env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
+  corsika::proposal::ContinuousProcess<SubWriter<decltype(energyloss)>> emContinuous(
+      env, energyloss);
+
+  //  NOT possible right now, due to interface differenc in PROPOSAL
+  //  InteractionCounter emCascadeCounted(emCascade);
 
   OutputManager output("em_shower_outputs");
-  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
-
-  // setup processes, decays and interactions
 
-  ParticleCut cut(10_GeV, 10_GeV, 100_PeV, 100_PeV, true);
-  corsika::proposal::Interaction emCascade(env);
-  corsika::proposal::ContinuousProcess emContinuous(env);
-  InteractionCounter emCascadeCounted(emCascade);
+  output.add("energyloss", energyloss);
 
-  TrackWriter trackWriter;
-  output.add("tracks", trackWriter); // register TrackWriter
+  TrackWriter tracks;
+  output.add("tracks", tracks);
 
-  // long. profile; columns for photon, e+, e- still need to be added
-  LongitudinalProfile longprof{showerAxis};
+  LongitudinalWriter profile{showerAxis, dX};
+  output.add("profile", profile);
+  LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
 
   Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
-  ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
-                                    "particles.dat");
-  output.add("obsplane", observationLevel);
+  ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
+      obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
+  output.add("particles", observationLevel);
 
-  auto sequence = make_sequence(emCascadeCounted, emContinuous, longprof, cut,
-                                observationLevel, trackWriter);
+  PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
+  output.add("primary", primaryWriter);
+
+  auto sequence = make_sequence(emCascade, emContinuous, longprof, observationLevel, cut);
   // define air shower object, run simulation
-  setup::Tracking tracking;
+  TrackingType tracking;
+
+  output.startOfLibrary();
+
+  auto const primaryProperties = std::make_tuple(
+      beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+      plab.normalized(), injectionPos, 0_ns);
+
+  // setup particle stack, and add primary particle
+  StackType stack;
+  stack.clear();
+  stack.addParticle(primaryProperties);
+  primaryWriter.recordPrimary(primaryProperties);
+
   Cascade EAS(env, tracking, sequence, output, stack);
 
   // to fix the point of first interaction, uncomment the following two lines:
-  //  EAS.setNodes();
   //  EAS.forceInteraction();
 
   EAS.run();
 
-  cut.showResults();
-  emContinuous.showResults();
-  observationLevel.showResults();
-  const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
-                               cut.getEmEnergy() + emContinuous.getEnergyLost() +
-                               observationLevel.getEnergyGround();
-  cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
-  observationLevel.reset();
-  cut.reset();
-  emContinuous.reset();
-
-  auto const hists = emCascadeCounted.getHistogram();
-  save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
-  save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
-  longprof.save("longprof_emShower.txt");
+  HEPEnergyType const Efinal =
+      energyloss.getEnergyLost() + observationLevel.getEnergyGround();
+
+  CORSIKA_LOG_INFO(
+      "total energy budget (GeV): {}, "
+      "relative difference (%): {}",
+      Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
 
   output.endOfLibrary();
+
+  return EXIT_SUCCESS;
 }

examples/cascade_examples/mars.cpp

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/* clang-format off */
+// InteractionCounter used boost/histogram, which
+// fails if boost/type_traits have been included before. Thus, we have
+// to include it first...
+#include <corsika/framework/process/InteractionCounter.hpp>
+/* clang-format on */
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/process/InteractionCounter.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/output/OutputManager.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+#include <corsika/modules/BetheBlochPDG.hpp>
+#include <corsika/modules/LongitudinalProfile.hpp>
+#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
+#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/Pythia8.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/Sophia.hpp>
+#include <corsika/modules/StackInspector.hpp>
+#include <corsika/modules/TrackWriter.hpp>
+#include <corsika/modules/FLUKA.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
+
+#include <CLI/App.hpp>
+#include <CLI/Config.hpp>
+#include <CLI/Formatter.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+
+using namespace corsika;
+using namespace std;
+
+using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using EnvType = Environment<EnvironmentInterface>;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
+using Particle = StackType::particle_type;
+
+typedef decltype(1 * pascal) PressureType;
+typedef decltype(1 * degree_celsius) TemperatureType;
+
+class MarsAtmModel {
+public:
+  MarsAtmModel() = delete;
+  MarsAtmModel(PressureType a, InverseLengthType b, TemperatureType c,
+               decltype(1 * degree_celsius / 1_m) d)
+      : a_(a)
+      , b_(b)
+      , c_(c)
+      , d_(d) {}
+
+  MassDensityType operator()(LengthType height) const {
+    PressureType const pressure = a_ * exp(-b_ * height);
+    TemperatureType const temperature = -c_ - d_ * height + 273.1_K; // in K
+    constexpr decltype(square(1_m) / (square(1_s) * 1_K)) constant =
+        1000 * 0.1921 * square(1_m) / (square(1_s) * 1_K);
+    return pressure / (constant * temperature);
+  }
+
+private:
+  PressureType a_;
+  InverseLengthType b_;
+  TemperatureType c_;
+  decltype(1_K / 1_m) d_;
+};
+
+void registerRandomStreams(int seed) {
+  RNGManager<>::getInstance().registerRandomStream("cascade");
+  RNGManager<>::getInstance().registerRandomStream("qgsjet");
+  RNGManager<>::getInstance().registerRandomStream("sibyll");
+  RNGManager<>::getInstance().registerRandomStream("sophia");
+  RNGManager<>::getInstance().registerRandomStream("pythia");
+  RNGManager<>::getInstance().registerRandomStream("fluka");
+  RNGManager<>::getInstance().registerRandomStream("proposal");
+  if (seed == 0) {
+    std::random_device rd;
+    seed = rd();
+    CORSIKA_LOG_INFO("random seed (auto) {} ", seed);
+  } else {
+    CORSIKA_LOG_INFO("random seed {} ", seed);
+  }
+  RNGManager<>::getInstance().setSeed(seed);
+}
+
+template <typename T>
+using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
+
+int main(int argc, char** argv) {
+
+  // the main command line description
+  CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."};
+
+  // some options that we want to fill in
+  int A, Z, nevent = 0;
+
+  // the following section adds the options to the parser
+
+  // we start by definining a sub-group for the primary ID
+  auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary")
+                   ->check(CLI::Range(0, 26))
+                   ->group("Primary");
+  auto opt_A = app.add_option("-A", A, "Atomic mass number for primary")
+                   ->needs(opt_Z)
+                   ->check(CLI::Range(1, 58))
+                   ->group("Primary");
+  app.add_option("-p,--pdg", "PDG code for primary.")
+      ->excludes(opt_A)
+      ->excludes(opt_Z)
+      ->group("Primary");
+  // the remainding options
+  app.add_option("-E,--energy", "Primary energy in GeV")
+      ->required()
+      ->check(CLI::PositiveNumber)
+      ->group("Primary");
+  app.add_option("-z,--zenith", "Primary zenith angle (deg)")
+      ->required()
+      ->default_val(0.)
+      ->check(CLI::Range(0, 90))
+      ->group("Primary");
+  app.add_option("-a,--azimuth", "Primary azimuth angle (deg)")
+      ->default_val(0.)
+      ->check(CLI::Range(0, 360))
+      ->group("Primary");
+  app.add_option("-N,--nevent", nevent, "The number of events/showers to run.")
+      ->required()
+      ->check(CLI::PositiveNumber)
+      ->group("Library/Output");
+  app.add_option("-f,--filename", "Filename for output library.")
+      ->required()
+      ->default_val("corsika_library")
+      ->check(CLI::NonexistentPath)
+      ->group("Library/Output");
+  app.add_option("-s,--seed", "The random number seed.")
+      ->default_val(0)
+      ->check(CLI::NonNegativeNumber)
+      ->group("Misc.");
+  app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.")
+      ->default_val("info")
+      ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
+      ->group("Misc.");
+
+  // parse the command line options into the variables
+  CLI11_PARSE(app, argc, argv);
+
+  if (app.count("--verbosity")) {
+    auto const loglevel = app["--verbosity"]->as<std::string>();
+    if (loglevel == "warn") {
+      logging::set_level(logging::level::warn);
+    } else if (loglevel == "info") {
+      logging::set_level(logging::level::info);
+    } else if (loglevel == "debug") {
+      logging::set_level(logging::level::debug);
+    } else if (loglevel == "trace") {
+
+#ifndef C8_DEBUG
+      CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
+      return 1;
+#endif
+      logging::set_level(logging::level::trace);
+    }
+  }
+
+  // check that we got either PDG or A/Z
+  if (app.count("--pdg") == 0) {
+    if ((app.count("-A") == 0) || (app.count("-Z") == 0)) {
+      CORSIKA_LOG_ERROR("If --pdg is not provided, then both -A and -Z are required.");
+      return 1;
+    }
+  }
+
+  // initialize random number sequence(s)
+  registerRandomStreams(app["--seed"]->as<int>());
+
+  /* === START: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
+  EnvType env;
+  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+  Point const center{rootCS, 0_m, 0_m, 0_m};
+  LengthType const radiusMars = 3389.5_km;
+  auto builder =
+      make_layered_spherical_atmosphere_builder<EnvironmentInterface, MyExtraEnv>::create(
+          center,
+          radiusMars,                                   // Mars
+          Medium::AirDry1Atm,                           // Mars, close enough
+          MagneticFieldVector{rootCS, 0_T, 0_uT, 0_T}); // Mars
+
+  builder.setNuclearComposition(                   // Mars
+      {{Code::Carbon, Code::Oxygen,                // 95.97 CO2
+        Code::Nitrogen},                           // 1.89 N2 + 1.93 Argon + 0.146 O2
+       {0.9597 / 3, 0.9597 * 2 / 3, 1 - 0.9597}}); // values taken from AIRES manual
+
+  MarsAtmModel layer1(0.699e3 * pascal, 0.00009 / 1_m, 31.0 * degree_celsius,
+                      0.000998 * 1 * degree_celsius / 1_m);
+  MarsAtmModel layer2(0.699e3 * pascal, 0.00009 / 1_m, 23.4 * degree_celsius,
+                      0.00222 * 1 * degree_celsius / 1_m);
+
+  builder.addTabularLayer(layer1, 100, 100_m, 7_km);
+  builder.addTabularLayer(layer2, 300, 500_m, 100_km);
+  builder.addLinearLayer(1_g / square(1_cm), 1e9_cm, 112.8_km);
+  builder.assemble(env);
+  /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
+
+  ofstream atmout("mars.dat");
+  for (LengthType h = 0_m; h < 110_km; h += 10_m) {
+    Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h};
+    auto rho =
+        env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity(
+            ptest);
+    atmout << h / 1_m << " " << rho / 1_kg * cube(1_m) << "\n";
+  }
+  atmout.close();
+
+  /* === START: CONSTRUCT PRIMARY PARTICLE === */
+
+  // parse the primary ID as a PDG or A/Z code
+  Code beamCode;
+
+  // check if we want to use a PDG code instead
+  if (app.count("--pdg") > 0) {
+    beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
+  } else {
+    // check manually for proton and neutrons
+    if ((A == 0) && (Z == 1)) beamCode = Code::Proton;
+    if ((A == 1) && (Z == 1)) beamCode = Code::Neutron;
+  }
+  HEPEnergyType const mass = get_mass(beamCode);
+
+  // particle energy
+  HEPEnergyType const E0 = 1_GeV * app["--energy"]->as<double>();
+
+  // direction of the shower in (theta, phi) space
+  auto const thetaRad = app["--zenith"]->as<double>() / 180. * M_PI;
+  auto const phiRad = app["--azimuth"]->as<double>() / 180. * M_PI;
+
+  // convert Elab to Plab
+  HEPMomentumType P0 = sqrt((E0 - mass) * (E0 + mass));
+
+  // convert the momentum to the zenith and azimuth angle of the primary
+  auto const [px, py, pz] =
+      std::make_tuple(P0 * sin(thetaRad) * cos(phiRad), P0 * sin(thetaRad) * sin(phiRad),
+                      -P0 * cos(thetaRad));
+  auto plab = MomentumVector(rootCS, {px, py, pz});
+  /* === END: CONSTRUCT PRIMARY PARTICLE === */
+
+  /* === START: CONSTRUCT GEOMETRY === */
+  auto const observationHeight = 0_km + builder.getPlanetRadius();
+  auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
+  auto const t = -observationHeight * cos(thetaRad) +
+                 sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+                      static_pow<2>(injectionHeight));
+  Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+  Point const injectionPos =
+      showerCore + DirectionVector{rootCS,
+                                   {-sin(thetaRad) * cos(phiRad),
+                                    -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
+                       t;
+
+  // we make the axis much longer than the inj-core distance since the
+  // profile will go beyond the core, depending on zenith angle
+  /* === END: CONSTRUCT GEOMETRY === */
+
+  // create the output manager that we then register outputs with
+  OutputManager output(app["--filename"]->as<std::string>());
+
+  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
+  auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
+
+  EnergyLossWriter dEdX{showerAxis, dX};
+  output.add("energyloss", dEdX);
+
+  HEPEnergyType const emcut = 1_GeV;
+  HEPEnergyType const hadcut = 1_GeV;
+  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, hadcut, true,
+                                             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}));
+  set_energy_production_threshold(Code::Positron, std::min({emcut, hadcut}));
+  set_energy_production_threshold(Code::Photon, std::min({emcut, hadcut}));
+  set_energy_production_threshold(Code::MuMinus, std::min({emcut, hadcut}));
+  set_energy_production_threshold(Code::MuPlus, std::min({emcut, hadcut}));
+  set_energy_production_threshold(Code::TauMinus, std::min({emcut, hadcut}));
+  set_energy_production_threshold(Code::TauPlus, std::min({emcut, hadcut}));
+
+  /* === START: SETUP PROCESS LIST === */
+  auto const all_elements = corsika::get_all_elements_in_universe(env);
+  corsika::sibyll::Interaction sibyll(all_elements, corsika::setup::C7trackedParticles);
+  InteractionCounter sibyllCounted(sibyll);
+
+  corsika::pythia8::Decay decayPythia;
+
+  // energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
+  // interactions and the hadronic photon model in proposal
+  HEPEnergyType heHadronModelThreshold = 63.1_GeV;
+
+  corsika::sophia::InteractionModel sophia;
+  corsika::proposal::Interaction emCascade(
+      env, sophia, sibyll.getHadronInteractionModel(), heHadronModelThreshold);
+
+  // use BetheBlochPDG for hadronic continuous losses, and proposal otherwise
+  corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuousProposal(
+      env, dEdX);
+  BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuousBethe{dEdX};
+  struct EMHadronSwitch {
+    EMHadronSwitch() = default;
+    bool operator()(const Particle& p) const { return is_hadron(p.getPID()); }
+  };
+  auto emContinuous =
+      make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
+
+  LongitudinalWriter longprof{showerAxis, dX};
+  output.add("profile", longprof);
+  LongitudinalProfile<SubWriter<decltype(longprof)>> profile{longprof};
+
+  corsika::fluka::Interaction leIntModel{all_elements};
+  InteractionCounter leIntCounted{leIntModel};
+  StackInspector<StackType> stackInspect(5000, false, E0);
+
+  // assemble all processes into an ordered process list
+  struct EnergySwitch {
+    HEPEnergyType cutE_;
+    EnergySwitch(HEPEnergyType cutE)
+        : cutE_(cutE) {}
+    bool operator()(Particle const& p) const { return (p.getKineticEnergy() < cutE_); }
+  };
+  auto hadronSequence =
+      make_select(EnergySwitch(heHadronModelThreshold), leIntCounted, sibyllCounted);
+
+  // track writer
+  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter
+
+  // observation plane
+  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+  ObservationPlane<TrackingType> observationLevel(obsPlane,
+                                                  DirectionVector(rootCS, {1., 0., 0.}));
+  // register the observation plane with the output
+  output.add("particles", observationLevel);
+
+  PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
+  output.add("primary", primaryWriter);
+
+  // assemble the final process sequence
+  auto sequence =
+      make_sequence(stackInspect, hadronSequence, decayPythia, emCascade, emContinuous,
+                    trackWriter, profile, observationLevel, cut);
+  /* === END: SETUP PROCESS LIST === */
+
+  // create the cascade object using the default stack and tracking
+  // implementation
+  TrackingType tracking;
+  StackType stack;
+  Cascade EAS(env, tracking, sequence, output, stack);
+
+  // print our primary parameters all in one place
+  if (app["--pdg"]->count() > 0) {
+    CORSIKA_LOG_INFO("Primary PDG ID:     {}", app["--pdg"]->as<int>());
+  } else {
+    CORSIKA_LOG_INFO("Primary Z/A:        {}/{}", Z, A);
+  }
+  CORSIKA_LOG_INFO("Primary Energy:     {}", E0);
+  CORSIKA_LOG_INFO("Primary Momentum:   {}", P0);
+  CORSIKA_LOG_INFO("Primary Direction:  {}", plab.getNorm());
+  CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
+  CORSIKA_LOG_INFO("Shower Axis Length: {}", (showerCore - injectionPos).getNorm() * 1.2);
+
+  // trigger the output manager to open the library for writing
+  output.startOfLibrary();
+
+  // loop over each shower
+  for (int i_shower = 1; i_shower < nevent + 1; i_shower++) {
+    CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, nevent);
+
+    // trigger the start of the outputs for this shower
+    output.startOfShower();
+
+    // directory for outputs
+    string const outdir(app["--filename"]->as<std::string>());
+    string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
+    string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
+
+    // setup particle stack, and add primary particle
+    stack.clear();
+
+    // add the desired particle to the stack
+    auto const primaryProperties = std::make_tuple(
+        beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+        plab.normalized(), injectionPos, 0_ns);
+    stack.addParticle(primaryProperties);
+
+    primaryWriter.recordPrimary(primaryProperties);
+
+    // run the shower
+    EAS.run();
+
+    HEPEnergyType const Efinal =
+        dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+
+    CORSIKA_LOG_INFO(
+        "total energy budget (GeV): {}, "
+        "relative difference (%): {}",
+        Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
+
+    // trigger the output manager to save this shower to disk
+    output.endOfShower();
+  }
+
+  // and finalize the output on disk
+  output.endOfLibrary();
+}