CMakeLists.txt

@@ -262,6 +262,13 @@ add_subdirectory (modules/conex)
 add_subdirectory (modules/epos)
 add_subdirectory (modules/fluka)
 
+#
+#+++++++++++++++++++++++++++++++
+# standard applications
+#
+
+add_subdirectory(applications)
+
 #+++++++++++++++++++++++++++++++
 # unit testing
 #

README.md

@@ -102,8 +102,6 @@ make install
 
 ### FLUKA support
 
-Warning: may only work when the next version of FLUKA is released (as of 2023.06.15)
-
 For legal reasons we do not distribute/bundle FLUKA together with CORSIKA 8.
 If you want to use FLUKA as low-energy hadronic interaction model, you have to download
 it separately from (http://www.fluka.org/), which requires registering there as FLUKA user.
@@ -142,22 +140,40 @@ make install
 To run the Unit Tests, just type `ctest` in your build area.
 
 
-## Running examples
+## Running applications and examples
+
+### Standard applications
+
+Applications for standard use-cases are located in the `applications` directory.
+These are example scripts that can be used directly or slightly modified for your use case.
+See [applications/README.md] for more.
+The applications are compiled automatically after running `make` and will appear your `corsika-build/bin` directory.
+After running `make install` the binaries will also be copied into your `corsika-install/bin` directory as well.
+
+
+For example, from inside your `corsika-install/bin` directory, run 
+```shell
+c8_air_shower --pdg 2212 -E 1e5 -f my_shower
+```
+This will run a vertical 100 TeV proton shower and will create and put the output into `./my_shower`.
+
 
 ### Building the examples
 
-From your top corsika build directory, (the one that includes `corsika-build` and `corsika-install`) type
+Unlike the applications, the examples must be compiled as a second step.
+From your top corsika directory, (the one that includes `corsika-build` and `corsika-install`) run
 ```shell
 cmake -Dcorsika_DIR=$PWD/corsika-build -S ./corsika/examples -B ./corsika-build-examples
 cd corsika-build-examples
 make -j4 #The number should match the number of available cores on your machine
 ```
 
-From any directory, run the program `corsika-build-examples/bin/corsika`. As an example, you can run with the following flags:
+You can run the examples by using the binaries in `corsika-build-examples/bin/`.
+For example:
 ```shell
-corsika-build-examples/bin/corsika  --pdg 2212 -E 1e5 -f my_shower
+corsika-build-examples/bin/known_particles
 ```
-This will run a vertical 100 TeV proton shower and will create and put the output into `./my_shower`.
+This will print out all of the particles that are known by CORSIKA.
 
 
 ### Generating doxygen documentation

applications/CMakeLists.txt

+SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+
+if(WITH_FLUKA)
+  message("compiling c8_air_shower.cpp with FLUKA")
+else()
+  message("compiling c8_air_shower.cpp with UrQMD")
+endif()
+add_executable (c8_air_shower c8_air_shower.cpp)
+target_link_libraries (c8_air_shower CORSIKA8)
+
+install (
+  TARGETS c8_air_shower DESTINATION bin
+)

applications/README.md

+# CORSIKA 8 Applications
+
+This directory contains standard applications which are typical for astroparticle physics solutions.
+They are "physics-complete" and are suitable for generating simulations that can be used in publications.
+
+For example, `c8_air_shower.cpp` would be a similar binary to what would be built by CORSIKA 7 and will simulate
+air showers in a curved atmosphere.

examples/corsika.cpp → applications/c8_air_shower.cpp

@@ -20,7 +20,6 @@
 #include <corsika/framework/geometry/Plane.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
 #include <corsika/framework/process/DynamicInteractionProcess.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>
@@ -35,7 +34,6 @@
 
 #include <corsika/media/CORSIKA7Atmospheres.hpp>
 #include <corsika/media/Environment.hpp>
-#include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/GeomagneticModel.hpp>
 #include <corsika/media/GladstoneDaleRefractiveIndex.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
@@ -50,7 +48,6 @@
 #include <corsika/modules/Epos.hpp>
 #include <corsika/modules/LongitudinalProfile.hpp>
 #include <corsika/modules/ObservationPlane.hpp>
-#include <corsika/modules/OnShellCheck.hpp>
 #include <corsika/modules/PROPOSAL.hpp>
 #include <corsika/modules/ParticleCut.hpp>
 #include <corsika/modules/Pythia8.hpp>
@@ -94,8 +91,17 @@ using namespace std;
 using EnvironmentInterface =
     IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
 using EnvType = Environment<EnvironmentInterface>;
-
-using Particle = setup::Stack<EnvType>::particle_type;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
+using Particle = StackType::particle_type;
+
+//
+// This is the main example script which runs EAS with fairly standard settings
+// w.r.t. what was implemented in CORSIKA 7. Users may want to change some of the
+// specifics (observation altitude, magnetic field, energy cuts, etc.), but this
+// example is the most physics-complete one and should be used for full simulations
+// of particle cascades in air
+//
 
 long registerRandomStreams(long seed) {
   RNGManager<>::getInstance().registerRandomStream("cascade");
@@ -111,9 +117,9 @@ long registerRandomStreams(long seed) {
   if (seed == 0) {
     std::random_device rd;
     seed = rd();
-    std::cout << "random seed (auto)  " << seed << std::endl;
+    CORSIKA_LOG_INFO("random seed (auto) {}", seed);
   } else {
-    std::cout << "random seed " << seed << std::endl;
+    CORSIKA_LOG_INFO("random seed {}", seed);
   }
   RNGManager<>::getInstance().setSeed(seed);
   return seed;
@@ -346,6 +352,8 @@ int main(int argc, char** argv) {
   // we make the axis much longer than the inj-core distance since the
   // profile will go beyond the core, depending on zenith angle
   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
+  uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
   /* === END: CONSTRUCT GEOMETRY === */
 
   double const emthinfrac = app["--emthin"]->as<double>();
@@ -364,22 +372,22 @@ int main(int argc, char** argv) {
   OutputManager output(app["--filename"]->as<std::string>(), seed, args.str(), outputDir);
 
   // register energy losses as output
-  EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
+  EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
   output.add("energyloss", dEdX);
 
-  DynamicInteractionProcess<setup::Stack<EnvType>> heModel;
+  DynamicInteractionProcess<StackType> heModel;
 
   // have SIBYLL always for PROPOSAL photo-hadronic interactions
   auto sibyll = std::make_shared<corsika::sibyll::Interaction>(env);
 
   if (auto const modelStr = app["--hadronModel"]->as<std::string>();
       modelStr == "SIBYLL-2.3d") {
-    heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{sibyll};
+    heModel = DynamicInteractionProcess<StackType>{sibyll};
   } else if (modelStr == "QGSJet-II.04") {
-    heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+    heModel = DynamicInteractionProcess<StackType>{
         std::make_shared<corsika::qgsjetII::Interaction>()};
   } else if (modelStr == "EPOS-LHC") {
-    heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+    heModel = DynamicInteractionProcess<StackType>{
         std::make_shared<corsika::epos::Interaction>()};
   } else {
     CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr);
@@ -426,7 +434,7 @@ int main(int argc, char** argv) {
   auto emContinuous =
       make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
 
-  LongitudinalWriter profile{showerAxis, 200, 10_g / square(1_cm)};
+  LongitudinalWriter profile{showerAxis, nAxisBins, dX};
   output.add("profile", profile);
   LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
 
@@ -438,7 +446,7 @@ int main(int argc, char** argv) {
 #endif
   InteractionCounter leIntCounted{leIntModel};
 
-  StackInspector<setup::Stack<EnvType>> stackInspect(10000, false, E0);
+  StackInspector<StackType> stackInspect(10000, false, E0);
 
   // assemble all processes into an ordered process list
   struct EnergySwitch {
@@ -452,30 +460,20 @@ int main(int argc, char** argv) {
 
   // observation plane
   Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
-  ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+  ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
       obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
       true,   // plane should "absorb" particles
       false}; // do not print z-coordinate
   // register ground particle output
   output.add("particles", observationLevel);
 
-  PrimaryWriter<setup::Tracking, ParticleWriterParquet> primaryWriter(observationLevel);
+  PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
   output.add("primary", primaryWriter);
 
   int ring_number{app["--ring"]->as<int>()};
-  std::cout << "Ring number : " << ring_number << std::endl;
   auto const radius_{ring_number * 25_m};
-  //  std::cout << "Radius = " << radius_ << std::endl;
   const int rr_ = static_cast<int>(radius_ / 1_m);
 
-  // if (ring_number == 0) {
-  //     // assemble the final process sequence without radio
-  //     auto sequence = make_sequence(stackInspect, hadronSequence,
-  //     decaySequence, emCascade,
-  //                                   emContinuous, longprof, observationLevel,
-  //                                   thinning, cut);
-  // } else {
-
   // Radio antennas and relevant information
   // the antenna time variables
   const TimeType duration_{4e-7_s};
@@ -491,14 +489,11 @@ int main(int argc, char** argv) {
   auto const injectionPosY_{injectionPos.getCoordinates().getY()};
   auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
   auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
-  std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
 
   if (ring_number != 0) {
     // setup CoREAS antennas - use the for loop for star shape pattern
     for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
       auto phiRad_1 = phi_1 / 180. * M_PI;
-      // auto phi_1 = 0;
-      // auto phiRad_1 = phi_1 / 180. * M_PI;
       auto const point_1{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_1),
                                showerCoreY_ + radius_ * sin(phiRad_1),
                                constants::EarthRadius::Mean)};
@@ -514,8 +509,6 @@ int main(int argc, char** argv) {
     // setup ZHS antennas - use the for loop for star shape pattern
     for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
       auto phiRad_ = phi_ / 180. * M_PI;
-      // auto phi_ = 0; phi_;
-      // auto phiRad_ = phi_ / 180. * M_PI;
       auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
                               showerCoreY_ + radius_ * sin(phiRad_),
                               constants::EarthRadius::Mean)};
@@ -557,18 +550,19 @@ int main(int argc, char** argv) {
 
   // create the cascade object using the default stack and tracking
   // implementation
-  setup::Tracking tracking;
-  setup::Stack<EnvType> stack;
+  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>());
+    CORSIKA_LOG_INFO("Primary PDG ID:     {}", app["--pdg"]->as<int>());
   } else {
-    CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
+    CORSIKA_LOG_INFO("Primary Z/A:        {}/{}", Z, A);
   }
-  CORSIKA_LOG_INFO("Primary Energy: {}", E0);
-  CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+  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);
 

corsika/detail/framework/process/SwitchProcessSequence.inl

@@ -77,10 +77,10 @@ namespace corsika {
   template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
             int IndexProcess1, int IndexProcess2>
   template <typename TParticle>
-  inline ProcessReturn SwitchProcessSequence<
-      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
-      IndexProcess2>::doContinuous(Step<TParticle>& step,
-                                   ContinuousProcessIndex const idLimit) {
+  inline ProcessReturn SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+                                             IndexProcess1, IndexProcess2>::
+      doContinuous(Step<TParticle>& step,
+                   [[maybe_unused]] ContinuousProcessIndex const idLimit) {
     if (select_(step.getParticlePre())) {
       if constexpr (process1_type::is_process_sequence) {
         return A_.doContinuous(step, idLimit);
@@ -218,10 +218,10 @@ namespace corsika {
   template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
             int IndexProcess1, int IndexProcess2>
   template <typename TParticle>
-  CrossSectionType SwitchProcessSequence<
-      TCondition, TSequence, USequence, IndexStart, IndexProcess1,
-      IndexProcess2>::getCrossSection(TParticle const& projectile, Code const targetId,
-                                      FourMomentum const& targetP4) const {
+  CrossSectionType SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+                                         IndexProcess1, IndexProcess2>::
+      getCrossSection(TParticle const& projectile, [[maybe_unused]] Code const targetId,
+                      [[maybe_unused]] FourMomentum const& targetP4) const {
 
     if (select_(projectile)) {
       if constexpr (is_interaction_process_v<process1_type>) {

corsika/detail/modules/energy_loss/BetheBlochPDG.inl

@@ -162,10 +162,8 @@ namespace corsika {
     HEPEnergyType const dE = getTotalEnergyLoss(step.getParticlePre(), dX);
     //    if (dE > HEPEnergyType::zero())
     //      dE = -dE;
-    [[maybe_unused]] const auto Ekin = step.getEkinPre();
-    auto EkinNew = Ekin + dE;
     CORSIKA_LOG_TRACE("EnergyLoss  dE={} MeV, Ekin={} GeV, EkinNew={} GeV", dE / 1_MeV,
-                      Ekin / 1_GeV, EkinNew / 1_GeV);
+                      step.getEkinPre() / 1_GeV, (step.getEkinPre() + dE) / 1_GeV);
     step.add_dEkin(dE);
 
     // also send to output

corsika/detail/modules/proposal/ContinuousProcess.inl

@@ -48,8 +48,8 @@ namespace corsika::proposal {
   template <typename TEnvironment, typename... TOutputArgs>
   inline ContinuousProcess<TOutput>::ContinuousProcess(TEnvironment const& _env,
                                                        TOutputArgs&&... args)
-      : TOutput(args...)
-      , ProposalProcessBase(_env) {}
+      : ProposalProcessBase(_env)
+      , TOutput(args...) {}
 
   template <typename TOutput>
   template <typename TParticle>

corsika/detail/modules/proposal/InteractionModel.inl

@@ -56,7 +56,8 @@ namespace corsika::proposal {
   template <typename TStackView>
   inline ProcessReturn
   InteractionModel<THadronicLEModel, THadronicHEModel>::doInteraction(
-      TStackView& view, Code const projectileId, FourMomentum const& projectileP4) {
+      TStackView& view, Code const projectileId,
+      [[maybe_unused]] FourMomentum const& projectileP4) {
 
     auto const projectile = view.getProjectile();
 

corsika/detail/modules/proposal/ProposalProcessBase.inl

@@ -46,7 +46,7 @@ namespace corsika::proposal {
     }
 
     return lowest_table_value;
-  };
+  }
 
   template <typename TEnvironment>
   inline ProposalProcessBase::ProposalProcessBase(TEnvironment const& _env) {

corsika/detail/modules/radio/RadioProcess.inl

@@ -57,7 +57,8 @@ namespace corsika {
   template <typename Particle, typename Track>
   inline LengthType
   RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::getMaxStepLength(
-      const Particle& vParticle, const Track& vTrack) const {
+      [[maybe_unused]] const Particle& vParticle,
+      [[maybe_unused]] const Track& vTrack) const {
     return meter * std::numeric_limits<double>::infinity();
   }
 
@@ -176,4 +177,4 @@ namespace corsika {
     return config;
   }
 
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika

corsika/modules/Random.hpp

@@ -23,7 +23,7 @@ namespace corsika {
     inline void rng_func(corsika::default_prng_type& rng, double* dest, std::size_t N) {
       std::uniform_real_distribution<double> udist(0.0, 1.0);
       std::generate(dest, std::next(dest, N), std::bind(udist, std::ref(rng)));
-    };
+    }
   } // namespace detail
 
   inline void connect_random_stream(corsika::default_prng_type& rng,

corsika/modules/writers/EnergyLossWriter.hpp

@@ -156,7 +156,7 @@ namespace corsika {
     /**
      * Return the configuration of this output.
      */
-    YAML::Node getConfig() const;
+    YAML::Node getConfig() const override;
 
   private:
     ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage

corsika/modules/writers/LongitudinalWriter.hpp

@@ -122,12 +122,12 @@ namespace corsika {
     /**
      * Return a summary.
      */
-    YAML::Node getSummary() const;
+    YAML::Node getSummary() const override;
 
     /**
      * Return the configuration of this output.
      */
-    YAML::Node getConfig() const;
+    YAML::Node getConfig() const override;
 
     number_profile::ProfileData const& getProfile(
         number_profile::ProfileIndex index) const {

corsika/modules/writers/WriterOff.hpp

@@ -43,7 +43,7 @@ namespace corsika {
     template <typename... TArgs>
     void write(TArgs&&...) {}
 
-    virtual YAML::Node getConfig() const { return YAML::Node(); }
+    virtual YAML::Node getConfig() const override { return YAML::Node(); }
 
   }; // class WriterOff
 

examples/CMakeLists.txt

@@ -9,88 +9,88 @@ 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)
 
-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 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 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  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 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.)
+if (WITH_FLUKA)
+  add_executable (water cascade_examples/water.cpp)
+  target_link_libraries (water 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 (particle_list_example particle_list_example.cpp)
-target_link_libraries (particle_list_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (particle_list_example)
+add_executable (boundary_crossing framework_examples/boundary_crossing.cpp)
+target_link_libraries (boundary_crossing CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (boundary_crossing)
 
-add_executable (em_shower em_shower.cpp)
-target_link_libraries (em_shower CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS 100 8472)
+add_executable (environment framework_examples/environment.cpp)
+target_link_libraries (environment CORSIKA8::CORSIKA8)
 
-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 (geometry framework_examples/geometry.cpp)
+target_link_libraries (geometry CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (geometry)
 
+add_executable (helix_trajectory framework_examples/helix_trajectory.cpp)
+target_link_libraries (helix_trajectory CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (helix_trajectory)
 
-add_executable (corsika corsika.cpp)
-# for ICRC 2023
-option(WITH_FLUKA "use FLUKA instead of UrQMD in corsika.cpp")
-if (WITH_FLUKA)
-  target_compile_definitions(corsika PRIVATE WITH_FLUKA)
-  message("compiling corsika.cpp with FLUKA")
-else()
-  message("compiling corsika.cpp with UrQMD")
-endif()
-target_link_libraries (corsika CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (corsika RUN_OPTIONS --energy=10 --zenith=10 --nevent=2 --filename=corsika_test --pdg=2212)
+add_executable (known_particles framework_examples/known_particles.cpp)
+target_link_libraries (known_particles CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (known_particles)
 
-add_executable (mars mars.cpp)
-target_link_libraries (mars CORSIKA8::CORSIKA8)
+add_executable (stack framework_examples/stack.cpp)
+target_link_libraries (stack CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (stack)
 
-add_executable (water water.cpp)
-target_link_libraries (water CORSIKA8::CORSIKA8)
+add_executable (static_sequence framework_examples/static_sequence.cpp)
+target_link_libraries (static_sequence CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (static_sequence)
 
-add_executable (environment environment.cpp)
 
-target_link_libraries (environment CORSIKA8::CORSIKA8)
-add_executable (synchrotron_test_manual_tracking synchrotron_test_manual_tracking.cpp)
-target_link_libraries (synchrotron_test_manual_tracking CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (synchrotron_test_manual_tracking)
+###################
+## physics_examples
+###################
 
-add_executable (synchrotron_test_C8tracking synchrotron_test_C8tracking.cpp)
+add_executable (stopping_power physics_examples/stopping_power.cpp)
+target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (stopping_power)
+
+add_executable (synchrotron_clover_leaf physics_examples/synchrotron_clover_leaf.cpp)
+target_link_libraries (synchrotron_clover_leaf CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_clover_leaf)
+
+add_executable (synchrotron_test_C8tracking physics_examples/synchrotron_test_C8tracking.cpp)
 target_link_libraries (synchrotron_test_C8tracking CORSIKA8::CORSIKA8)
 CORSIKA_REGISTER_EXAMPLE (synchrotron_test_C8tracking)
 
-add_executable (clover_leaf clover_leaf.cpp)
-target_link_libraries (clover_leaf CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (clover_leaf)
-
-add_executable (radio_em_shower radio_em_shower.cpp)
-target_link_libraries (radio_em_shower CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (radio_em_shower RUN_OPTIONS 10 1121673489)
+add_executable (synchrotron_test_manual_tracking physics_examples/synchrotron_test_manual_tracking.cpp)
+target_link_libraries (synchrotron_test_manual_tracking 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/process/ProcessSequence.hpp>
-#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/core/EnergyMomentumOperations.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-
-#include <corsika/output/OutputManager.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriter.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/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>
-
-#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::warn);
-
-  CORSIKA_LOG_INFO("cascade_example");
-
-  LengthType const height_atmosphere = 112.8_km;
-
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  RNGManager<>::getInstance().registerRandomStream("cascade");
-
-  // setup environment, geometry
-  using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
-  using EnvType = Environment<EnvironmentInterface>;
-  EnvType env;
-  auto& universe = *(env.getUniverse());
-
-  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
-
-  auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
-
-  using MyHomogeneousModel =
-      MediumPropertyModel<UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>;
-
-  // fraction of oxygen
-  double 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({Code::Nitrogen, Code::Oxygen}, {1. - fox, fox}));
-
-  auto innerMedium = EnvType::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<EnvType> stack;
-  stack.clear();
-  const int nuclA = 4;
-  const int nuclZ = int(nuclA / 2.15 + 0.7);
-  const Code beamCode = get_nucleus_code(nuclA, nuclZ);
-  const HEPMassType mass = get_nucleus_mass(beamCode);
-  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");
-
-  theta *= M_PI / 180.;
-  phi *= M_PI / 180.;
-  DirectionVector const direction(
-      rootCS, {sin(theta) * cos(phi), sin(theta) * sin(phi), -cos(theta)});
-
-  ShowerAxis const showerAxis{injectionPos, direction * 100_km, env};
-  EnergyLossWriter dEdX{showerAxis};
-  output.add("energyloss", dEdX);
-
-  {
-    HEPMomentumType const P0 = calculate_momentum(E0, mass);
-    auto plab = direction * P0;
-    CORSIKA_LOG_INFO("input particle: {}", beamCode);
-    CORSIKA_LOG_INFO("input angles: theta={} phi={}", theta, phi);
-    CORSIKA_LOG_INFO("input momentum: {}", plab.getComponents() / 1_GeV);
-    stack.addParticle(std::make_tuple(
-        beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)), direction,
-        injectionPos, 0_ns));
-  }
-
-  // setup processes, decays and interactions
-  setup::Tracking tracking;
-  StackInspector<setup::Stack<EnvType>> stackInspect(100, true, E0);
-
-  RNGManager<>::getInstance().registerRandomStream("sibyll");
-
-  corsika::sibyll::Interaction sibyll{env};
-  corsika::sibyll::Decay decay;
-
-  // cascade with only HE model ==> HE cut
-  ParticleCut<SubWriter<decltype(dEdX)>> cut(80_GeV, true, dEdX);
-  BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
-
-  TrackWriter trackWriter;
-  output.add("tracks", trackWriter); // register TrackWriter
-
-  // assemble all processes into an ordered process list
-  auto sequence = make_sequence(stackInspect, sibyll, decay, eLoss, trackWriter, cut);
-
-  // define air shower object, run simulation
-  Cascade EAS(env, tracking, sequence, output, stack);
-
-  output.startOfLibrary();
-  EAS.run();
-  output.endOfLibrary();
-
-  const HEPEnergyType Efinal = dEdX.getEnergyLost();
-  CORSIKA_LOG_INFO(
-      "\n"
-      "total cut energy (GeV) : {}\n"
-      "relative difference (%): {}\n"
-      "total dEdX energy (GeV): {}\n"
-      "relative difference (%): {}\n",
-      Efinal / 1_GeV, (Efinal / E0 - 1) * 100, dEdX.getEnergyLost() / 1_GeV,
-      dEdX.getEnergyLost() / E0 * 100);
-}

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

@@ -18,10 +18,10 @@
 #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>
 
@@ -29,11 +29,9 @@
 #include <corsika/media/Environment.hpp>
 #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
 #include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/NuclearComposition.hpp>
 #include <corsika/media/ShowerAxis.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>
@@ -62,13 +60,19 @@ void registerRandomStreams(int seed) {
   if (seed == 0) {
     std::random_device rd;
     seed = rd();
-    cout << "new random seed (auto) " << seed << endl;
+    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) {
 
@@ -87,29 +91,23 @@ int main(int argc, char** argv) {
   registerRandomStreams(seed);
 
   // setup environment, geometry
-  using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
-  using EnvType = Environment<EnvironmentInterface>;
   EnvType env;
   CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
   Point const center{rootCS, 0_m, 0_m, 0_m};
 
   // 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,
-      MagneticFieldVector{rootCS, 20.4_uT, 0_T, -43.23_uT});
+      env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, bField);
 
   std::unordered_map<Code, HEPEnergyType> energy_resolution = {
-      {Code::Electron, 2_MeV},
-      {Code::Positron, 2_MeV},
-      {Code::Photon, 2_MeV},
+      {Code::Electron, 5_MeV},
+      {Code::Positron, 5_MeV},
+      {Code::Photon, 5_MeV},
   };
   for (auto [pcode, energy] : energy_resolution)
     set_energy_production_threshold(pcode, energy);
 
-  // setup particle stack, and add primary particle
-  setup::Stack<EnvType> stack;
-  stack.clear();
-
   const Code beamCode = Code::Electron;
   auto const mass = get_mass(beamCode);
   const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
@@ -123,10 +121,6 @@ int main(int argc, char** argv) {
 
   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 = 0.0_km + constants::EarthRadius::Mean;
   auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
@@ -137,54 +131,67 @@ 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;
-
-  stack.addParticle(std::make_tuple(
-      beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
-      plab.normalized(), injectionPos, 0_ns));
-
-  CORSIKA_LOG_INFO("shower axis length: {} ",
-                   (showerCore - injectionPos).getNorm() * 1.02);
-
   ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
                               false, 1000};
-
-  OutputManager output("em_shower_outputs");
-
-  EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
-  // register energy losses as output
-  output.add("dEdX", dEdX);
+  auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
+  uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
+
+  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);
 
   // setup processes, decays and interactions
+  EnergyLossWriter energyloss{showerAxis, dX, nAxisBins};
+  ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
+                                                   energyloss);
 
-  ParticleCut<SubWriter<decltype(dEdX)>> cut(2_MeV, 2_MeV, 100_GeV, 100_GeV, true, dEdX);
   corsika::sibyll::Interaction sibyll{env};
   corsika::sophia::InteractionModel sophia;
-  HEPEnergyType heThresholdNN = 60_GeV;
+  HEPEnergyType heThresholdNN = 80_GeV;
   corsika::proposal::Interaction emCascade(
       env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
-  corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
-  //  BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
+  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");
+
+  output.add("energyloss", energyloss);
+
   TrackWriter tracks;
   output.add("tracks", tracks);
 
-  // long. profile
-  LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
+  LongitudinalWriter profile{showerAxis, nAxisBins, dX};
   output.add("profile", profile);
   LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
 
   Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
-  ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+  ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
       obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
   output.add("particles", observationLevel);
 
+  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;
+
+  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);
 
   output.startOfLibrary();
   Cascade EAS(env, tracking, sequence, output, stack);
@@ -194,7 +201,8 @@ int main(int argc, char** argv) {
 
   EAS.run();
 
-  HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+  HEPEnergyType const Efinal =
+      energyloss.getEnergyLost() + observationLevel.getEnergyGround();
 
   CORSIKA_LOG_INFO(
       "total energy budget (GeV): {}, "
@@ -202,4 +210,6 @@ int main(int argc, char** argv) {
       Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
 
   output.endOfLibrary();
+
+  return EXIT_SUCCESS;
 }

examples/mars.cpp → examples/cascade_examples/mars.cpp

@@ -30,6 +30,7 @@
 
 #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>
 
@@ -55,7 +56,7 @@
 #include <corsika/modules/Sophia.hpp>
 #include <corsika/modules/StackInspector.hpp>
 #include <corsika/modules/TrackWriter.hpp>
-#include <corsika/modules/UrQMD.hpp>
+#include <corsika/modules/FLUKA.hpp>
 
 #include <corsika/setup/SetupStack.hpp>
 #include <corsika/setup/SetupTrajectory.hpp>
@@ -74,8 +75,9 @@ using namespace std;
 
 using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
 using EnvType = Environment<EnvironmentInterface>;
-
-using Particle = setup::Stack<EnvType>::particle_type;
+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;
@@ -111,7 +113,7 @@ void registerRandomStreams(int seed) {
   RNGManager<>::getInstance().registerRandomStream("sibyll");
   RNGManager<>::getInstance().registerRandomStream("sophia");
   RNGManager<>::getInstance().registerRandomStream("pythia");
-  RNGManager<>::getInstance().registerRandomStream("urqmd");
+  RNGManager<>::getInstance().registerRandomStream("fluka");
   RNGManager<>::getInstance().registerRandomStream("proposal");
   if (seed == 0) {
     std::random_device rd;
@@ -200,8 +202,6 @@ int main(int argc, char** argv) {
   }
 
   // check that we got either PDG or A/Z
-  // this can be done with option_groups but the ordering
-  // gets all messed up
   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.");
@@ -224,12 +224,10 @@ int main(int argc, char** argv) {
           Medium::AirDry1Atm,                           // Mars, close enough
           MagneticFieldVector{rootCS, 0_T, 0_uT, 0_T}); // Mars
 
-  builder.setNuclearComposition(                             // Mars
-      {{Code::Nitrogen, Code::Oxygen}, {1. / 3., 2. / 3.}}); // simplified
-  //{{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, Ar removed for now
+  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);
@@ -305,8 +303,10 @@ int main(int argc, char** argv) {
   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
+  uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
 
-  EnergyLossWriter dEdX{showerAxis};
+  EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
   output.add("energyloss", dEdX);
 
   HEPEnergyType const emcut = 1_GeV;
@@ -338,13 +338,13 @@ int main(int argc, char** argv) {
   auto emContinuous =
       make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
 
-  LongitudinalWriter longprof{showerAxis};
+  LongitudinalWriter longprof{showerAxis, nAxisBins, dX};
   output.add("profile", longprof);
   LongitudinalProfile<SubWriter<decltype(longprof)>> profile{longprof};
 
-  corsika::urqmd::UrQMD urqmd;
-  InteractionCounter urqmdCounted{urqmd};
-  StackInspector<setup::Stack<EnvType>> stackInspect(5000, false, E0);
+  corsika::fluka::Interaction leIntModel{env};
+  InteractionCounter leIntCounted{leIntModel};
+  StackInspector<StackType> stackInspect(5000, false, E0);
 
   // assemble all processes into an ordered process list
   struct EnergySwitch {
@@ -354,7 +354,7 @@ int main(int argc, char** argv) {
     bool operator()(Particle const& p) const { return (p.getKineticEnergy() < cutE_); }
   };
   auto hadronSequence =
-      make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyllCounted);
+      make_select(EnergySwitch(heHadronModelThreshold), leIntCounted, sibyllCounted);
 
   // track writer
   TrackWriter trackWriter;
@@ -362,11 +362,14 @@ int main(int argc, char** argv) {
 
   // observation plane
   Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
-  ObservationPlane<setup::Tracking> observationLevel(
-      obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
+  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,
@@ -375,18 +378,19 @@ int main(int argc, char** argv) {
 
   // create the cascade object using the default stack and tracking
   // implementation
-  setup::Tracking tracking;
-  setup::Stack<EnvType> stack;
+  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>());
+    CORSIKA_LOG_INFO("Primary PDG ID:     {}", app["--pdg"]->as<int>());
   } else {
-    CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
+    CORSIKA_LOG_INFO("Primary Z/A:        {}/{}", Z, A);
   }
-  CORSIKA_LOG_INFO("Primary Energy: {}", E0);
-  CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+  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);
 
@@ -409,9 +413,12 @@ int main(int argc, char** argv) {
     stack.clear();
 
     // add the desired particle to the stack
-    stack.addParticle(std::make_tuple(
+    auto const primaryProperties = std::make_tuple(
         beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
-        plab.normalized(), injectionPos, 0_ns));
+        plab.normalized(), injectionPos, 0_ns);
+    stack.addParticle(primaryProperties);
+
+    primaryWriter.recordPrimary(primaryProperties);
 
     // run the shower
     EAS.run();
@@ -424,11 +431,6 @@ int main(int argc, char** argv) {
         "relative difference (%): {}",
         Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
 
-    auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
-
-    save_hist(hists.labHist(), labHist_file, true);
-    save_hist(hists.CMSHist(), cMSHist_file, true);
-
     // trigger the output manager to save this shower to disk
     output.endOfShower();
   }