CONTRIBUTING.md

 # Guidelines for code development, structure, formating etc.
 
-The CORSIKA Project very much welcomes contributions. Here we outlined
+The CORSIKA Project very much welcomes contributions. Here we outline
 how you can find the right place to contribute, and how to do that.
-Connect to https://gitlab.ikp.kit.edu and corsika-devel@lists.kit.edu (self-register at https://www.lists.kit.edu/sympa/subscribe/corsika-devel) to get in touch with the project.
-The CORSIKA Project decides on the [GUIDELINES](CONTRIBUTING.md) and can decide to
+Connect to https://gitlab.iap.kit.edu and corsika-devel@lists.kit.edu (self-register at https://www.lists.kit.edu/sympa/subscribe/corsika-devel) to get in touch with the project.
+The CORSIKA Project decides on the [contributing guidelines](CONTRIBUTING.md) and can decide to
 change/improve them.
 
 # How to contribute
 
-  - We organize all development via `Issues` that may be feature requests,
+  - We organize all development via gitlab `Issues` that may be feature requests,
     ideas, discussions, or bugs fix requests. 
   - New issues can be created, or existing issues
     picked up or contributed to. 
@@ -19,7 +19,8 @@ change/improve them.
     created directly via the gitlab web interface. 
   - Proposed code to close one issue (located in a specific git
     branch) is reviewed, discussed, and eventually merged
-    into the master branch to close the issue.
+    into the master branch via a merge-request (MR) to close the issue.
+  - all merge request will undergo a code review, and must be approved before merge, in order to ensure high code qualtiy: [Code Approval Procedure](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/Code-Approval-Procedure)
 
 
 ## Code formatting
@@ -31,9 +32,9 @@ formatting. We provide a script `do-clang-format.sh`, which can be
 very useful. But we urge everybody to integrate `clang-format` already
 on the level of your source code editor. See [the official
 page](https://clang.llvm.org/docs/ClangFormat.html) for information
-about `clang-format` and its editor integration. At least: run
-`do-clang-format.sh` prior to any `git add` command. Code with
-improper formatting will not be accepted for merging. 
+about `clang-format` and its editor integrations. At least: run
+`do-clang-format.sh` prior to any `git add/commit` command. Code with
+improper formatting will not be accepted for merging. It will trigger automatic warnings by the continuous integration (CI) system.
 
 The definition of source code format is written down in the file
 [.clang-format](.clang-format) and can be changed, if the CORSIKA
@@ -42,87 +43,35 @@ e.g. [link1](https://clangformat.com/) or
 [link2](https://zed0.co.uk/clang-format-configurator/).
 
 
-## Naming conventions
-
-While `clang-format` does the structural formatting, we still need to agree on naming conventions:
-
-  - Classes and structs start with capital letters
-  - Class member variables start with "f"
-  - Any static variable has a "g" prefix. A static member variable starts with "fg"
-  - Class member functions start with capital letters
-  - Any class getter begins with "Get", and setter with "Set". Logical getters start with "Is" or "Has".
-  - enums should be "enum class" and start with a capital "E"
-  - Function parameter names start with "v"
-
-  - We use namespaces to avoid clashes and to structure code
-    - *Everything* is part of the corsika namespace
-    - All classes and objects are encapsulated into suited sub-namespaces,
-      thus corsika::geometry, corsika::processes, corsika::units, etc.
-    - Namespace names do not use capital letters.  
-    - Every header file is copied during build and install into
-      "include/corsika/[namespace]" which also means, each header file
-      can only provide definitions for a _single_ namespace. It is one
-      main purpose of namespaces to structure the location of header
-      files.
-    - Each header file uses an include protection that includes at
-      least the namespace name, and header file name, thus, `#ifndef
-      __include_geometry_Point_h__` or `#ifndef __geometry_Point_h__`,
-      or similar are acceptable.
-    - Header files should always be included with `<..>`, thus,
-      `#include <corsika/geometry/Point.h>` since the build system
-      will always provide the correct include directives (and files
-      anyway cannot be found in file-system paths, which generally do
-      not follow the namespace naming conventions outlined
-      here).
-
-  - Header files are named after the main class (or object) they
-    define. This also means each header file name starts with a
-    capital letter.
-  
-
-## Coding rules
-
-  - Code may not introduce any compiler errors, or warnings
-  - All unit tests must succeed at all times. If tests fail, code is not merged.
-  - We use C++17 concepts wherever useful and helpful
-  - On any major error or malfunction we throw an exception. This is needed and required for complex physics and shower debugging.
-  - We never catch exceptions for error handling, there might be very few special exceptions from this. We need to discuss such cases.
-  - Everything that should not change should be `const`
-  - Everything that does not need to be visible to the outside of a class/struct should be `private` or `protected`
-  - We prefer the use of references, wherever useful
-  - There cannot be any pointer in the interface of any class or object
-    exposed to outside users, there might be pointers for very special cases
-    inside of classes.
-  - When you contribute new code, or extend existing code, at the same time provide unit-tests for all functionality.
-  - When you contribute new physics algorithms, in addition you also need to provide a validation module
-  - Code must be documented with `doxygen` commands
-  
+## Coding rules, conventions and guidelines
+
+Please read the [Coding wiki page](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/Coding-Conventions-and-Guidelines).
+
+
+
 
 ## Release versioning scheme
 
-Releases of CORSIKA are thought to be the baseline for larger scale
-testing, and full production.  The releases are numbered as x.y.z,
-starting with x=8 form the gitlab c++ version. X will only be
-incremented for major design or physics changes. The y index is
-updated for new official releases that normally contain improved or
-enhanced physics performance, and may also contain normal interface
-changes to accomodate improvements. The z index can be updated more
+Releases of CORSIKA 8 are thought to be the baseline for larger scale
+validation, and full production.  The releases are numbered as CORSIKA 8 vx.y.z,
+starting with x=1. The x index is updated for new releases that normally contain improved or
+enhanced physics performance, and also interface
+changes to accomodate improvements. The y and z indices can be updated more
 frequently for bug fixes or new features. Changes in z will not
-contain (major) interface changes, thus, production code will remain
-fully compatible within changes of z.  Special releases of CORSIKA can
-also have a tag name from git, e.g. as in the "milestone1" release.
+contain interface changes, thus, production code will remain
+fully compatible within changes of z.  Special releases of CORSIKA might
+also have a release name.
 
 
 # How to become scientific author of the CORSIKA Project
 
 The CORSIKA Project decides on who becomes scientific author. The
-following conditions are clearly sufficient, but not all of them are
+following conditions are sufficient, but not all of them are
 required all the time:
   - responsibility for a particular functionality or software/management part 
-  - have read and follow these [GUIDELINES](CONTRIBUTING.md)
+  - have read and follow our [contributing guidelines](CONTRIBUTING.md)
   - active in the CORSIKA Project, that means responsive to
-    discussions and problems in corsika-devel@list.kit.edu or on https//gitlab.ikp.kit.edu, 
-    of relevant *Issues*,
-    or in (phone) meetings
-  - agreement to the [COLLABORATION_AGREEMENT](COLLABORATION_AGREEMENT.md) is strictly required
-  - the members of the CORSIKA Project panel agree
+    discussions and problems in corsika-devel@list.kit.edu or on https//gitlab.iap.kit.edu, 
+    of relevant *Issues*, in (phone) meetings or our Mattermost workspace
+  - agreement to the [using and collaborating agreement](USING_COLLABORATING.md) 
+  - the members of the CORSIKA Project must agree

Documentation/CMakeLists.txt

-
-add_subdirectory (Doxygen)
-add_subdirectory (Examples)
-

Documentation/Doxygen/CMakeLists.txt

-
-find_package (Doxygen OPTIONAL_COMPONENTS dot mscgen dia)
-
-if (DOXYGEN_FOUND)
-  if (NOT DOXYGEN_DOT_EXECUTABLE)
-    message (FATAL_ERROR "Found doxygen but not 'dot' command, please install graphviz or set DOXYGEN_DOT_EXECUTABLE")
-  endif()
-
-  set (DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
-  set (DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
-  set (DOXYGEN_GENERATE_HTML YES)
-  set (DOXYGEN_GENERATE_MAN YES)
-  configure_file (${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
-  message ("Start doxygen with \"make doxygen\"")
-  
-  # note the option ALL which allows to build the docs together with the application
-  add_custom_target (doxygen # ALL
-    COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
-    WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-    COMMENT "Generating API documentation with Doxygen"
-    VERBATIM)
-
-  add_custom_command(TARGET doxygen POST_BUILD
-    COMMAND cd ${CMAKE_CURRENT_BINARY_DIR}/latex; pdflatex refman.tex
-    )
-  
-  install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html DESTINATION share/doc OPTIONAL)
-  install(FILES ${CMAKE_CURRENT_BINARY_DIR}/latex/refman.pdf DESTINATION share/doc OPTIONAL)
-         
-else (DOXYGEN_FOUND)
-  
-  message ("Doxygen need to be installed to generate the doxygen documentation")
-  
-endif (DOXYGEN_FOUND)

Documentation/Doxygen/Doxyfile.in

-PROJECT_NAME           = CORSIKA 8
-PROJECT_NUMBER         = 0.0.0
-
-GENERATE_HTML          = YES
-GENERATE_LATEX         = YES
-GENERATE_XML           = YES
-
-OUTPUT_DIRECTORY       = @CMAKE_CURRENT_BINARY_DIR@/
-INPUT                  = @PROJECT_SOURCE_DIR@ @PROJECT_BINARY_DIR@/Framework
-EXCLUDE_PATTERNS       = */ThirdParty/*/*
-FULL_PATH_NAMES        = YES
-STRIP_FROM_PATH        = @PROJECT_SOURCE_DIR@
-
-FILE_PATTERNS          = *.cc *.cpp *.cxx *.h *.dox *.inc *.md
-EXTENSION_MAPPING      = inc=C++
-RECURSIVE              = YES
-
-SOURCE_BROWSER         = YES
-# INLINE_SOURCES
-GENERATE_TREEVIEW      = YES
-
-CLASS_DIAGRAMS         = NO
-HAVE_DOT               = YES
-CLASS_GRAPH            = YES
-COLLABORATION_GRAPH    = NO
-UML_LOOK               = NO
-TEMPLATE_RELATIONS     = YES
-INCLUDE_GRAPH          = YES
-GRAPHICAL_HIERARCHY    = YES
-DIRECTORY_GRAPH        = YES
-GENERATE_LEGEND        = YES
-USE_MATHJAX            = YES
-
-SEARCHENGINE           = YES

Documentation/Examples/CMakeLists.txt

-add_executable (helix_example helix_example.cc)
-target_compile_options(helix_example PRIVATE -g) # do not skip asserts
-target_link_libraries(helix_example CORSIKAgeometry CORSIKAunits)
-install (TARGETS helix_example DESTINATION share/examples)
-CORSIKA_ADD_TEST (helix_example)
-
-add_executable (geometry_example geometry_example.cc)
-target_compile_options(geometry_example PRIVATE -g) # do not skip asserts
-target_link_libraries (geometry_example CORSIKAgeometry CORSIKAunits)
-install (TARGETS geometry_example DESTINATION share/examples)
-CORSIKA_ADD_TEST (geometry_example)
-
-add_executable (logger_example logger_example.cc)
-target_compile_options(logger_example PRIVATE -g) # do not skip asserts
-target_link_libraries (logger_example CORSIKAunits CORSIKAlogging)
-install (TARGETS logger_example DESTINATION share/examples)
-CORSIKA_ADD_TEST (logger_example)
-
-add_executable (stack_example stack_example.cc)
-target_compile_options(stack_example PRIVATE -g) # do not skip asserts
-target_link_libraries (stack_example SuperStupidStack CORSIKAunits
-  CORSIKAlogging)
-CORSIKA_ADD_TEST (stack_example)
-
-if (Pythia8_FOUND)
-  add_executable (cascade_example cascade_example.cc)
-  target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
-  target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
-    CORSIKAlogging
-    CORSIKArandom
-    ProcessSibyll
-    ProcessPythia
-    CORSIKAcascade
-    ProcessEnergyLoss
-    ProcessStackInspector
-    ProcessTrackWriter
-    ProcessTrackingLine
-    CORSIKAprocesses
-    CORSIKAcascade
-    CORSIKAparticles
-    CORSIKAgeometry
-    CORSIKAenvironment
-    CORSIKAprocesssequence
-    )
-  install (TARGETS cascade_example DESTINATION share/examples)
-  CORSIKA_ADD_TEST (cascade_example)
-endif()
-
-add_executable (boundary_example boundary_example.cc)
-target_compile_options(boundary_example PRIVATE -g) # do not skip asserts
-target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlogging
-  CORSIKArandom
-  ProcessSibyll
-  CORSIKAcascade
-  ProcessTrackWriter
-  ProcessTrackingLine
-  CORSIKAprocesses
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  )
-
-install (TARGETS boundary_example DESTINATION share/examples)
-CORSIKA_ADD_TEST (boundary_example)
-
-if (Pythia8_FOUND)
-  add_executable (cascade_proton_example cascade_proton_example.cc)
-  target_compile_options(cascade_proton_example PRIVATE -g) # do not skip asserts
-  target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
-    CORSIKAlogging
-    CORSIKArandom
-    ProcessSibyll
-    ProcessPythia
-    CORSIKAcascade
-    ProcessEnergyLoss
-    ProcessTrackWriter
-    ProcessTrackingLine
-    ProcessHadronicElasticModel
-    CORSIKAprocesses
-    CORSIKAcascade
-    CORSIKAparticles
-    CORSIKAgeometry
-    CORSIKAenvironment
-    CORSIKAprocesssequence
-    )
-
-  install (TARGETS cascade_proton_example DESTINATION share/examples)
-  CORSIKA_ADD_TEST (cascade_proton_example)
-endif()
-
-if (Pythia8_FOUND)
-  add_executable (vertical_EAS vertical_EAS.cc)
-  target_compile_options(vertical_EAS PRIVATE -g) # do not skip asserts
-  target_link_libraries (vertical_EAS SuperStupidStack CORSIKAunits
-    CORSIKAlogging
-    CORSIKArandom
-    ProcessSibyll
-    ProcessPythia
-    CORSIKAcascade
-    ProcessEnergyLoss
-    ProcessTrackWriter
-    ProcessTrackingLine
-    ProcessHadronicElasticModel
-    CORSIKAprocesses
-    CORSIKAcascade
-    CORSIKAparticles
-    CORSIKAgeometry
-    CORSIKAenvironment
-    CORSIKAprocesssequence
-    )
-  install (TARGETS vertical_EAS DESTINATION share/examples)
-  # CORSIKA_ADD_TEST (vertical_EAS) not yet...
-endif()
-
-add_executable (staticsequence_example staticsequence_example.cc)
-target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
-target_link_libraries (staticsequence_example
-  CORSIKAprocesssequence
-  CORSIKAunits
-  CORSIKAgeometry
-  CORSIKAlogging)
-install (TARGETS staticsequence_example DESTINATION share/examples)
-CORSIKA_ADD_TEST (staticsequence_example)

Documentation/Examples/boundary_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/cascade/Cascade.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/setup/SetupEnvironment.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-#include <corsika/geometry/Sphere.h>
-
-#include <corsika/process/sibyll/Decay.h>
-#include <corsika/process/sibyll/Interaction.h>
-#include <corsika/process/sibyll/NuclearInteraction.h>
-
-#include <corsika/process/track_writer/TrackWriter.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <iostream>
-#include <limits>
-#include <typeinfo>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::particles;
-using namespace corsika::random;
-using namespace corsika::setup;
-using namespace corsika::geometry;
-using namespace corsika::environment;
-
-using namespace std;
-using namespace corsika::units::si;
-
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
-  HEPEnergyType fECut;
-
-  HEPEnergyType fEnergy = 0_GeV;
-  HEPEnergyType fEmEnergy = 0_GeV;
-  int fEmCount = 0;
-  HEPEnergyType fInvEnergy = 0_GeV;
-  int fInvCount = 0;
-
-public:
-  ProcessCut(const HEPEnergyType cut)
-      : fECut(cut) {}
-
-  template <typename Particle>
-  bool isBelowEnergyCut(Particle& p) const {
-    // nuclei
-    if (p.GetPID() == particles::Code::Nucleus) {
-      auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
-      auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
-      if (ElabNuc < fECut || EcmNN < 10_GeV)
-        return true;
-      else
-        return false;
-    } else {
-      // TODO: center-of-mass energy hard coded
-      const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
-      if (p.GetEnergy() < fECut || Ecm < 10_GeV)
-        return true;
-      else
-        return false;
-    }
-  }
-
-  bool isEmParticle(Code pCode) const {
-    bool is_em = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::Electron:
-        is_em = true;
-        break;
-      case Code::Positron:
-        is_em = true;
-        break;
-      case Code::Gamma:
-        is_em = true;
-        break;
-      default:
-        break;
-    }
-    return is_em;
-  }
-
-  void defineEmParticles() const {
-    // create bool array identifying em particles
-  }
-
-  bool isInvisible(Code pCode) const {
-    bool is_inv = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::NuE:
-        is_inv = true;
-        break;
-      case Code::NuEBar:
-        is_inv = true;
-        break;
-      case Code::NuMu:
-        is_inv = true;
-        break;
-      case Code::NuMuBar:
-        is_inv = true;
-        break;
-      case Code::MuPlus:
-        is_inv = true;
-        break;
-      case Code::MuMinus:
-        is_inv = true;
-        break;
-
-      case Code::Neutron:
-        is_inv = true;
-        break;
-
-      case Code::AntiNeutron:
-        is_inv = true;
-        break;
-
-      default:
-        break;
-    }
-    return is_inv;
-  }
-
-  template <typename TSecondaries>
-  EProcessReturn DoSecondaries(TSecondaries& vS) {
-    auto p = vS.begin();
-    while (p != vS.end()) {
-      const Code pid = p.GetPID();
-      HEPEnergyType energy = p.GetEnergy();
-      cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
-           << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
-           << endl;
-      if (isEmParticle(pid)) {
-        cout << "removing em. particle..." << endl;
-        fEmEnergy += energy;
-        fEmCount += 1;
-        p.Delete();
-      } else if (isInvisible(pid)) {
-        cout << "removing inv. particle..." << endl;
-        fInvEnergy += energy;
-        fInvCount += 1;
-        p.Delete();
-      } else if (isBelowEnergyCut(p)) {
-        cout << "removing low en. particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else if (p.GetTime() > 10_ms) {
-        cout << "removing OLD particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else {
-        ++p; // next entry in SecondaryView
-      }
-    }
-    return EProcessReturn::eOk;
-  }
-
-  void Init() {
-    fEmEnergy = 0. * 1_GeV;
-    fEmCount = 0;
-    fInvEnergy = 0. * 1_GeV;
-    fInvCount = 0;
-    fEnergy = 0. * 1_GeV;
-    // defineEmParticles();
-  }
-
-  void ShowResults() {
-    cout << " ******************************" << endl
-         << " ParticleCut: " << endl
-         << " energy in em.  component (GeV):  " << fEmEnergy / 1_GeV << endl
-         << " no. of em.  particles injected:  " << fEmCount << endl
-         << " energy in inv. component (GeV):  " << fInvEnergy / 1_GeV << endl
-         << " no. of inv. particles injected:  " << fInvCount << endl
-         << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
-         << " ******************************" << endl;
-  }
-
-  HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
-  HEPEnergyType GetCutEnergy() const { return fEnergy; }
-  HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
-
-template <bool deleteParticle>
-struct MyBoundaryCrossingProcess
-    : public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
-
-  MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
-
-  template <typename Particle>
-  EProcessReturn DoBoundaryCrossing(Particle& p,
-                                    typename Particle::BaseNodeType const& from,
-                                    typename Particle::BaseNodeType const& to) {
-    std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
-
-    auto const& name = particles::GetName(p.GetPID());
-    auto const start = p.GetPosition().GetCoordinates();
-
-    fFile << name << "    " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
-          << start[2] / 1_m << '\n';
-
-    if constexpr (deleteParticle) { p.Delete(); }
-
-    return EProcessReturn::eOk;
-  }
-
-  void Init() {}
-
-private:
-  std::ofstream fFile;
-};
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
-
-  // setup environment, geometry
-  using EnvType = Environment<setup::IEnvironmentModel>;
-  EnvType env;
-  auto& universe = *(env.GetUniverse());
-
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  auto outerMedium = EnvType::CreateNode<Sphere>(
-      Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
-  // fraction of oxygen
-  auto const props =
-      outerMedium
-          ->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
-              1_kg / (1_m * 1_m * 1_m),
-              environment::NuclearComposition(
-                  std::vector<particles::Code>{particles::Code::Proton},
-                  std::vector<float>{1.f}));
-
-  auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
-
-  innerMedium->SetModelProperties(props);
-
-  outerMedium->AddChild(std::move(innerMedium));
-
-  universe.AddChild(std::move(outerMedium));
-
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-
-  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-  process::sibyll::Interaction sibyll;
-  process::sibyll::Decay decay{{particles::Code::PiPlus, particles::Code::PiMinus,
-                                particles::Code::KPlus, particles::Code::KMinus,
-                                particles::Code::K0Long, particles::Code::K0Short}};
-  ProcessCut cut(20_GeV);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-  MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
-
-  // assemble all processes into an ordered process list
-  auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter;
-
-  // setup particle stack, and add primary particles
-  setup::Stack stack;
-  stack.Clear();
-  const Code beamCode = Code::Proton;
-  const HEPMassType mass = particles::GetMass(Code::Proton);
-  const HEPEnergyType E0 = 50_TeV;
-
-  std::uniform_real_distribution distTheta(0., 180.);
-  std::uniform_real_distribution distPhi(0., 360.);
-  std::mt19937 rng;
-
-  for (int i = 0; i < 100; ++i) {
-    auto const theta = distTheta(rng);
-    auto const phi = distPhi(rng);
-
-    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 = corsika::stack::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;
-    Point pos(rootCS, 0_m, 0_m, 0_m);
-    stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            beamCode, E0, plab, pos, 0_ns});
-  }
-
-  // define air shower object, run simulation
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-  EAS.Init();
-  EAS.Run();
-
-  cout << "Result: E0=" << E0 / 1_GeV << endl;
-  cut.ShowResults();
-  const HEPEnergyType Efinal =
-      cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
-  cout << "total energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
-}

Documentation/Examples/cascade_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/cascade/Cascade.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/energy_loss/EnergyLoss.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/setup/SetupEnvironment.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-#include <corsika/geometry/Sphere.h>
-
-#include <corsika/process/sibyll/Decay.h>
-#include <corsika/process/sibyll/Interaction.h>
-#include <corsika/process/sibyll/NuclearInteraction.h>
-
-#include <corsika/process/pythia/Decay.h>
-
-#include <corsika/process/track_writer/TrackWriter.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <iostream>
-#include <limits>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::particles;
-using namespace corsika::random;
-using namespace corsika::setup;
-using namespace corsika::geometry;
-using namespace corsika::environment;
-
-using namespace std;
-using namespace corsika::units::si;
-
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
-  HEPEnergyType fECut;
-
-  HEPEnergyType fEnergy = 0_GeV;
-  HEPEnergyType fEmEnergy = 0_GeV;
-  int fEmCount = 0;
-  HEPEnergyType fInvEnergy = 0_GeV;
-  int fInvCount = 0;
-
-public:
-  ProcessCut(const HEPEnergyType cut)
-      : fECut(cut) {}
-
-  template <typename Particle>
-  bool isBelowEnergyCut(Particle& p) const {
-    // nuclei
-    if (p.GetPID() == particles::Code::Nucleus) {
-      auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
-      auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
-      if (ElabNuc < fECut || EcmNN < 10_GeV)
-        return true;
-      else
-        return false;
-    } else {
-      // TODO: center-of-mass energy hard coded
-      const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
-      if (p.GetEnergy() < fECut || Ecm < 10_GeV)
-        return true;
-      else
-        return false;
-    }
-  }
-
-  bool isEmParticle(Code pCode) const {
-    bool is_em = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::Electron:
-        is_em = true;
-        break;
-      case Code::Positron:
-        is_em = true;
-        break;
-      case Code::Gamma:
-        is_em = true;
-        break;
-      default:
-        break;
-    }
-    return is_em;
-  }
-
-  void defineEmParticles() const {
-    // create bool array identifying em particles
-  }
-
-  bool isInvisible(Code pCode) const {
-    bool is_inv = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::NuE:
-        is_inv = true;
-        break;
-      case Code::NuEBar:
-        is_inv = true;
-        break;
-      case Code::NuMu:
-        is_inv = true;
-        break;
-      case Code::NuMuBar:
-        is_inv = true;
-        break;
-      case Code::MuPlus:
-        is_inv = true;
-        break;
-      case Code::MuMinus:
-        is_inv = true;
-        break;
-
-      case Code::Neutron:
-        is_inv = true;
-        break;
-
-      case Code::AntiNeutron:
-        is_inv = true;
-        break;
-
-      default:
-        break;
-    }
-    return is_inv;
-  }
-
-  template <typename TSecondaries>
-  EProcessReturn DoSecondaries(TSecondaries& vS) {
-    auto p = vS.begin();
-    while (p != vS.end()) {
-      const Code pid = p.GetPID();
-      HEPEnergyType energy = p.GetEnergy();
-      cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
-           << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
-           << endl;
-      if (isEmParticle(pid)) {
-        cout << "removing em. particle..." << endl;
-        fEmEnergy += energy;
-        fEmCount += 1;
-        p.Delete();
-      } else if (isInvisible(pid)) {
-        cout << "removing inv. particle..." << endl;
-        fInvEnergy += energy;
-        fInvCount += 1;
-        p.Delete();
-      } else if (isBelowEnergyCut(p)) {
-        cout << "removing low en. particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else if (p.GetTime() > 10_ms) {
-        cout << "removing OLD particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else {
-        ++p; // next entry in SecondaryView
-      }
-    }
-    return EProcessReturn::eOk;
-  }
-
-  void Init() {
-    fEmEnergy = 0. * 1_GeV;
-    fEmCount = 0;
-    fInvEnergy = 0. * 1_GeV;
-    fInvCount = 0;
-    fEnergy = 0. * 1_GeV;
-    // defineEmParticles();
-  }
-
-  void ShowResults() {
-    cout << " ******************************" << endl
-         << " ParticleCut: " << endl
-         << " energy in em.  component (GeV):  " << fEmEnergy / 1_GeV << endl
-         << " no. of em.  particles injected:  " << fEmCount << endl
-         << " energy in inv. component (GeV):  " << fInvEnergy / 1_GeV << endl
-         << " no. of inv. particles injected:  " << fInvCount << endl
-         << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
-         << " ******************************" << endl;
-  }
-
-  HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
-  HEPEnergyType GetCutEnergy() const { return fEnergy; }
-  HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-
-  const LengthType height_atmosphere = 112.8_km;
-
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
-
-  // setup environment, geometry
-  using EnvType = environment::Environment<setup::IEnvironmentModel>;
-  EnvType env;
-  auto& universe = *(env.GetUniverse());
-
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  auto outerMedium = EnvType::CreateNode<Sphere>(
-      Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
-  // fraction of oxygen
-  const float fox = 0.20946;
-  auto const props =
-      outerMedium
-          ->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
-              1_kg / (1_m * 1_m * 1_m),
-              environment::NuclearComposition(
-                  std::vector<particles::Code>{particles::Code::Nitrogen,
-                                               particles::Code::Oxygen},
-                  std::vector<float>{1.f - fox, fox}));
-
-  auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
-
-  innerMedium->SetModelProperties(props);
-
-  outerMedium->AddChild(std::move(innerMedium));
-
-  universe.AddChild(std::move(outerMedium));
-
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-  stack_inspector::StackInspector<setup::Stack> stackInspect(1, true);
-
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-  process::sibyll::Interaction sibyll;
-  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay(trackedHadrons);
-  ProcessCut cut(20_GeV);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-  process::energy_loss::EnergyLoss eLoss;
-
-  // assemble all processes into an ordered process list
-  auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
-                               << trackWriter;
-
-  // 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 = GetNucleusMass(nuclA, nuclZ);
-  const HEPEnergyType E0 = nuclA * 1_TeV;
-  double theta = 0.;
-  double phi = 0.;
-
-  {
-    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 = corsika::stack::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;
-    Point pos(rootCS, 0_m, 0_m,
-              height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
-    stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                 corsika::stack::MomentumVector, geometry::Point,
-                                 units::si::TimeType, unsigned short, unsigned short>{
-        beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
-  }
-
-  // define air shower object, run simulation
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-  EAS.Init();
-  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;
-}

Documentation/Examples/cascade_proton_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/cascade/Cascade.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-#include <corsika/geometry/Sphere.h>
-
-#include <corsika/process/sibyll/Decay.h>
-#include <corsika/process/sibyll/Interaction.h>
-#include <corsika/process/sibyll/NuclearInteraction.h>
-
-#include <corsika/process/pythia/Decay.h>
-#include <corsika/process/pythia/Interaction.h>
-
-#include <corsika/process/track_writer/TrackWriter.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <boost/type_index.hpp>
-using boost::typeindex::type_id_with_cvr;
-
-#include <iostream>
-#include <limits>
-#include <typeinfo>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::particles;
-using namespace corsika::random;
-using namespace corsika::setup;
-using namespace corsika::geometry;
-using namespace corsika::environment;
-
-using namespace std;
-using namespace corsika::units::si;
-
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
-  HEPEnergyType fECut;
-
-  HEPEnergyType fEnergy = 0_GeV;
-  HEPEnergyType fEmEnergy = 0_GeV;
-  int fEmCount = 0;
-  HEPEnergyType fInvEnergy = 0_GeV;
-  int fInvCount = 0;
-
-public:
-  ProcessCut(const HEPEnergyType cut)
-      : fECut(cut) {}
-
-  template <typename Particle>
-  bool isBelowEnergyCut(Particle& p) const {
-    // nuclei
-    if (p.GetPID() == particles::Code::Nucleus) {
-      auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
-      auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
-      if (ElabNuc < fECut || EcmNN < 10_GeV)
-        return true;
-      else
-        return false;
-    } else {
-      // TODO: center-of-mass energy hard coded
-      const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
-      if (p.GetEnergy() < fECut || Ecm < 10_GeV)
-        return true;
-      else
-        return false;
-    }
-  }
-
-  bool isEmParticle(Code pCode) const {
-    bool is_em = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::Electron:
-        is_em = true;
-        break;
-      case Code::Positron:
-        is_em = true;
-        break;
-      case Code::Gamma:
-        is_em = true;
-        break;
-      default:
-        break;
-    }
-    return is_em;
-  }
-
-  void defineEmParticles() const {
-    // create bool array identifying em particles
-  }
-
-  bool isInvisible(Code pCode) const {
-    bool is_inv = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::NuE:
-        is_inv = true;
-        break;
-      case Code::NuEBar:
-        is_inv = true;
-        break;
-      case Code::NuMu:
-        is_inv = true;
-        break;
-      case Code::NuMuBar:
-        is_inv = true;
-        break;
-      case Code::MuPlus:
-        is_inv = true;
-        break;
-      case Code::MuMinus:
-        is_inv = true;
-        break;
-
-      case Code::Neutron:
-        is_inv = true;
-        break;
-
-      case Code::AntiNeutron:
-        is_inv = true;
-        break;
-
-      default:
-        break;
-    }
-    return is_inv;
-  }
-
-  template <typename TSecondaries>
-  EProcessReturn DoSecondaries(TSecondaries& vS) {
-    auto p = vS.begin();
-    while (p != vS.end()) {
-      const Code pid = p.GetPID();
-      HEPEnergyType energy = p.GetEnergy();
-      cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
-           << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
-           << endl;
-      if (isEmParticle(pid)) {
-        cout << "removing em. particle..." << endl;
-        fEmEnergy += energy;
-        fEmCount += 1;
-        p.Delete();
-      } else if (isInvisible(pid)) {
-        cout << "removing inv. particle..." << endl;
-        fInvEnergy += energy;
-        fInvCount += 1;
-        p.Delete();
-      } else if (isBelowEnergyCut(p)) {
-        cout << "removing low en. particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else if (p.GetTime() > 10_ms) {
-        cout << "removing OLD particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else {
-        ++p; // next entry in SecondaryView
-      }
-    }
-    return EProcessReturn::eOk;
-  }
-
-  void Init() {
-    fEmEnergy = 0. * 1_GeV;
-    fEmCount = 0;
-    fInvEnergy = 0. * 1_GeV;
-    fInvCount = 0;
-    fEnergy = 0. * 1_GeV;
-    // defineEmParticles();
-  }
-
-  void ShowResults() {
-    cout << " ******************************" << endl
-         << " ParticleCut: " << endl
-         << " energy in em.  component (GeV):  " << fEmEnergy / 1_GeV << endl
-         << " no. of em.  particles injected:  " << fEmCount << endl
-         << " energy in inv. component (GeV):  " << fInvEnergy / 1_GeV << endl
-         << " no. of inv. particles injected:  " << fInvCount << endl
-         << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
-         << " ******************************" << endl;
-  }
-
-  HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
-  HEPEnergyType GetCutEnergy() const { return fEnergy; }
-  HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
-
-  // setup environment, geometry
-  using EnvType = Environment<setup::IEnvironmentModel>;
-  EnvType env;
-  auto& universe = *(env.GetUniverse());
-
-  auto theMedium =
-      EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
-                                  1_km * std::numeric_limits<double>::infinity());
-
-  using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
-  theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_kg / (1_m * 1_m * 1_m),
-      NuclearComposition(std::vector<particles::Code>{particles::Code::Hydrogen},
-                         std::vector<float>{(float)1.}));
-
-  universe.AddChild(std::move(theMedium));
-
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-  //  process::sibyll::Interaction sibyll(env);
-  process::pythia::Interaction pythia;
-  //  process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
-  //  process::sibyll::Decay decay(trackedHadrons);
-  process::pythia::Decay decay(trackedHadrons);
-  ProcessCut cut(20_GeV);
-
-  // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
-  // process::HadronicElasticModel::HadronicElasticInteraction
-  // hadronicElastic(env);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-
-  // assemble all processes into an ordered process list
-  // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
-  auto sequence = pythia << decay << cut << trackWriter;
-
-  // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
-  // << "\n";
-
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.Clear();
-  const Code beamCode = Code::Proton;
-  const HEPMassType mass = particles::Proton::GetMass();
-  const HEPEnergyType E0 = 100_GeV;
-  double theta = 0.;
-  double phi = 0.;
-
-  {
-    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-      return sqrt(Elab * Elab - m * 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 = corsika::stack::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;
-    Point pos(rootCS, 0_m, 0_m, 0_m);
-    stack.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            beamCode, E0, plab, pos, 0_ns});
-  }
-
-  // define air shower object, run simulation
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-  EAS.Init();
-  EAS.Run();
-
-  cout << "Result: E0=" << E0 / 1_GeV << endl;
-  cut.ShowResults();
-  const HEPEnergyType Efinal =
-      cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
-  cout << "total energy (GeV): " << Efinal / 1_GeV << endl
-       << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
-}

Documentation/Examples/geometry_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cstdlib>
-#include <iostream>
-#include <typeinfo>
-
-using namespace corsika;
-using namespace corsika::geometry;
-using namespace corsika::units::si;
-
-int main() {
-  // define the root coordinate system
-  geometry::CoordinateSystem& root =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  // another CS defined by a translation relative to the root CS
-  CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});
-
-  // rotations are possible, too; parameters are axis vector and angle
-  CoordinateSystem cs3 =
-      root.rotate(QuantityVector<length_d>{1_m, 0_m, 0_m}, 90 * degree_angle);
-
-  // now let's define some geometrical objects:
-  Point const p1(root, {0_m, 0_m, 0_m}); // the origin of the root CS
-  Point const p2(cs2, {0_m, 0_m, 0_m});  // the origin of cs2
-
-  Vector<length_d> const diff =
-      p2 -
-      p1; // the distance between the points, basically the translation vector given above
-  auto const norm = diff.squaredNorm(); // squared length with the right dimension
-
-  // print the components of the vector as given in the different CS
-  std::cout << "p2-p1 components in root: " << diff.GetComponents(root) << std::endl;
-  std::cout << "p2-p1 components in cs2: " << diff.GetComponents(cs2)
-            << std::endl; // by definition invariant under translations
-  std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3)
-            << std::endl; // but not under rotations
-  std::cout << "p2-p1 norm^2: " << norm << std::endl;
-  assert(norm == 1 * meter * meter);
-
-  Sphere s(p1, 10_m); // define a sphere around a point with a radius
-  std::cout << "p1 inside s:  " << s.Contains(p2) << std::endl;
-  assert(s.Contains(p2) == 1);
-
-  Sphere s2(p1, 3_um); // another sphere
-  std::cout << "p1 inside s2: " << s2.Contains(p2) << std::endl;
-  assert(s2.Contains(p2) == 0);
-
-  // let's try parallel projections:
-  auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m});
-  auto const v2 = Vector<length_d>(root, {1_m, 0_m, 0_m});
-
-  auto const v3 = v1.parallelProjectionOnto(v2);
-
-  // cross product
-  auto const cross =
-      v1.cross(v2).normalized(); // normalized() returns dimensionless, normalized vectors
-
-  // if a CS is not given as parameter for getComponents(), the components
-  // in the "home" CS are returned
-  std::cout << "v1: " << v1.GetComponents() << std::endl;
-  std::cout << "v2: " << v2.GetComponents() << std::endl;
-  std::cout << "parallel projection of v1 onto v2: " << v3.GetComponents() << std::endl;
-  std::cout << "normalized cross product of v1 x v2" << cross.GetComponents()
-            << std::endl;
-
-  return EXIT_SUCCESS;
-}

Documentation/Examples/helix_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/geometry/Helix.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <array>
-#include <cstdlib>
-#include <iostream>
-
-using namespace corsika;
-using namespace corsika::geometry;
-using namespace corsika::units::si;
-
-int main() {
-  geometry::CoordinateSystem& root =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  Point const r0(root, {0_m, 0_m, 0_m});
-  auto const omegaC = 2 * M_PI * 1_Hz;
-  Vector<speed_d> vPar(root, {0_m / second, 0_m / second, 10_cm / second});
-  Vector<speed_d> vPerp(root, {1_m / second, 0_m / second, 0_m / second});
-
-  Helix h(r0, omegaC, vPar, vPerp);
-
-  auto constexpr t0 = 0_s;
-  auto constexpr t1 = 1_s;
-  auto constexpr dt = 1_ms;
-  auto constexpr n = long((t1 - t0) / dt) + 1;
-
-  auto arr = std::make_unique<std::array<std::array<double, 4>, n>>();
-  auto& positions = *arr;
-
-  for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i) {
-    auto const r = h.GetPosition(t).GetCoordinates();
-
-    positions[i][0] = t / 1_s;
-    positions[i][1] = r[0] / 1_m;
-    positions[i][2] = r[1] / 1_m;
-    positions[i][3] = r[2] / 1_m;
-  }
-
-  std::cout << positions[n - 2][0] << " " << positions[n - 2][1] << " "
-            << positions[n - 2][2] << " " << positions[n - 2][3] << std::endl;
-
-  return EXIT_SUCCESS;
-}

Documentation/Examples/logger_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/logging/Logger.h>
-#include <boost/format.hpp>
-#include <fstream>
-#include <iostream>
-#include <string>
-
-using namespace std;
-using namespace corsika::logging;
-
-int main() {
-  {
-    cout << "writing to \"another.log\"" << endl;
-    ofstream logfile("another.log");
-    sink::SinkStream unbuffered_sink(logfile);
-    sink::BufferedSinkStream sink(logfile, sink::StdBuffer(10000));
-    Logger<MessageOn, sink::BufferedSinkStream> info("\033[32m", "info", sink);
-    Logger<MessageOn, sink::BufferedSinkStream> err("\033[31m", "error", sink);
-    // logger<ostream,messageconst,StdBuffer> info(std::cout, StdBuffer(10000));
-
-    /*
-      Logging& logs = Logging::GetInstance();
-      logs.AddLogger<>("info", info);
-      auto& log_1 = logs.GetLogger("info"); // no so useful, since type of log_1 is
-      std::any
-    */
-
-    for (int i = 0; i < 10000; ++i) {
-      LOG(info, "irgendwas", " ", string("and more"), " ",
-          boost::format("error: %i message: %s. done."), i, "stupido");
-      LOG(err, "Fehler");
-    }
-  }
-
-  {
-    sink::NoSink off;
-    Logger<MessageOff> info("", "", off);
-
-    for (int i = 0; i < 10000; ++i) {
-      LOG(info, "irgendwas", string("and more"),
-          boost::format("error: %i message: %s. done."), i, "stupido", "a-number:", 8.99,
-          "ENDE");
-    }
-  }
-
-  return 0;
-}

Documentation/Examples/stack_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/particles/ParticleProperties.h>
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-
-#include <cassert>
-#include <iomanip>
-#include <iostream>
-
-using namespace corsika;
-using namespace corsika::units::si;
-using namespace corsika::stack;
-using namespace corsika::geometry;
-using namespace std;
-
-void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
-  const geometry::CoordinateSystem& rootCS =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-  for (int i = 0; i < 11; ++i) {
-    s.AddParticle(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            particles::Code::Electron, 1.5_GeV * i,
-            corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
-            geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-  }
-}
-
-void read(corsika::stack::super_stupid::SuperStupidStack& s) {
-  assert(s.GetSize() == 11); // stack has 11 particles
-
-  HEPEnergyType total_energy;
-  int i = 0;
-  for (auto& p : s) {
-    total_energy += p.GetEnergy();
-    // particles are electrons with 1.5 GeV energy times i
-    assert(p.GetPID() == particles::Code::Electron);
-    assert(p.GetEnergy() == 1.5_GeV * (i++));
-  }
-}
-
-int main() {
-  corsika::stack::super_stupid::SuperStupidStack s;
-  fill(s);
-  read(s);
-  return 0;
-}

Documentation/Examples/staticsequence_example.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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 <array>
-#include <iomanip>
-#include <iostream>
-
-#include <corsika/process/ProcessSequence.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-using namespace corsika;
-using namespace corsika::units::si;
-using namespace corsika::process;
-using namespace std;
-
-const int nData = 10;
-
-class Process1 : public BaseProcess<Process1> {
-public:
-  Process1() {}
-  template <typename D, typename T, typename S>
-  EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] += 1;
-    return EProcessReturn::eOk;
-  }
-};
-
-class Process2 : public BaseProcess<Process2> {
-public:
-  Process2() {}
-
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
-    return EProcessReturn::eOk;
-  }
-};
-
-class Process3 : public BaseProcess<Process3> {
-public:
-  Process3() {}
-
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D&, T&, S&) const {
-    return EProcessReturn::eOk;
-  }
-};
-
-class Process4 : public BaseProcess<Process4> {
-public:
-  Process4(const double v)
-      : fV(v) {}
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
-    for (int i = 0; i < nData; ++i) d.p[i] *= fV;
-    return EProcessReturn::eOk;
-  }
-
-private:
-  double fV;
-};
-
-struct DummyData {
-  double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-};
-struct DummyStack {
-  void clear() {}
-};
-struct DummyTrajectory {};
-
-void modular() {
-
-  Process1 m1;
-  Process2 m2;
-  Process3 m3;
-  Process4 m4(0.9);
-
-  auto sequence = m1 << m2 << m3 << m4;
-
-  DummyData particle;
-  DummyTrajectory track;
-
-  const int n = 1000;
-  for (int i = 0; i < n; ++i) { sequence.DoContinuous(particle, track); }
-
-  for (int i = 0; i < nData; ++i) {
-    // cout << p.p[i] << endl;
-    // assert(p.p[i] == n-i*100);
-  }
-
-  cout << " done (nothing...) " << endl;
-}
-
-int main() {
-  modular();
-  return 0;
-}

Documentation/Examples/vertical_EAS.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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/cascade/Cascade.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/energy_loss/EnergyLoss.h>
-#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-#include <corsika/geometry/Sphere.h>
-
-#include <corsika/process/sibyll/Decay.h>
-#include <corsika/process/sibyll/Interaction.h>
-#include <corsika/process/sibyll/NuclearInteraction.h>
-
-#include <corsika/process/pythia/Decay.h>
-
-#include <corsika/process/track_writer/TrackWriter.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/utl/CorsikaFenv.h>
-
-#include <boost/type_index.hpp>
-using boost::typeindex::type_id_with_cvr;
-
-#include <iostream>
-#include <limits>
-#include <typeinfo>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::particles;
-using namespace corsika::random;
-using namespace corsika::setup;
-using namespace corsika::geometry;
-using namespace corsika::environment;
-
-using namespace std;
-using namespace corsika::units::si;
-
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
-  HEPEnergyType fECut;
-
-  HEPEnergyType fEnergy = 0_GeV;
-  HEPEnergyType fEmEnergy = 0_GeV;
-  int fEmCount = 0;
-  HEPEnergyType fInvEnergy = 0_GeV;
-  int fInvCount = 0;
-
-public:
-  ProcessCut(const HEPEnergyType cut)
-      : fECut(cut) {}
-
-  template <typename Particle>
-  bool isBelowEnergyCut(Particle& p) const {
-    // nuclei
-    if (p.GetPID() == particles::Code::Nucleus) {
-      auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
-      auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
-      if (ElabNuc < fECut || EcmNN < 10_GeV)
-        return true;
-      else
-        return false;
-    } else {
-      // TODO: center-of-mass energy hard coded
-      const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
-      if (p.GetEnergy() < fECut || Ecm < 10_GeV)
-        return true;
-      else
-        return false;
-    }
-  }
-
-  bool isEmParticle(Code pCode) const {
-    bool is_em = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::Electron:
-        is_em = true;
-        break;
-      case Code::Positron:
-        is_em = true;
-        break;
-      case Code::Gamma:
-        is_em = true;
-        break;
-      default:
-        break;
-    }
-    return is_em;
-  }
-
-  void defineEmParticles() const {
-    // create bool array identifying em particles
-  }
-
-  bool isInvisible(Code pCode) const {
-    bool is_inv = false;
-    // FOR NOW: switch
-    switch (pCode) {
-      case Code::NuE:
-        is_inv = true;
-        break;
-      case Code::NuEBar:
-        is_inv = true;
-        break;
-      case Code::NuMu:
-        is_inv = true;
-        break;
-      case Code::NuMuBar:
-        is_inv = true;
-        break;
-      case Code::MuPlus:
-        is_inv = true;
-        break;
-      case Code::MuMinus:
-        is_inv = true;
-        break;
-
-      case Code::Neutron:
-        is_inv = true;
-        break;
-
-      case Code::AntiNeutron:
-        is_inv = true;
-        break;
-
-      default:
-        break;
-    }
-    return is_inv;
-  }
-
-  template <typename TSecondaries>
-  EProcessReturn DoSecondaries(TSecondaries& vS) {
-    auto p = vS.begin();
-    while (p != vS.end()) {
-      const Code pid = p.GetPID();
-      HEPEnergyType energy = p.GetEnergy();
-      cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
-           << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
-           << endl;
-      if (isEmParticle(pid)) {
-        cout << "removing em. particle..." << endl;
-        fEmEnergy += energy;
-        fEmCount += 1;
-        p.Delete();
-      } else if (isInvisible(pid)) {
-        cout << "removing inv. particle..." << endl;
-        fInvEnergy += energy;
-        fInvCount += 1;
-        p.Delete();
-      } else if (isBelowEnergyCut(p)) {
-        cout << "removing low en. particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else if (p.GetTime() > 10_ms) {
-        cout << "removing OLD particle..." << endl;
-        fEnergy += energy;
-        p.Delete();
-      } else {
-        ++p; // next entry in SecondaryView
-      }
-    }
-    return EProcessReturn::eOk;
-  }
-
-  void Init() {
-    fEmEnergy = 0. * 1_GeV;
-    fEmCount = 0;
-    fInvEnergy = 0. * 1_GeV;
-    fInvCount = 0;
-    fEnergy = 0. * 1_GeV;
-    // defineEmParticles();
-  }
-
-  void ShowResults() {
-    cout << " ******************************" << endl
-         << " ParticleCut: " << endl
-         << " energy in em.  component (GeV):  " << fEmEnergy / 1_GeV << endl
-         << " no. of em.  particles injected:  " << fEmCount << endl
-         << " energy in inv. component (GeV):  " << fInvEnergy / 1_GeV << endl
-         << " no. of inv. particles injected:  " << fInvCount << endl
-         << " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
-         << " ******************************" << endl;
-  }
-
-  HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
-  HEPEnergyType GetCutEnergy() const { return fEnergy; }
-  HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
-};
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-  feenableexcept(FE_INVALID);
-  // initialize random number sequence(s)
-  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
-
-  // setup environment, geometry
-  using EnvType = Environment<setup::IEnvironmentModel>;
-  EnvType env;
-  auto& universe = *(env.GetUniverse());
-
-  auto theMedium =
-      EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
-                                  1_km * std::numeric_limits<double>::infinity());
-
-  // fraction of oxygen
-  const float fox = 0.20946;
-  using MyHomogeneousModel = HomogeneousMedium<environment::IMediumModel>;
-  theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_kg / (1_m * 1_m * 1_m),
-      environment::NuclearComposition(
-          std::vector<particles::Code>{particles::Code::Nitrogen,
-                                       particles::Code::Oxygen},
-          std::vector<float>{(float)1. - fox, fox}));
-
-  universe.AddChild(std::move(theMedium));
-
-  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
-
-  // setup processes, decays and interactions
-  tracking_line::TrackingLine tracking;
-
-  const std::vector<particles::Code> trackedHadrons = {
-      particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
-      particles::Code::KMinus, particles::Code::K0Long,  particles::Code::K0Short};
-
-  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
-  process::sibyll::Interaction sibyll;
-  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  process::sibyll::Decay decay(trackedHadrons);
-  // random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-  // process::pythia::Decay decay(trackedHadrons);
-  ProcessCut cut(20_GeV);
-
-  // random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
-  // process::HadronicElasticModel::HadronicElasticInteraction
-  // hadronicElastic(env);
-
-  process::track_writer::TrackWriter trackWriter("tracks.dat");
-  process::energy_loss::EnergyLoss eLoss;
-
-  // assemble all processes into an ordered process list
-  // cut << trackWriter;
-  auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
-
-  // cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
-  // << "\n";
-
-  // 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 = GetNucleusMass(nuclA, nuclZ);
-  const HEPEnergyType E0 = nuclA * 10_TeV;
-  double theta = 0.;
-  double phi = 0.;
-
-  {
-    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 = corsika::stack::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;
-    Point pos(rootCS, 0_m, 0_m,
-              112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
-    stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
-                                 corsika::stack::MomentumVector, geometry::Point,
-                                 units::si::TimeType, unsigned short, unsigned short>{
-        beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
-  }
-
-  // define air shower object, run simulation
-  cascade::Cascade EAS(env, tracking, sequence, stack);
-  EAS.Init();
-  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;
-}

Environment/BaseExponential.h

-/*
- * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_Environment_BaseExponential_h_
-#define _include_Environment_BaseExponential_h_
-
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <limits>
-
-namespace corsika::environment {
-
-  /**
-   * This class provides the grammage/length conversion functionality for
-   * (locally) flat exponential atmospheres.
-   */
-  template <class TDerived>
-  class BaseExponential {
-  protected:
-    units::si::MassDensityType const fRho0;
-    units::si::LengthType const fLambda;
-    units::si::InverseLengthType const fInvLambda;
-    geometry::Point const fP0;
-
-    auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
-
-    // clang-format off
-    /**
-     * For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
-     * direction \f$ \vec{u} \f$ is given by
-     * \f[
-     *   X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right)
-     * \f], where \f$ \varrho_0 \f$ is the density at the starting point.
-     * 
-     * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
-     * \f[
-     *   X = \varrho_0 l;
-     * \f]
-     */
-    // clang-format on
-    units::si::GrammageType IntegratedGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine, units::si::LengthType vL,
-        geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
-      auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
-      auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
-
-      if (uDotA == 0) {
-        return vL * rhoStart;
-      } else {
-        return rhoStart * (fLambda / uDotA) * (exp(uDotA * vL * fInvLambda) - 1);
-      }
-    }
-
-    // clang-format off
-    /**
-     * For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized)
-     * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by
-     * \f[
-     *   l = \begin{cases}
-     *   \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y :=  0 > 1 +
-     *     \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda} 
-     *   \infty & \text{else,}
-     *   \end{cases}
-     * \f] where \f$ \varrho_0 \f$ is the density at the starting point.
-     * 
-     * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
-     * \f[
-     *   l =  \frac{X}{\varrho_0}
-     * \f]
-     */
-    // clang-format on
-    units::si::LengthType ArclengthFromGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::GrammageType vGrammage,
-        geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
-      auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
-      auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
-
-      if (uDotA == 0) {
-        return vGrammage / rhoStart;
-      } else {
-        auto const logArg = vGrammage * fInvLambda * uDotA / rhoStart + 1;
-        if (logArg > 0) {
-          return fLambda / uDotA * log(logArg);
-        } else {
-          return std::numeric_limits<typename decltype(
-                     vGrammage)::value_type>::infinity() *
-                 units::si::meter;
-        }
-      }
-    }
-
-  public:
-    BaseExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
-                    units::si::LengthType vLambda)
-        : fRho0(vRho)
-        , fLambda(vLambda)
-        , fInvLambda(1 / vLambda)
-        , fP0(vP0) {}
-  };
-
-} // namespace corsika::environment
-#endif

Environment/CMakeLists.txt

-set (
-  ENVIRONMENT_HEADERS
-  VolumeTreeNode.h
-  IMediumModel.h
-  NuclearComposition.h
-  HomogeneousMedium.h
-  InhomogeneousMedium.h
-  HomogeneousMedium.h
-  LinearApproximationIntegrator.h
-  DensityFunction.h
-  Environment.h
-  NameModel.h
-  BaseExponential.h
-  FlatExponential.h
-  SlidingPlanarExponential.h
-  )
-
-set (
-  ENVIRONMENT_NAMESPACE
-  corsika/environment
-  )
-
-add_library (CORSIKAenvironment INTERFACE)
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAenvironment
-  INTERFACE
-  CORSIKAgeometry
-  CORSIKAparticles
-  CORSIKAunits
-  CORSIKArandom
-  )
-
-target_include_directories (
-  CORSIKAenvironment
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (
-  TARGETS CORSIKAenvironment
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
-  )
-
-# --------------------
-# code unit testing
-add_executable (testEnvironment testEnvironment.cc)
-
-target_link_libraries (
-  testEnvironment
-  CORSIKAsetup
-  CORSIKAenvironment
-  CORSIKAthirdparty # for catch2
-  )
-CORSIKA_ADD_TEST(testGeometry)

Environment/DensityFunction.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_environment_DensityFunction_h_
-#define _include_environment_DensityFunction_h_
-
-#include <corsika/environment/LinearApproximationIntegrator.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-
-namespace corsika::environment {
-
-  template <class TDerivableRho,
-            template <typename> class TIntegrator = LinearApproximationIntegrator>
-  class DensityFunction
-      : public TIntegrator<DensityFunction<TDerivableRho, TIntegrator>> {
-    friend class TIntegrator<DensityFunction<TDerivableRho, TIntegrator>>;
-
-    TDerivableRho fRho; //!< functor for density
-
-  public:
-    DensityFunction(TDerivableRho rho)
-        : fRho(rho) {}
-
-    corsika::units::si::MassDensityType EvaluateAt(
-        corsika::geometry::Point const& p) const {
-      return fRho(p);
-    }
-  };
-} // namespace corsika::environment
-
-#endif

Environment/Environment.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_environment_Environment_h
-#define _include_environment_Environment_h
-
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Sphere.h>
-#include <limits>
-
-namespace corsika::environment {
-  struct Universe : public corsika::geometry::Sphere {
-    Universe(corsika::geometry::CoordinateSystem const& pCS)
-        : corsika::geometry::Sphere(
-              corsika::geometry::Point{pCS, 0 * corsika::units::si::meter,
-                                       0 * corsika::units::si::meter,
-                                       0 * corsika::units::si::meter},
-              corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
-
-    bool Contains(corsika::geometry::Point const&) const override { return true; }
-  };
-
-  template <typename IEnvironmentModel>
-  class Environment {
-  public:
-    using BaseNodeType = VolumeTreeNode<IEnvironmentModel>;
-
-    Environment()
-        : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
-                                .GetRootCoordinateSystem()}
-        , fUniverse(std::make_unique<BaseNodeType>(
-              std::make_unique<Universe>(fCoordinateSystem))) {}
-
-    // using IEnvironmentModel = corsika::setup::IEnvironmentModel;
-
-    auto& GetUniverse() { return fUniverse; }
-    auto const& GetUniverse() const { return fUniverse; }
-
-    auto const& GetCoordinateSystem() const { return fCoordinateSystem; }
-
-    // factory method for creation of VolumeTreeNodes
-    template <class TVolumeType, typename... TVolumeArgs>
-    static auto CreateNode(TVolumeArgs&&... args) {
-      static_assert(std::is_base_of_v<corsika::geometry::Volume, TVolumeType>,
-                    "unusable type provided, needs to be derived from "
-                    "\"corsika::geometry::Volume\"");
-
-      return std::make_unique<BaseNodeType>(
-          std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
-    }
-
-  private:
-    corsika::geometry::CoordinateSystem const& fCoordinateSystem;
-    typename BaseNodeType::VTNUPtr fUniverse;
-  };
-
-  // using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
-  // using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
-
-} // namespace corsika::environment
-
-#endif

Environment/FlatExponential.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_Environment_FlatExponential_h_
-#define _include_Environment_FlatExponential_h_
-
-#include <corsika/environment/BaseExponential.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  //clang-format off
-  /**
-   * flat exponential density distribution with
-   * \f[
-   *  \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
-   *    \vec{a} \right).
-   * \f]
-   * \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
-   * with the scale parameter \f$ \lambda \f$.
-   */
-  //clang-format on
-  template <class T>
-  class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
-    geometry::Vector<units::si::dimensionless_d> const fAxis;
-    NuclearComposition const fNuclComp;
-
-    using Base = BaseExponential<FlatExponential<T>>;
-
-  public:
-    FlatExponential(geometry::Point const& vP0,
-                    geometry::Vector<units::si::dimensionless_d> const& vAxis,
-                    units::si::MassDensityType vRho, units::si::LengthType vLambda,
-                    NuclearComposition vNuclComp)
-        : Base(vP0, vRho, vLambda)
-        , fAxis(vAxis)
-        , fNuclComp(vNuclComp) {}
-
-    units::si::MassDensityType GetMassDensity(geometry::Point const& vP) const override {
-      return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis));
-    }
-
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    units::si::GrammageType IntegratedGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::LengthType vTo) const override {
-      return Base::IntegratedGrammage(vLine, vTo, fAxis);
-    }
-
-    units::si::LengthType ArclengthFromGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::GrammageType vGrammage) const override {
-      return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
-    }
-  };
-} // namespace corsika::environment
-#endif

Environment/HomogeneousMedium.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_HomogeneousMedium_h_
-#define _include_HomogeneousMedium_h_
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cassert>
-
-/**
- * a homogeneous medium
- */
-
-namespace corsika::environment {
-
-  template <class T>
-  class HomogeneousMedium : public T {
-    corsika::units::si::MassDensityType const fDensity;
-    NuclearComposition const fNuclComp;
-
-  public:
-    HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
-                      NuclearComposition pNuclComp)
-        : fDensity(pDensity)
-        , fNuclComp(pNuclComp) {}
-
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const override {
-      return fDensity;
-    }
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::LengthType pTo) const override {
-      using namespace corsika::units::si;
-      return pTo * fDensity;
-    }
-
-    corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::GrammageType pGrammage) const override {
-      return pGrammage / fDensity;
-    }
-  };
-
-} // namespace corsika::environment
-#endif