SCIENTIFIC_AUTHORS

-This is currently derived from
-https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/wikis/ICRC2019
-and will be revised whenever needed:
-
-- Luisa Arrabito (1)
-- Dominik Baack (2)
-- Johan Bregeon (1)
-- Hans Dembinski (3)
-- Ralph Engel (4)
-- Dieter Heck (4)
-- Tim Huege (4,5)
-- Lukas Nellen (6)
-- Tanguy Pierog (4)
-- Maximilian Reininghaus (4,7)
-- Felix Riehn (8)
-- Ralf Ulrich (4)
-- Darko Veberic (4)
-
-Affiliations:
-
-- University of Montpellier, France
-- TU Dortmund University, Germany
-- Max Planck Institute for Nuclear Physics, Heidelberg, Germany
-- Karlsruhe Institute of Technology (KIT), Germany
-- Vrije Universiteit Brussel, Belgium
-- National Autonomous University of Mexico (UNAM)
-- Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina
-- Laboratory of Instrumentation and Experimental Particles (LIP), Portugal

Stack/DummyStack/testDummyStack.cc

-/*
- * (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/stack/dummy/DummyStack.h>
-
-using namespace corsika;
-using namespace corsika::stack;
-
-#include <catch2/catch.hpp>
-
-#include <tuple>
-
-TEST_CASE("DummyStack", "[stack]") {
-
-  using TestStack = dummy::DummyStack;
-
-  dummy::NoData noData;
-
-  SECTION("write node") {
-
-    TestStack s;
-    s.AddParticle(std::tuple<dummy::NoData>{noData});
-    CHECK(s.getEntries() == 1);
-  }
-
-  SECTION("stack fill and cleanup") {
-
-    TestStack s;
-    // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
-    for (int i = 0; i < 99; ++i) { s.AddParticle(std::tuple<dummy::NoData>{noData}); }
-
-    CHECK(s.getEntries() == 99);
-    for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete();
-    CHECK(s.getEntries() == 0);
-  }
-}

Stack/GeometryNodeStackExtension/CMakeLists.txt

-set (GeometryNodeStackExtension_HEADERS GeometryNodeStackExtension.h)
-set (GeometryNodeStackExtension_NAMESPACE corsika/stack/node)
-
-add_library (GeometryNodeStackExtension INTERFACE)
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (GeometryNodeStackExtension ${GeometryNodeStackExtension_NAMESPACE} ${GeometryNodeStackExtension_HEADERS})
-
-target_link_libraries (
-  GeometryNodeStackExtension
-  INTERFACE
-  CORSIKAstackinterface
-  CORSIKAunits
-  CORSIKAparticles
-  CORSIKAgeometry
-  )
-
-target_include_directories (
-  GeometryNodeStackExtension
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (
-  FILES
-  ${GeometryNodeStackExtension_HEADERS}
-  DESTINATION
-  include/${GeometryNodeStackExtension_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
- CORSIKA_ADD_TEST(testGeometryNodeStackExtension)
- target_link_libraries (
-   testGeometryNodeStackExtension
-   GeometryNodeStackExtension
-   CORSIKAtesting
-   )

Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h

-/*
- * (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.
- */
-
-#pragma once
-
-#include <corsika/logging/Logging.h>
-#include <corsika/stack/Stack.h>
-
-#include <tuple>
-#include <utility>
-#include <vector>
-
-namespace corsika::stack::node {
-
-  /**
-   * @class GeometryDataInterface
-   *
-   * corresponding defintion of a stack-readout object, the iteractor
-   * dereference operator will deliver access to these function
-  // defintion of a stack-readout object, the iteractor dereference
-  // operator will deliver access to these function
-   */
-  template <typename T, typename TEnvType>
-  class GeometryDataInterface : public T {
-
-  protected:
-    using T::GetStack;
-    using T::GetStackData;
-
-  public:
-    using T::GetIndex;
-    using BaseNodeType = typename TEnvType::BaseNodeType;
-
-  public:
-    // default version for particle-creation from input data
-    void SetParticleData(const std::tuple<BaseNodeType const*> v) {
-      SetNode(std::get<0>(v));
-    }
-    void SetParticleData(GeometryDataInterface& parent,
-                         const std::tuple<BaseNodeType const*>) {
-      SetNode(parent.GetNode()); // copy Node from parent particle!
-    }
-    void SetParticleData() { SetNode(nullptr); }
-    void SetParticleData(GeometryDataInterface& parent) {
-      SetNode(parent.GetNode()); // copy Node from parent particle!
-    }
-
-    std::string as_string() const { return fmt::format("node={}", fmt::ptr(GetNode())); }
-
-    void SetNode(BaseNodeType const* v) { GetStackData().SetNode(GetIndex(), v); }
-    BaseNodeType const* GetNode() const { return GetStackData().GetNode(GetIndex()); }
-  };
-
-  // definition of stack-data object to store geometry information
-  template <typename TEnvType>
-
-  /**
-   * @class GeometryData
-   *
-   * definition of stack-data object to store geometry information
-   */
-  class GeometryData {
-
-  public:
-    using BaseNodeType = typename TEnvType::BaseNodeType;
-
-    // these functions are needed for the Stack interface
-    void Clear() { fNode.clear(); }
-    unsigned int GetSize() const { return fNode.size(); }
-    unsigned int GetCapacity() const { return fNode.size(); }
-    void Copy(const int i1, const int i2) { fNode[i2] = fNode[i1]; }
-    void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); }
-
-    // custom data access function
-    void SetNode(const int i, BaseNodeType const* v) { fNode[i] = v; }
-    BaseNodeType const* GetNode(const int i) const { return fNode[i]; }
-
-    // these functions are also needed by the Stack interface
-    void IncrementSize() { fNode.push_back(nullptr); }
-    void DecrementSize() {
-      if (fNode.size() > 0) { fNode.pop_back(); }
-    }
-
-    // custom private data section
-  private:
-    std::vector<const BaseNodeType*> fNode;
-  };
-
-  template <typename T, typename TEnv>
-  struct MakeGeometryDataInterface {
-    typedef GeometryDataInterface<T, TEnv> type;
-  };
-
-} // namespace corsika::stack::node

Stack/GeometryNodeStackExtension/testGeometryNodeStackExtension.cc

-/*
- * (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/stack/CombinedStack.h>
-#include <corsika/stack/dummy/DummyStack.h>
-#include <corsika/stack/node/GeometryNodeStackExtension.h>
-
-using namespace corsika;
-using namespace corsika::stack;
-
-#include <catch2/catch.hpp>
-
-#include <iostream>
-using namespace std;
-
-// this is our dummy environment, it only knows its trivial BaseNodeType
-class DummyEnv {
-public:
-  typedef int BaseNodeType;
-};
-
-// the GeometryNode stack needs to know the type of geometry-nodes from the DummyEnv:
-template <typename TStackIter>
-using DummyGeometryDataInterface =
-    typename corsika::stack::node::MakeGeometryDataInterface<TStackIter, DummyEnv>::type;
-
-// combine dummy stack with geometry information for tracking
-template <typename TStackIter>
-using StackWithGeometryInterface =
-    corsika::stack::CombinedParticleInterface<dummy::DummyStack::MPIType,
-                                              DummyGeometryDataInterface, TStackIter>;
-
-using TestStack =
-    corsika::stack::CombinedStack<typename stack::dummy::DummyStack::StackImpl,
-                                  stack::node::GeometryData<DummyEnv>,
-                                  StackWithGeometryInterface>;
-
-TEST_CASE("GeometryNodeStackExtension", "[stack]") {
-
-  dummy::NoData noData;
-
-  SECTION("write node") {
-
-    const int data = 5;
-
-    TestStack s;
-    s.AddParticle(std::make_tuple(noData), std::tuple<const int*>{&data});
-
-    CHECK(s.getEntries() == 1);
-  }
-
-  SECTION("write/read node") {
-    const int data = 15;
-
-    TestStack s;
-    auto p = s.AddParticle(std::make_tuple(noData));
-    p.SetNode(&data);
-    CHECK(s.getEntries() == 1);
-
-    const auto pout = s.GetNextParticle();
-    CHECK(*(pout.GetNode()) == 15);
-  }
-
-  SECTION("stack fill and cleanup") {
-
-    const int data = 16;
-
-    TestStack s;
-    // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
-    for (int i = 0; i < 99; ++i) {
-      auto p = s.AddParticle(std::tuple<dummy::NoData>{noData});
-      p.SetNode(&data);
-    }
-
-    CHECK(s.getEntries() == 99);
-    double v = 0;
-    for (int i = 0; i < 99; ++i) {
-      auto p = s.GetNextParticle();
-      v += *(p.GetNode());
-      p.Delete();
-    }
-    CHECK(v == 99 * data);
-    CHECK(s.getEntries() == 0);
-  }
-}

Stack/History/CMakeLists.txt

-set (
-  HISTORY_HEADERS
-  EventType.hpp
-  Event.hpp
-  HistorySecondaryProducer.hpp
-  HistoryObservationPlane.hpp
-  HistoryStackExtension.hpp
-  SecondaryParticle.hpp
-  )
-
-set (
-  HISTORY_NAMESPACE
-  corsika/stack/history
-  )
-
-if (WITH_HISTORY)
-  set (
-    HISTORY_SOURCES
-    HistoryObservationPlane.cpp
-    )
-  add_library (CORSIKAhistory STATIC ${HISTORY_SOURCES})
-  set (IS_INTERFACE "")
-else (WITH_HISTORY)
-  add_library (CORSIKAhistory INTERFACE)
-  set (IS_INTERFACE "INTERFACE")
-endif (WITH_HISTORY)
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAhistory ${HISTORY_NAMESPACE} ${HISTORY_HEADERS})
-
-target_link_libraries (
-  CORSIKAhistory
-  ${IS_INTERFACE}
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAprocesssequence # for HistoryObservationPlane
-  CORSIKAsetup # for HistoryObservationPlane
-  CORSIKAlogging
-  SuperStupidStack
-  NuclearStackExtension # for testHistoryView.cc
-  C8::ext::boost # for HistoryObservationPlane
-  )
-
-target_include_directories (
-  CORSIKAhistory
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  FILES ${HISTORY_HEADERS} 
-  DESTINATION include/${HISTORY_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-CORSIKA_ADD_TEST(testHistoryStack)
-target_link_libraries (
-  testHistoryStack
-  CORSIKAhistory
-  CORSIKAtesting
-  DummyStack
-  )
-CORSIKA_ADD_TEST(testHistoryView)
-target_link_libraries (
-  testHistoryView
-  CORSIKAhistory
-  CORSIKAtesting
-  )

Stack/History/HistoryObservationPlane.cpp

-/*
- * (c) Copyright 2020 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/logging/Logging.h>
-#include <corsika/stack/history/HistoryObservationPlane.hpp>
-
-#include <boost/histogram/ostream.hpp>
-
-#include <iomanip>
-#include <iostream>
-
-using namespace corsika::units::si;
-using namespace corsika::history;
-using namespace corsika;
-
-HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
-                                                 geometry::Plane const& obsPlane,
-                                                 bool deleteOnHit)
-    : stack_{stack}
-    , plane_{obsPlane}
-    , deleteOnHit_{deleteOnHit} {}
-
-corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
-    setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
-  TimeType const timeOfIntersection =
-      (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
-      trajectory.GetLine().GetV0().dot(plane_.GetNormal());
-
-  if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
-
-  if (plane_.IsAbove(trajectory.GetLine().GetR0()) ==
-      plane_.IsAbove(trajectory.GetPosition(1))) {
-    return process::EProcessReturn::eOk;
-  }
-
-  C8LOG_DEBUG(fmt::format("HistoryObservationPlane: Particle detected: pid={}",
-                          particle.GetPID()));
-
-  auto const pid = particle.GetPID();
-  if (particles::IsMuon(pid)) { fillHistoryHistogram(particle); }
-
-  if (deleteOnHit_) {
-    return process::EProcessReturn::eParticleAbsorbed;
-  } else {
-    return process::EProcessReturn::eOk;
-  }
-}
-
-LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
-                                                  setup::Trajectory const& trajectory) {
-  TimeType const timeOfIntersection =
-      (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
-      trajectory.GetLine().GetV0().dot(plane_.GetNormal());
-
-  if (timeOfIntersection < TimeType::zero()) {
-    return std::numeric_limits<double>::infinity() * 1_m;
-  }
-
-  auto const pointOfIntersection = trajectory.GetLine().GetPosition(timeOfIntersection);
-  return (trajectory.GetLine().GetR0() - pointOfIntersection).norm() * 1.0001;
-}
-
-void HistoryObservationPlane::fillHistoryHistogram(
-    setup::Stack::ParticleType const& muon) {
-  double const muon_energy = muon.GetEnergy() / 1_GeV;
-
-  int genctr{0};
-  Event const* event = muon.GetEvent().get();
-  while (event) {
-    auto const projectile = stack_.cfirst() + event->projectileIndex();
-    if (event->eventType() == EventType::Interaction) {
-      genctr++;
-      double const projEnergy = projectile.GetEnergy() / 1_GeV;
-      int const pdg = static_cast<int>(particles::GetPDG(projectile.GetPID()));
-
-      histogram_(muon_energy, projEnergy, pdg);
-    }
-    event = event->parentEvent().get(); // projectile.GetEvent().get();
-  }
-}

Stack/History/testHistoryView.cc

-/*
- * (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/stack/history/Event.hpp>
-#include <corsika/stack/history/HistorySecondaryProducer.hpp>
-#include <corsika/stack/history/HistoryStackExtension.hpp>
-
-#include <corsika/stack/CombinedStack.h>
-#include <corsika/stack/dummy/DummyStack.h>
-#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
-
-#include <corsika/logging/Logging.h>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::geometry;
-using namespace corsika::units::si;
-
-/**
-   Need to replicate setup::SetupStack in a maximally simplified
-   way, but with real particle data
- */
-
-// combine dummy stack with geometry information for tracking
-template <typename TStackIter>
-using StackWithHistoryInterface = corsika::stack::CombinedParticleInterface<
-    stack::nuclear_extension::ParticleDataStack::MPIType,
-    history::HistoryEventDataInterface, TStackIter>;
-
-using TestStack = corsika::stack::CombinedStack<
-    typename stack::nuclear_extension::ParticleDataStack::StackImpl,
-    history::HistoryEventData, StackWithHistoryInterface>;
-
-/*
-    See Issue 161
-
-    unfortunately clang does not support this in the same way (yet) as
-    gcc, so we have to distinguish here. If clang cataches up, we
-    could remove the clang branch here and also in
-    corsika::Cascade. The gcc code is much more generic and
-    universal. If we could do the gcc version, we won't had to define
-    StackView globally, we could do it with MakeView whereever it is
-    actually needed. Keep an eye on this!
-  */
-#if defined(__clang__)
-using TheTestStackView = corsika::stack::SecondaryView<typename TestStack::StackImpl,
-                                                       StackWithHistoryInterface,
-                                                       history::HistorySecondaryProducer>;
-#elif defined(__GNUC__) || defined(__GNUG__)
-using TheTestStackView =
-    corsika::stack::MakeView<TestStack, history::HistorySecondaryProducer>::type;
-#endif
-
-using TestStackView = TheTestStackView;
-
-template <typename Event>
-int count_generations(Event const* event) {
-  int genCounter = 0;
-  while (event) {
-    event = event->parentEvent().get();
-    genCounter++;
-  }
-
-  return genCounter;
-}
-
-TEST_CASE("HistoryStackExtension", "[stack]") {
-
-  logging::SetLevel(logging::level::debug);
-
-  geometry::CoordinateSystem& dummyCS =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  // in this test we only use one singel stack !
-  TestStack stack;
-
-  // add primary particle
-  auto p0 = stack.AddParticle(
-      std::make_tuple(particles::Code::Electron, 1.5_GeV,
-                      corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
-                      Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
-
-  CHECK(stack.getEntries() == 1);
-  corsika::history::EventPtr evt = p0.GetEvent();
-  CHECK(evt == nullptr);
-  CHECK(count_generations(evt.get()) == 0);
-
-  SECTION("interface test, view") {
-
-    // add secondaries, 1st generation
-    TestStackView hview0(p0);
-
-    auto const ev0 = p0.GetEvent();
-    CHECK(ev0 == nullptr);
-
-    C8LOG_DEBUG("loop VIEW");
-
-    // add 5 secondaries
-    for (int i = 0; i < 5; ++i) {
-      auto sec = hview0.AddSecondary(
-          std::make_tuple(particles::Code::Electron, 1.5_GeV,
-                          corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
-                          Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
-
-      CHECK(sec.GetParentEventIndex() == i);
-      CHECK(sec.GetEvent() != nullptr);
-      CHECK(sec.GetEvent()->parentEvent() == nullptr);
-      CHECK(count_generations(sec.GetEvent().get()) == 1);
-    }
-
-    // read 1st genertion particle particle
-    auto p1 = stack.GetNextParticle();
-    CHECK(count_generations(p1.GetEvent().get()) == 1);
-
-    TestStackView hview1(p1);
-
-    auto const ev1 = p1.GetEvent();
-
-    // add second generation of secondaries
-    // add 10 secondaries
-    for (int i = 0; i < 10; ++i) {
-      auto sec = hview1.AddSecondary(
-          std::make_tuple(particles::Code::Electron, 1.5_GeV,
-                          corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
-                          Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
-
-      CHECK(sec.GetParentEventIndex() == i);
-      CHECK(sec.GetEvent()->parentEvent() == ev1);
-      CHECK(sec.GetEvent()->parentEvent()->parentEvent() == ev0);
-
-      CHECK(count_generations(sec.GetEvent().get()) == 2);
-
-      const auto org_projectile = stack.at(sec.GetEvent()->projectileIndex());
-      CHECK(org_projectile.GetEvent() == sec.GetEvent()->parentEvent());
-    }
-
-    // read 2nd genertion particle particle
-    auto p2 = stack.GetNextParticle();
-
-    TestStackView hview2(p2);
-
-    auto const ev2 = p2.GetEvent();
-
-    // add third generation of secondaries
-    // add 15 secondaries
-    for (int i = 0; i < 15; ++i) {
-      C8LOG_TRACE("loop, view: " + std::to_string(i));
-
-      auto sec = hview2.AddSecondary(
-          std::make_tuple(particles::Code::Electron, 1.5_GeV,
-                          corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
-                          Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
-      C8LOG_TRACE("loop, ---- ");
-
-      CHECK(sec.GetParentEventIndex() == i);
-      CHECK(sec.GetEvent()->parentEvent() == ev2);
-      CHECK(sec.GetEvent()->parentEvent()->parentEvent() == ev1);
-      CHECK(sec.GetEvent()->parentEvent()->parentEvent()->parentEvent() == ev0);
-
-      CHECK(count_generations(sec.GetEvent().get()) == 3);
-    }
-  }
-
-  SECTION("also test projectile access") {
-
-    C8LOG_TRACE("projectile test");
-
-    // add secondaries, 1st generation
-    TestStackView hview0(p0);
-    auto proj0 = hview0.GetProjectile();
-    auto const ev0 = p0.GetEvent();
-    CHECK(ev0 == nullptr);
-
-    C8LOG_TRACE("loop");
-
-    // add 5 secondaries
-    for (int i = 0; i < 5; ++i) {
-      C8LOG_TRACE("loop " + std::to_string(i));
-      auto sec = proj0.AddSecondary(
-          std::make_tuple(particles::Code::Electron, 1.5_GeV,
-                          corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
-                          Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
-
-      CHECK(sec.GetParentEventIndex() == i);
-      CHECK(sec.GetEvent() != nullptr);
-      CHECK(sec.GetEvent()->parentEvent() == nullptr);
-    }
-    CHECK(stack.getEntries() == 6);
-
-    // read 1st genertion particle particle
-    auto p1 = stack.GetNextParticle();
-
-    TestStackView hview1(p1);
-    auto proj1 = hview1.GetProjectile();
-    auto const ev1 = p1.GetEvent();
-
-    // add second generation of secondaries
-    // add 10 secondaries
-    for (int i = 0; i < 10; ++i) {
-      auto sec = proj1.AddSecondary(
-          std::make_tuple(particles::Code::Electron, 1.5_GeV,
-                          corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
-                          Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
-
-      CHECK(sec.GetParentEventIndex() == i);
-      CHECK(sec.GetEvent()->parentEvent() == ev1);
-      CHECK(sec.GetEvent()->secondaries().size() == i + 1);
-      CHECK(sec.GetEvent()->parentEvent()->parentEvent() == ev0);
-
-      const auto org_projectile = stack.at(sec.GetEvent()->projectileIndex());
-      CHECK(org_projectile.GetEvent() == sec.GetEvent()->parentEvent());
-    }
-    CHECK(stack.getEntries() == 16);
-  }
-}

ThirdParty/cnpy/CMakeLists.txt

-set (
-  CNPY_SOURCES
-  cnpy.cpp
-  )
-
-set (
-  CNPY_HEADERS
-  cnpy.hpp
-  )
-
-set (
-  CNPY_NAMESPACE
-  corsika/third_party/cnpy
-  )
-
-add_library (cnpy STATIC ${CNPY_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (cnpy ${CNPY_NAMESPACE} ${CNPY_HEADERS})
-
-set_target_properties (
-  cnpy
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  cnpy
-  PUBLIC
-  ZLIB::ZLIB
-  )
-
-target_include_directories (
-  cnpy 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS cnpy
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  )
-

ThirdParty/cnpy/cnpy.cpp

-//Copyright (C) 2011  Carl Rogers
-//Released under MIT License
-//license available in LICENSE file, or at http://www.opensource.org/licenses/mit-license.php
-
-#include "cnpy.hpp"
-#include <complex>
-#include <cstdlib>
-#include <algorithm>
-#include <cstring>
-#include <iomanip>
-#include <stdint.h>
-#include <stdexcept>
-#include <regex>
-
-char cnpy::BigEndianTest() {
-    int x = 1;
-    return (((char *)&x)[0]) ? '<' : '>';
-}
-
-char cnpy::map_type(const std::type_info& t)
-{
-    if(t == typeid(float) ) return 'f';
-    if(t == typeid(double) ) return 'f';
-    if(t == typeid(long double) ) return 'f';
-
-    if(t == typeid(int) ) return 'i';
-    if(t == typeid(char) ) return 'i';
-    if(t == typeid(short) ) return 'i';
-    if(t == typeid(long) ) return 'i';
-    if(t == typeid(long long) ) return 'i';
-
-    if(t == typeid(unsigned char) ) return 'u';
-    if(t == typeid(unsigned short) ) return 'u';
-    if(t == typeid(unsigned long) ) return 'u';
-    if(t == typeid(unsigned long long) ) return 'u';
-    if(t == typeid(unsigned int) ) return 'u';
-
-    if(t == typeid(bool) ) return 'b';
-
-    if(t == typeid(std::complex<float>) ) return 'c';
-    if(t == typeid(std::complex<double>) ) return 'c';
-    if(t == typeid(std::complex<long double>) ) return 'c';
-
-    else return '?';
-}
-
-template<> std::vector<char>& cnpy::operator+=(std::vector<char>& lhs, const std::string rhs) {
-    lhs.insert(lhs.end(),rhs.begin(),rhs.end());
-    return lhs;
-}
-
-template<> std::vector<char>& cnpy::operator+=(std::vector<char>& lhs, const char* rhs) {
-    //write in little endian
-    size_t len = strlen(rhs);
-    lhs.reserve(len);
-    for(size_t byte = 0; byte < len; byte++) {
-        lhs.push_back(rhs[byte]);
-    }
-    return lhs;
-}
-
-void cnpy::parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector<size_t>& shape, bool& fortran_order) {
-    //std::string magic_string(buffer,6);
-    [[maybe_unused]] uint8_t major_version = *reinterpret_cast<uint8_t*>(buffer+6);
-    [[maybe_unused]] uint8_t minor_version = *reinterpret_cast<uint8_t*>(buffer+7);
-    uint16_t header_len = *reinterpret_cast<uint16_t*>(buffer+8);
-    std::string header(reinterpret_cast<char*>(buffer+9),header_len);
-
-    size_t loc1, loc2;
-
-    //fortran order
-    loc1 = header.find("fortran_order")+16;
-    fortran_order = (header.substr(loc1,4) == "True" ? true : false);
-
-    //shape
-    loc1 = header.find("(");
-    loc2 = header.find(")");
-
-    std::regex num_regex("[0-9][0-9]*");
-    std::smatch sm;
-    shape.clear();
-
-    std::string str_shape = header.substr(loc1+1,loc2-loc1-1);
-    while(std::regex_search(str_shape, sm, num_regex)) {
-        shape.push_back(std::stoi(sm[0].str()));
-        str_shape = sm.suffix().str();
-    }
-
-    //endian, word size, data type
-    //byte order code | stands for not applicable. 
-    //not sure when this applies except for byte array
-    loc1 = header.find("descr")+9;
-    bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false);
-    assert(littleEndian);
-
-    //char type = header[loc1+1];
-    //assert(type == map_type(T));
-
-    std::string str_ws = header.substr(loc1+2);
-    loc2 = str_ws.find("'");
-    word_size = atoi(str_ws.substr(0,loc2).c_str());
-}
-
-void cnpy::parse_npy_header(FILE* fp, size_t& word_size, std::vector<size_t>& shape, bool& fortran_order) {  
-    char buffer[256];
-    size_t res = fread(buffer,sizeof(char),11,fp);       
-    if(res != 11)
-        throw std::runtime_error("parse_npy_header: failed fread");
-    std::string header = fgets(buffer,256,fp);
-    assert(header[header.size()-1] == '\n');
-
-    size_t loc1, loc2;
-
-    //fortran order
-    loc1 = header.find("fortran_order");
-    if (loc1 == std::string::npos)
-        throw std::runtime_error("parse_npy_header: failed to find header keyword: 'fortran_order'");
-    loc1 += 16;
-    fortran_order = (header.substr(loc1,4) == "True" ? true : false);
-
-    //shape
-    loc1 = header.find("(");
-    loc2 = header.find(")");
-    if (loc1 == std::string::npos || loc2 == std::string::npos)
-        throw std::runtime_error("parse_npy_header: failed to find header keyword: '(' or ')'");
-
-    std::regex num_regex("[0-9][0-9]*");
-    std::smatch sm;
-    shape.clear();
-
-    std::string str_shape = header.substr(loc1+1,loc2-loc1-1);
-    while(std::regex_search(str_shape, sm, num_regex)) {
-        shape.push_back(std::stoi(sm[0].str()));
-        str_shape = sm.suffix().str();
-    }
-
-    //endian, word size, data type
-    //byte order code | stands for not applicable. 
-    //not sure when this applies except for byte array
-    loc1 = header.find("descr");
-    if (loc1 == std::string::npos)
-        throw std::runtime_error("parse_npy_header: failed to find header keyword: 'descr'");
-    loc1 += 9;
-    bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false);
-    assert(littleEndian);
-
-    //char type = header[loc1+1];
-    //assert(type == map_type(T));
-
-    std::string str_ws = header.substr(loc1+2);
-    loc2 = str_ws.find("'");
-    word_size = atoi(str_ws.substr(0,loc2).c_str());
-}
-
-void cnpy::parse_zip_footer(FILE* fp, uint16_t& nrecs, size_t& global_header_size, size_t& global_header_offset)
-{
-    std::vector<char> footer(22);
-    fseek(fp,-22,SEEK_END);
-    size_t res = fread(&footer[0],sizeof(char),22,fp);
-    if(res != 22)
-        throw std::runtime_error("parse_zip_footer: failed fread");
-
-    uint16_t disk_no, disk_start, nrecs_on_disk, comment_len;
-    disk_no = *(uint16_t*) &footer[4];
-    disk_start = *(uint16_t*) &footer[6];
-    nrecs_on_disk = *(uint16_t*) &footer[8];
-    nrecs = *(uint16_t*) &footer[10];
-    global_header_size = *(uint32_t*) &footer[12];
-    global_header_offset = *(uint32_t*) &footer[16];
-    comment_len = *(uint16_t*) &footer[20];
-
-    assert(disk_no == 0);
-    assert(disk_start == 0);
-    assert(nrecs_on_disk == nrecs);
-    assert(comment_len == 0);
-}
-
-cnpy::NpyArray load_the_npy_file(FILE* fp) {
-    std::vector<size_t> shape;
-    size_t word_size;
-    bool fortran_order;
-    cnpy::parse_npy_header(fp,word_size,shape,fortran_order);
-
-    cnpy::NpyArray arr(shape, word_size, fortran_order);
-    size_t nread = fread(arr.data<char>(),1,arr.num_bytes(),fp);
-    if(nread != arr.num_bytes())
-        throw std::runtime_error("load_the_npy_file: failed fread");
-    return arr;
-}
-
-cnpy::NpyArray load_the_npz_array(FILE* fp, uint32_t compr_bytes, uint32_t uncompr_bytes) {
-
-    std::vector<unsigned char> buffer_compr(compr_bytes);
-    std::vector<unsigned char> buffer_uncompr(uncompr_bytes);
-    size_t nread = fread(&buffer_compr[0],1,compr_bytes,fp);
-    if(nread != compr_bytes)
-        throw std::runtime_error("load_the_npy_file: failed fread");
-
-    [[maybe_unused]] int err;
-    z_stream d_stream;
-
-    d_stream.zalloc = Z_NULL;
-    d_stream.zfree = Z_NULL;
-    d_stream.opaque = Z_NULL;
-    d_stream.avail_in = 0;
-    d_stream.next_in = Z_NULL;
-    err = inflateInit2(&d_stream, -MAX_WBITS);
-
-    d_stream.avail_in = compr_bytes;
-    d_stream.next_in = &buffer_compr[0];
-    d_stream.avail_out = uncompr_bytes;
-    d_stream.next_out = &buffer_uncompr[0];
-
-    err = inflate(&d_stream, Z_FINISH);
-    err = inflateEnd(&d_stream);
-
-    std::vector<size_t> shape;
-    size_t word_size;
-    bool fortran_order;
-    cnpy::parse_npy_header(&buffer_uncompr[0],word_size,shape,fortran_order);
-
-    cnpy::NpyArray array(shape, word_size, fortran_order);
-
-    size_t offset = uncompr_bytes - array.num_bytes();
-    memcpy(array.data<unsigned char>(),&buffer_uncompr[0]+offset,array.num_bytes());
-
-    return array;
-}
-
-cnpy::npz_t cnpy::npz_load(std::string fname) {
-    FILE* fp = fopen(fname.c_str(),"rb");
-
-    if(!fp) {
-        throw std::runtime_error("npz_load: Error! Unable to open file "+fname+"!");
-    }
-
-    cnpy::npz_t arrays;  
-
-    while(1) {
-        std::vector<char> local_header(30);
-        size_t headerres = fread(&local_header[0],sizeof(char),30,fp);
-        if(headerres != 30)
-            throw std::runtime_error("npz_load: failed fread");
-
-        //if we've reached the global header, stop reading
-        if(local_header[2] != 0x03 || local_header[3] != 0x04) break;
-
-        //read in the variable name
-        uint16_t name_len = *(uint16_t*) &local_header[26];
-        std::string varname(name_len,' ');
-        size_t vname_res = fread(&varname[0],sizeof(char),name_len,fp);
-        if(vname_res != name_len)
-            throw std::runtime_error("npz_load: failed fread");
-
-        //erase the lagging .npy        
-        varname.erase(varname.end()-4,varname.end());
-
-        //read in the extra field
-        uint16_t extra_field_len = *(uint16_t*) &local_header[28];
-        if(extra_field_len > 0) {
-            std::vector<char> buff(extra_field_len);
-            size_t efield_res = fread(&buff[0],sizeof(char),extra_field_len,fp);
-            if(efield_res != extra_field_len)
-                throw std::runtime_error("npz_load: failed fread");
-        }
-
-        uint16_t compr_method = *reinterpret_cast<uint16_t*>(&local_header[0]+8);
-        uint32_t compr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+18);
-        uint32_t uncompr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+22);
-
-        if(compr_method == 0) {arrays[varname] = load_the_npy_file(fp);}
-        else {arrays[varname] = load_the_npz_array(fp,compr_bytes,uncompr_bytes);}
-    }
-
-    fclose(fp);
-    return arrays;  
-}
-
-cnpy::NpyArray cnpy::npz_load(std::string fname, std::string varname) {
-    FILE* fp = fopen(fname.c_str(),"rb");
-
-    if(!fp) throw std::runtime_error("npz_load: Unable to open file "+fname);
-
-    while(1) {
-        std::vector<char> local_header(30);
-        size_t header_res = fread(&local_header[0],sizeof(char),30,fp);
-        if(header_res != 30)
-            throw std::runtime_error("npz_load: failed fread");
-
-        //if we've reached the global header, stop reading
-        if(local_header[2] != 0x03 || local_header[3] != 0x04) break;
-
-        //read in the variable name
-        uint16_t name_len = *(uint16_t*) &local_header[26];
-        std::string vname(name_len,' ');
-        size_t vname_res = fread(&vname[0],sizeof(char),name_len,fp);      
-        if(vname_res != name_len)
-            throw std::runtime_error("npz_load: failed fread");
-        vname.erase(vname.end()-4,vname.end()); //erase the lagging .npy
-
-        //read in the extra field
-        uint16_t extra_field_len = *(uint16_t*) &local_header[28];
-        fseek(fp,extra_field_len,SEEK_CUR); //skip past the extra field
-        
-        uint16_t compr_method = *reinterpret_cast<uint16_t*>(&local_header[0]+8);
-        uint32_t compr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+18);
-        uint32_t uncompr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+22);
-
-        if(vname == varname) {
-            NpyArray array  = (compr_method == 0) ? load_the_npy_file(fp) : load_the_npz_array(fp,compr_bytes,uncompr_bytes);
-            fclose(fp);
-            return array;
-        }
-        else {
-            //skip past the data
-            uint32_t size = *(uint32_t*) &local_header[22];
-            fseek(fp,size,SEEK_CUR);
-        }
-    }
-
-    fclose(fp);
-
-    //if we get here, we haven't found the variable in the file
-    throw std::runtime_error("npz_load: Variable name "+varname+" not found in "+fname);
-}
-
-cnpy::NpyArray cnpy::npy_load(std::string fname) {
-
-    FILE* fp = fopen(fname.c_str(), "rb");
-
-    if(!fp) throw std::runtime_error("npy_load: Unable to open file "+fname);
-
-    NpyArray arr = load_the_npy_file(fp);
-
-    fclose(fp);
-    return arr;
-}
-
-
-

ThirdParty/cnpy/cnpy.hpp

-//Copyright (C) 2011  Carl Rogers
-//Released under MIT License
-//license available in LICENSE file, or at http://www.opensource.org/licenses/mit-license.php
-
-#pragma once
-
-#include<string>
-#include<stdexcept>
-#include<sstream>
-#include<vector>
-#include<cstdio>
-#include<typeinfo>
-#include<iostream>
-#include<cassert>
-#include<zlib.h>
-#include<map>
-#include<memory>
-#include<stdint.h>
-#include<numeric>
-
-namespace cnpy {
-
-    struct NpyArray {
-        NpyArray(const std::vector<size_t>& _shape, size_t _word_size, bool _fortran_order) :
-            shape(_shape), word_size(_word_size), fortran_order(_fortran_order)
-        {
-            num_vals = 1;
-            for(size_t i = 0;i < shape.size();i++) num_vals *= shape[i];
-            data_holder = std::shared_ptr<std::vector<char>>(
-                new std::vector<char>(num_vals * word_size));
-        }
-
-        NpyArray() : shape(0), word_size(0), fortran_order(0), num_vals(0) { }
-
-        template<typename T>
-        T* data() {
-            return reinterpret_cast<T*>(&(*data_holder)[0]);
-        }
-
-        template<typename T>
-        const T* data() const {
-            return reinterpret_cast<T*>(&(*data_holder)[0]);
-        }
-
-        template<typename T>
-        std::vector<T> as_vec() const {
-            const T* p = data<T>();
-            return std::vector<T>(p, p+num_vals);
-        }
-
-        size_t num_bytes() const {
-            return data_holder->size();
-        }
-
-        std::shared_ptr<std::vector<char>> data_holder;
-        std::vector<size_t> shape;
-        size_t word_size;
-        bool fortran_order;
-        size_t num_vals;
-    };
-   
-    using npz_t = std::map<std::string, NpyArray>; 
-
-    char BigEndianTest();
-    char map_type(const std::type_info& t);
-    template<typename T> std::vector<char> create_npy_header(const std::vector<size_t>& shape);
-    void parse_npy_header(FILE* fp,size_t& word_size, std::vector<size_t>& shape, bool& fortran_order);
-    void parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector<size_t>& shape, bool& fortran_order);
-    void parse_zip_footer(FILE* fp, uint16_t& nrecs, size_t& global_header_size, size_t& global_header_offset);
-    npz_t npz_load(std::string fname);
-    NpyArray npz_load(std::string fname, std::string varname);
-    NpyArray npy_load(std::string fname);
-
-    template<typename T> std::vector<char>& operator+=(std::vector<char>& lhs, const T rhs) {
-        //write in little endian
-        for(size_t byte = 0; byte < sizeof(T); byte++) {
-            char val = *((char*)&rhs+byte); 
-            lhs.push_back(val);
-        }
-        return lhs;
-    }
-
-    template<> std::vector<char>& operator+=(std::vector<char>& lhs, const std::string rhs);
-    template<> std::vector<char>& operator+=(std::vector<char>& lhs, const char* rhs);
-
-
-    template<typename T> void npy_save(std::string fname, const T* data, const std::vector<size_t> shape, std::string mode = "w") {
-        FILE* fp = NULL;
-        std::vector<size_t> true_data_shape; //if appending, the shape of existing + new data
-
-        if(mode == "a") fp = fopen(fname.c_str(),"r+b");
-
-        if(fp) {
-            //file exists. we need to append to it. read the header, modify the array size
-            size_t word_size;
-            bool fortran_order;
-            parse_npy_header(fp,word_size,true_data_shape,fortran_order);
-            assert(!fortran_order);
-
-            if(word_size != sizeof(T)) {
-                std::cout<<"libnpy error: "<<fname<<" has word size "<<word_size<<" but npy_save appending data sized "<<sizeof(T)<<"\n";
-                assert( word_size == sizeof(T) );
-            }
-            if(true_data_shape.size() != shape.size()) {
-                std::cout<<"libnpy error: npy_save attempting to append misdimensioned data to "<<fname<<"\n";
-                assert(true_data_shape.size() != shape.size());
-            }
-
-            for(size_t i = 1; i < shape.size(); i++) {
-                if(shape[i] != true_data_shape[i]) {
-                    std::cout<<"libnpy error: npy_save attempting to append misshaped data to "<<fname<<"\n";
-                    assert(shape[i] == true_data_shape[i]);
-                }
-            }
-            true_data_shape[0] += shape[0];
-        }
-        else {
-            fp = fopen(fname.c_str(),"wb");
-            true_data_shape = shape;
-        }
-
-        std::vector<char> header = create_npy_header<T>(true_data_shape);
-        size_t nels = std::accumulate(shape.begin(),shape.end(),1,std::multiplies<size_t>());
-
-        fseek(fp,0,SEEK_SET);
-        fwrite(&header[0],sizeof(char),header.size(),fp);
-        fseek(fp,0,SEEK_END);
-        fwrite(data,sizeof(T),nels,fp);
-        fclose(fp);
-    }
-
-    template<typename T> void npz_save(std::string zipname, std::string fname, const T* data, const std::vector<size_t>& shape, std::string mode = "w")
-    {
-        //first, append a .npy to the fname
-        fname += ".npy";
-
-        //now, on with the show
-        FILE* fp = NULL;
-        uint16_t nrecs = 0;
-        size_t global_header_offset = 0;
-        std::vector<char> global_header;
-
-        if(mode == "a") fp = fopen(zipname.c_str(),"r+b");
-
-        if(fp) {
-            //zip file exists. we need to add a new npy file to it.
-            //first read the footer. this gives us the offset and size of the global header
-            //then read and store the global header.
-            //below, we will write the the new data at the start of the global header then append the global header and footer below it
-            size_t global_header_size;
-            parse_zip_footer(fp,nrecs,global_header_size,global_header_offset);
-            fseek(fp,global_header_offset,SEEK_SET);
-            global_header.resize(global_header_size);
-            size_t res = fread(&global_header[0],sizeof(char),global_header_size,fp);
-            if(res != global_header_size){
-                throw std::runtime_error("npz_save: header read error while adding to existing zip");
-            }
-            fseek(fp,global_header_offset,SEEK_SET);
-        }
-        else {
-            fp = fopen(zipname.c_str(),"wb");
-        }
-
-        std::vector<char> npy_header = create_npy_header<T>(shape);
-
-        size_t nels = std::accumulate(shape.begin(),shape.end(),1,std::multiplies<size_t>());
-        size_t nbytes = nels*sizeof(T) + npy_header.size();
-
-        //get the CRC of the data to be added
-        uint32_t crc = crc32(0L,(uint8_t*)&npy_header[0],npy_header.size());
-        crc = crc32(crc,(uint8_t*)data,nels*sizeof(T));
-
-        //build the local header
-        std::vector<char> local_header;
-        local_header += "PK"; //first part of sig
-        local_header += (uint16_t) 0x0403; //second part of sig
-        local_header += (uint16_t) 20; //min version to extract
-        local_header += (uint16_t) 0; //general purpose bit flag
-        local_header += (uint16_t) 0; //compression method
-        local_header += (uint16_t) 0; //file last mod time
-        local_header += (uint16_t) 0;     //file last mod date
-        local_header += (uint32_t) crc; //crc
-        local_header += (uint32_t) nbytes; //compressed size
-        local_header += (uint32_t) nbytes; //uncompressed size
-        local_header += (uint16_t) fname.size(); //fname length
-        local_header += (uint16_t) 0; //extra field length
-        local_header += fname;
-
-        //build global header
-        global_header += "PK"; //first part of sig
-        global_header += (uint16_t) 0x0201; //second part of sig
-        global_header += (uint16_t) 20; //version made by
-        global_header.insert(global_header.end(),local_header.begin()+4,local_header.begin()+30);
-        global_header += (uint16_t) 0; //file comment length
-        global_header += (uint16_t) 0; //disk number where file starts
-        global_header += (uint16_t) 0; //internal file attributes
-        global_header += (uint32_t) 0; //external file attributes
-        global_header += (uint32_t) global_header_offset; //relative offset of local file header, since it begins where the global header used to begin
-        global_header += fname;
-
-        //build footer
-        std::vector<char> footer;
-        footer += "PK"; //first part of sig
-        footer += (uint16_t) 0x0605; //second part of sig
-        footer += (uint16_t) 0; //number of this disk
-        footer += (uint16_t) 0; //disk where footer starts
-        footer += (uint16_t) (nrecs+1); //number of records on this disk
-        footer += (uint16_t) (nrecs+1); //total number of records
-        footer += (uint32_t) global_header.size(); //nbytes of global headers
-        footer += (uint32_t) (global_header_offset + nbytes + local_header.size()); //offset of start of global headers, since global header now starts after newly written array
-        footer += (uint16_t) 0; //zip file comment length
-
-        //write everything
-        fwrite(&local_header[0],sizeof(char),local_header.size(),fp);
-        fwrite(&npy_header[0],sizeof(char),npy_header.size(),fp);
-        fwrite(data,sizeof(T),nels,fp);
-        fwrite(&global_header[0],sizeof(char),global_header.size(),fp);
-        fwrite(&footer[0],sizeof(char),footer.size(),fp);
-        fclose(fp);
-    }
-
-    template<typename T> void npy_save(std::string fname, const std::vector<T> data, std::string mode = "w") {
-        std::vector<size_t> shape;
-        shape.push_back(data.size());
-        npy_save(fname, &data[0], shape, mode);
-    }
-
-    template<typename T> void npz_save(std::string zipname, std::string fname, const std::vector<T> data, std::string mode = "w") {
-        std::vector<size_t> shape;
-        shape.push_back(data.size());
-        npz_save(zipname, fname, &data[0], shape, mode);
-    }
-
-    template<typename T> std::vector<char> create_npy_header(const std::vector<size_t>& shape) {  
-
-        std::vector<char> dict;
-        dict += "{'descr': '";
-        dict += BigEndianTest();
-        dict += map_type(typeid(T));
-        dict += std::to_string(sizeof(T));
-        dict += "', 'fortran_order': False, 'shape': (";
-        dict += std::to_string(shape[0]);
-        for(size_t i = 1;i < shape.size();i++) {
-            dict += ", ";
-            dict += std::to_string(shape[i]);
-        }
-        if(shape.size() == 1) dict += ",";
-        dict += "), }";
-        //pad with spaces so that preamble+dict is modulo 16 bytes. preamble is 10 bytes. dict needs to end with \n
-        int remainder = 16 - (10 + dict.size()) % 16;
-        dict.insert(dict.end(),remainder,' ');
-        dict.back() = '\n';
-
-        std::vector<char> header;
-        header += (char) 0x93;
-        header += "NUMPY";
-        header += (char) 0x01; //major version of numpy format
-        header += (char) 0x00; //minor version of numpy format
-        header += (uint16_t) dict.size();
-        header.insert(header.end(),dict.begin(),dict.end());
-
-        return header;
-    }
-
-
-}

spdlog @ 58875bdc

-Subproject commit 58875bdcd790dd2f5cb8d95ef1da7970de0a4530

Tools/plot_tracks.py

-#!/usr/bin/env python3
-
-import os
-import sys, getopt
-import re
-
-# (c) Copyright 2020 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.
-
-# with this script you can plot an animation of output of TrackWriter
-
-tracks_dat="$1"
-if [ -z "$tracks_dat" ]; then
-  echo "usage: $0 <tracks.dat> [output.gif]" >&2
-  exit 1
-fi
-
-directory=$(dirname "$tracks_dat")
-hadrons_dat="$directory/hadrons.dat"
-muons_dat="$directory/muons.dat"
-electrons_dat="$directory/electrons.dat"
-gammas_dat="$directory/gammas.dat"
-
-#if [ "$muons_dat" -ot "$tracks_dat" -o "$hadrons_dat" -ot "$muons_dat" ]; then
-  echo "splitting $tracks_dat into $muons_dat and $hadrons_dat."
-  cat "$tracks_dat" | egrep '^\s+-*13\s' > "$muons_dat"
-  cat "$tracks_dat" | egrep '^\s+-*11\s' > "$electrons_dat"
-  cat "$tracks_dat" | egrep '^\s+-*22\s' > "$gammas_dat"
-  cat "$tracks_dat" | egrep -v '^\s+-*13\s' | egrep -v '^\s+-*11\s' | egrep -v '^\s+-*22\s' > "$hadrons_dat"
-#fi
-
-output="$2"
-if [ -z "$output" ]; then
-  output="$tracks_dat.gif"
-fi
-
-echo "creating $output..."
-
-cat <<EOF | gnuplot
-set term gif animate size 900,900
-set output "$output"
-
-#set zrange [0:40e3]
-#set xrange [-10:10]
-#set yrange [-10:10]
-set xlabel "x / m"
-set ylabel "y / m"
-set zlabel "z / m"
-set title "CORSIKA 8 preliminary"
-
-do for [t=0:359:1] {
-# for separate file per angle:
-#	set output sprintf("%03d_$output", t)
-
-	set view 80, t
-        splot "$gammas_dat" u 3:4:5:6:7:8 w vectors nohead lt rgb "orange" t "", "$electrons_dat" u 3:4:5:6:7:8 w vectors nohead lt rgb "blue" t "", "$muons_dat" u 3:4:5:6:7:8 w vectors nohead lt rgb "red" t "", "$hadrons_dat" u 3:4:5:6:7:8 w vectors nohead  lc rgb "black" t ""
-}
-EOF
-
-exit $?

MCNET_GUIDELINES → USING_COLLABORATING.md

-These guidelines are copied from: http://www.montecarlonet.org/
----------------------------------------------------------------
+# Usage and collaboration agreement
 
+The CORSIKA 8 project very much welcomes all collaboration and
+contributions. The aim of the project is to create a
+scientific software framework as a fundamental tool for research.
 
-               MCNET GUIDELINES
+The project consists of the contributions from the scientific
+community and individuals in a best effort to deliver the best
+possible performance and physics output. 
 
-      for Event Generator Authors and Users
 
+## The software license of the CORSIKA project
 
-    PREAMBLE
+CORSIKA 8 is by default released under the BSD 3-Clause License, as copied in full in the file
+[LICENSE](LICENSE). Each source file of the CORSIKA project contains a
+short statement of the copyright and this license. Each binary or
+source code release of CORSIKA contains the file LICENSE.
 
-This generator has been developed as part of an academic research
-project and is the result of many years of work by the authors.
-Proper academic recognition is a requirement for its continued
-development.
+The code, documentation and content in the folder [./externals](./externals)
+is not integral part of the CORSIKA project and can be based on, or
+include, other licenses, which must be compatible with the CORSIKA 8 license. 
 
-The components of the program have been developed to work together
-as a coherent physics framework. We believe that the creation of
-separately maintained forks or piecewise distribution of individual
-parts would diminish their scientific value.
+The folder [./modules](./modules) contains the code of several
+external physics models for your convenience. They each come with
+their own license which we ask you to honor.  Please also make sure to cite the
+adequate reference papers when using their models in scientific work
+and publications.
 
-The authors are convinced that software development in a scientific
-context requires full availability of all source code, to further
-progress and to allow local modifications to meet the specific
-requirements of the individual user.
+Of course, we have the authors' consent to
+distribute their code together with CORSIKA 8. 
 
-Therefore we have decided to release this program under the GNU
-General Public License (GPL) version 2 (with the option to instead
-follow the terms and conditions of any later version of GPL). This
-ensures that the source code will be available to you and grants you
-the freedom to use and modify the program. You can redistribute your
-modified versions as long as you retain the GPL and respect existing
-copyright notices (see the file 'COPYING' for details).
+Check the content of these folders carefully for details and additional 
+license information. It depends on the configuration of
+the build system to what extent this code is used to build CORSIKA.
 
-By using the GPL, we entrust you with considerable freedom and expect
-you to use it wisely, since the GPL does not address the issues in
-the first two paragraphs. To remedy this shortcoming, we have
-formulated the following guidelines relevant for the distribution
-and usage of event generator software in an academic setting.
 
+## Contributing
 
-    GUIDELINES
+If you want to contribute, you need to read
+[the contributing GUIDELINES](CONTRIBUTING.md) and comply with these rules, or help to
+improve them. 
+
+
+## General guidelines
+
+We reproduce below some guidelines copied from  http://www.montecarlonet.org/ that we also ask you to follow in the spirit of academic collaboration.
 
 1) The integrity of the program should be respected.
     -------------------------------------------------
@@ -90,23 +94,7 @@ and usage of event generator software in an academic setting.
     modifications should be spelled out.
 
 
-    POSTSCRIPT
-
-The copyright license of the software is the GPL v2 alone, therefore
-the above guidelines are not legally binding. However, we reserve the
-right to criticize offenders. The guidelines should always be combined
-with common sense, for interpretation and for issues not covered.
-Enquiries regarding the guidelines and related issues are encouraged
-and should be directed to the authors of the program.
-
-Please note that the program, including all its code and documentation,
-is intended for academic use and is delivered "as is" to be used at 
-your own risk, without any guarantees.
-
-----------------------------------------------------------------------
-
-These guidelines were edited by Nils Lavesson and David Grellscheid
+(These guidelines were originally edited by Nils Lavesson and David Grellscheid
 for the MCnet collaboration, which has approved and agreed to respect
 them.  MCnet is a Marie Curie Research Training Network funded under
-Framework Programme 6 contract MRTN-CT-2006-035606.
-
+Framework Programme 6 contract MRTN-CT-2006-035606.)

applications/CMakeLists.txt

+SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+
+add_executable (c8_air_shower c8_air_shower.cpp)
+target_link_libraries (c8_air_shower CORSIKA8)
+
+if(WITH_FLUKA)
+  message("compiling c8_air_shower.cpp with FLUKA")
+  target_compile_definitions(c8_air_shower PRIVATE WITH_FLUKA)
+else()
+  message("compiling c8_air_shower.cpp with UrQMD")
+endif()
+
+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.

applications/c8_air_shower.cpp

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/* clang-format off */
+// InteractionCounter used boost/histogram, which
+// fails if boost/type_traits have been included before. Thus, we have
+// to include it first...
+#include <corsika/framework/process/InteractionCounter.hpp>
+/* clang-format on */
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/process/DynamicInteractionProcess.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/random/PowerLawDistribution.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+
+#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/InteractionWriter.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/output/OutputManager.hpp>
+
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/GeomagneticModel.hpp>
+#include <corsika/media/GladstoneDaleRefractiveIndex.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+#include <corsika/modules/BetheBlochPDG.hpp>
+#include <corsika/modules/Epos.hpp>
+#include <corsika/modules/LongitudinalProfile.hpp>
+#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
+#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/Pythia8.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/Sophia.hpp>
+#include <corsika/modules/StackInspector.hpp>
+#include <corsika/modules/thinning/EMThinning.hpp>
+// for ICRC2023
+#ifdef WITH_FLUKA
+#include <corsika/modules/FLUKA.hpp>
+#else
+#include <corsika/modules/UrQMD.hpp>
+#endif
+#include <corsika/modules/TAUOLA.hpp>
+
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/RadioProcess.hpp>
+#include <corsika/modules/radio/ZHS.hpp>
+#include <corsika/modules/radio/observers/Observer.hpp>
+#include <corsika/modules/radio/observers/TimeDomainObserver.hpp>
+#include <corsika/modules/radio/detectors/ObserverCollection.hpp>
+#include <corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/setup/SetupC7trackedParticles.hpp>
+
+#include <boost/filesystem.hpp>
+
+#include <CLI/App.hpp>
+#include <CLI/Config.hpp>
+#include <CLI/Formatter.hpp>
+
+#include <cstdlib>
+#include <iomanip>
+#include <limits>
+#include <string>
+
+using namespace corsika;
+using namespace std;
+
+using EnvironmentInterface =
+    IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+using EnvType = Environment<EnvironmentInterface>;
+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");
+  RNGManager<>::getInstance().registerRandomStream("qgsjet");
+  RNGManager<>::getInstance().registerRandomStream("sibyll");
+  RNGManager<>::getInstance().registerRandomStream("sophia");
+  RNGManager<>::getInstance().registerRandomStream("epos");
+  RNGManager<>::getInstance().registerRandomStream("pythia");
+  RNGManager<>::getInstance().registerRandomStream("urqmd");
+  RNGManager<>::getInstance().registerRandomStream("fluka");
+  RNGManager<>::getInstance().registerRandomStream("proposal");
+  RNGManager<>::getInstance().registerRandomStream("thinning");
+  RNGManager<>::getInstance().registerRandomStream("primary_particle");
+  if (seed == 0) {
+    std::random_device rd;
+    seed = rd();
+    CORSIKA_LOG_INFO("random seed (auto) {}", seed);
+  } else {
+    CORSIKA_LOG_INFO("random seed {}", seed);
+  }
+  RNGManager<>::getInstance().setSeed(seed);
+  return seed;
+}
+
+template <typename T>
+using MyExtraEnv =
+    GladstoneDaleRefractiveIndex<MediumPropertyModel<UniformMagneticField<T>>>;
+
+int main(int argc, char** argv) {
+
+  // the main command line description
+  CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."};
+
+  CORSIKA_LOG_INFO(
+      "Please cite the following papers when using CORSIKA 8:\n"
+      " - \"Towards a Next Generation of CORSIKA: A Framework for the Simulation of "
+      "Particle Cascades in Astroparticle Physics\", Comput. Softw. Big Sci. 3 (2019) "
+      "2, https://doi.org/10.1007/s41781-018-0013-0\n"
+      " - \"Simulating radio emission from particle cascades with CORSIKA 8\", "
+      "Astropart. Phys. 166 (2025) 103072, "
+      "https://doi.org/10.1016/j.astropartphys.2024.103072");
+
+  //////// Primary options ////////
+
+  // some options that we want to fill in
+  int A, Z, nevent = 0;
+  std::vector<double> cli_energy_range;
+
+  // the following section adds the options to the parser
+
+  // we start by definining a sub-group for the primary ID
+  auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary")
+                   ->check(CLI::Range(0, 26))
+                   ->group("Primary");
+  auto opt_A = app.add_option("-A", A, "Atomic mass number for primary")
+                   ->needs(opt_Z)
+                   ->check(CLI::Range(1, 58))
+                   ->group("Primary");
+  app.add_option("-p,--pdg",
+                 "PDG code for primary (p=2212, gamma=22, e-=11, nu_e=12, mu-=13, "
+                 "nu_mu=14, tau=15, nu_tau=16).")
+      ->excludes(opt_A)
+      ->excludes(opt_Z)
+      ->group("Primary");
+  app.add_option("-E,--energy", "Primary energy in GeV")->default_val(0);
+  app.add_option("--energy_range", cli_energy_range,
+                 "Low and high values that define the range of the primary energy in GeV")
+      ->expected(2)
+      ->check(CLI::PositiveNumber)
+      ->group("Primary");
+  app.add_option("--eslope", "Spectral index for sampling energies, dN/dE = E^eSlope")
+      ->default_val(-1.0)
+      ->group("Primary");
+  app.add_option("-z,--zenith", "Primary zenith angle (deg)")
+      ->default_val(0.)
+      ->check(CLI::Range(0., 90.))
+      ->group("Primary");
+  app.add_option("-a,--azimuth", "Primary azimuth angle (deg)")
+      ->default_val(0.)
+      ->check(CLI::Range(0., 360.))
+      ->group("Primary");
+
+  //////// Config options ////////
+
+  app.add_option("--emcut",
+                 "Min. kin. energy of photons, electrons and "
+                 "positrons in tracking (GeV)")
+      ->default_val(0.5e-3)
+      ->check(CLI::Range(0.000001, 1.e13))
+      ->group("Config");
+  app.add_option("--hadcut", "Min. kin. energy of hadrons in tracking (GeV)")
+      ->default_val(0.3)
+      ->check(CLI::Range(0.02, 1.e13))
+      ->group("Config");
+  app.add_option("--mucut", "Min. kin. energy of muons in tracking (GeV)")
+      ->default_val(0.3)
+      ->check(CLI::Range(0.000001, 1.e13))
+      ->group("Config");
+  app.add_option("--taucut", "Min. kin. energy of tau leptons in tracking (GeV)")
+      ->default_val(0.3)
+      ->check(CLI::Range(0.000001, 1.e13))
+      ->group("Config");
+  app.add_option("--max-deflection-angle",
+                 "maximal deflection angle in tracking in radians")
+      ->default_val(0.2)
+      ->check(CLI::Range(1.e-8, 1.))
+      ->group("Config");
+  bool track_neutrinos = false;
+  app.add_flag("--track-neutrinos", track_neutrinos, "switch on tracking of neutrinos")
+      ->group("Config");
+
+  //////// Misc options ////////
+
+  app.add_option("--neutrino-interaction-type",
+                 "charged (CC) or neutral current (NC) or both")
+      ->default_val("both")
+      ->check(CLI::IsMember({"neutral", "NC", "charged", "CC", "both"}))
+      ->group("Misc.");
+  app.add_option("--observation-level",
+                 "Height above earth radius of the observation level (in m)")
+      ->default_val(0.)
+      ->check(CLI::Range(-1.e3, 1.e5))
+      ->group("Config");
+  app.add_option("--injection-height",
+                 "Height above earth radius of the injection point (in m)")
+      ->default_val(112.75e3)
+      ->check(CLI::Range(-1.e3, 1.e6))
+      ->group("Config");
+  app.add_option("-N,--nevent", nevent, "The number of events/showers to run.")
+      ->default_val(1)
+      ->check(CLI::PositiveNumber)
+      ->group("Library/Output");
+  app.add_option("-f,--filename", "Filename for output library.")
+      ->required()
+      ->default_val("corsika_library")
+      ->check(CLI::NonexistentPath)
+      ->group("Library/Output");
+  bool compressOutput = false;
+  app.add_flag("--compress", compressOutput, "Compress the output directory to a tarball")
+      ->group("Library/Output");
+  app.add_option("-s,--seed", "The random number seed.")
+      ->default_val(0)
+      ->check(CLI::NonNegativeNumber)
+      ->group("Misc.");
+  bool force_interaction = false;
+  app.add_flag("--force-interaction", force_interaction,
+               "Force the location of the first interaction.")
+      ->group("Misc.");
+  bool force_decay = false;
+  app.add_flag("--force-decay", force_decay, "Force the primary to immediately decay")
+      ->group("Misc.");
+  bool disable_interaction_hists = false;
+  app.add_flag("--disable-interaction-histograms", disable_interaction_hists,
+               "Store interaction histograms")
+      ->group("Misc.");
+  app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.")
+      ->default_val("info")
+      ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
+      ->group("Misc.");
+  app.add_option("-M,--hadronModel", "High-energy hadronic interaction model")
+      ->default_val("SIBYLL-2.3d")
+      ->check(CLI::IsMember({"SIBYLL-2.3d", "QGSJet-II.04", "EPOS-LHC", "Pythia8"}))
+      ->group("Misc.");
+  app.add_option("-T,--hadronModelTransitionEnergy",
+                 "Transition between high-/low-energy hadronic interaction "
+                 "model in GeV")
+      ->default_val(std::pow(10, 1.9)) // 79.4 GeV
+      ->check(CLI::NonNegativeNumber)
+      ->group("Misc.");
+
+  //////// Thinning options ////////
+
+  app.add_option("--emthin",
+                 "fraction of primary energy at which thinning of EM particles starts")
+      ->default_val(1.e-6)
+      ->check(CLI::Range(0., 1.))
+      ->group("Thinning");
+  app.add_option("--max-weight",
+                 "maximum weight for thinning of EM particles (0 to select Kobal's "
+                 "optimum times 0.5)")
+      ->default_val(0)
+      ->check(CLI::NonNegativeNumber)
+      ->group("Thinning");
+  bool multithin = false;
+  app.add_flag("--multithin", multithin, "keep thinned particles (with weight=0)")
+      ->group("Thinning");
+  app.add_option("--ring", "concentric ring of star shape pattern of observers")
+      ->default_val(0)
+      ->check(CLI::Range(0, 20))
+      ->group("Radio");
+  // parse the command line options into the variables
+  CLI11_PARSE(app, argc, argv);
+
+  if (app.count("--verbosity")) {
+    auto const loglevel = app["--verbosity"]->as<std::string>();
+    if (loglevel == "warn") {
+      logging::set_level(logging::level::warn);
+    } else if (loglevel == "info") {
+      logging::set_level(logging::level::info);
+    } else if (loglevel == "debug") {
+      logging::set_level(logging::level::debug);
+    } else if (loglevel == "trace") {
+#ifndef _C8_DEBUG_
+      CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
+      return 1;
+#endif
+      logging::set_level(logging::level::trace);
+    }
+  }
+
+  // check that we got either PDG or A/Z
+  // 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.");
+      return 1;
+    }
+  }
+
+  // initialize random number sequence(s)
+  auto seed = registerRandomStreams(app["--seed"]->as<long>());
+
+  /* === START: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
+  EnvType env;
+  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+  Point const center{rootCS, 0_m, 0_m, 0_m};
+  Point const surface_{rootCS, 0_m, 0_m, constants::EarthRadius::Mean};
+  GeomagneticModel wmm(center, corsika_data("GeoMag/WMM.COF"));
+
+  // build an atmosphere with Keilhauer's parametrization of the
+  // US standard atmosphere into `env`
+  create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+      env, AtmosphereId::USStdBK, center, 1.000327, surface_, Medium::AirDry1Atm,
+      MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
+
+  /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
+
+  /* === START: CONSTRUCT PRIMARY PARTICLE === */
+
+  // parse the primary ID as a PDG or A/Z code
+  Code beamCode;
+
+  // check if we want to use a PDG code instead
+  if (app.count("--pdg") > 0) {
+    beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
+  } else {
+    // check manually for proton and neutrons
+    if ((A == 1) && (Z == 1))
+      beamCode = Code::Proton;
+    else if ((A == 1) && (Z == 0))
+      beamCode = Code::Neutron;
+    else
+      beamCode = get_nucleus_code(A, Z);
+  }
+
+  HEPEnergyType eMin = 0_GeV;
+  HEPEnergyType eMax = 0_GeV;
+  // check the particle energy parameters
+  if (app["--energy"]->as<double>() > 0.0) {
+    eMin = app["--energy"]->as<double>() * 1_GeV;
+    eMax = app["--energy"]->as<double>() * 1_GeV;
+  } else if (cli_energy_range.size()) {
+    if (cli_energy_range[0] > cli_energy_range[1]) {
+      CORSIKA_LOG_WARN(
+          "Energy range lower bound is greater than upper bound. swapping...");
+      eMin = cli_energy_range[1] * 1_GeV;
+      eMax = cli_energy_range[0] * 1_GeV;
+    } else {
+      eMin = cli_energy_range[0] * 1_GeV;
+      eMax = cli_energy_range[1] * 1_GeV;
+    }
+  } else {
+    CORSIKA_LOG_CRITICAL(
+        "Must set either the (--energy) flag or the (--energy_range) flag to "
+        "positive value(s)");
+    return 0;
+  }
+
+  // direction of the shower in (theta, phi) space
+  auto const thetaRad = app["--zenith"]->as<double>() / 180. * M_PI;
+  auto const phiRad = app["--azimuth"]->as<double>() / 180. * M_PI;
+
+  auto const [nx, ny, nz] = std::make_tuple(sin(thetaRad) * cos(phiRad),
+                                            sin(thetaRad) * sin(phiRad), -cos(thetaRad));
+  auto propDir = DirectionVector(rootCS, {nx, ny, nz});
+  /* === END: CONSTRUCT PRIMARY PARTICLE === */
+
+  /* === START: CONSTRUCT GEOMETRY === */
+  auto const observationHeight =
+      app["--observation-level"]->as<double>() * 1_m + constants::EarthRadius::Mean;
+  auto const injectionHeight =
+      app["--injection-height"]->as<double>() * 1_m + constants::EarthRadius::Mean;
+  auto const t = -observationHeight * cos(thetaRad) +
+                 sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+                      static_pow<2>(injectionHeight));
+  Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+  Point const injectionPos =
+      showerCore + DirectionVector{rootCS,
+                                   {-sin(thetaRad) * cos(phiRad),
+                                    -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
+                       t;
+
+  // we make the axis much longer than the inj-core distance since the
+  // profile will go beyond the core, depending on zenith angle
+  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
+  /* === END: CONSTRUCT GEOMETRY === */
+
+  std::stringstream args;
+  for (int i = 0; i < argc; ++i) { args << argv[i] << " "; }
+  // create the output manager that we then register outputs with
+  OutputManager output(app["--filename"]->as<std::string>(), seed, args.str(),
+                       compressOutput);
+
+  // register energy losses as output
+  EnergyLossWriter dEdX{showerAxis, dX};
+  output.add("energyloss", dEdX);
+
+  DynamicInteractionProcess<StackType> heModel;
+
+  auto const all_elements = corsika::get_all_elements_in_universe(env);
+  // have SIBYLL always for PROPOSAL photo-hadronic interactions
+  auto sibyll = std::make_shared<corsika::sibyll::Interaction>(
+      all_elements, corsika::setup::C7trackedParticles);
+
+  if (auto const modelStr = app["--hadronModel"]->as<std::string>();
+      modelStr == "SIBYLL-2.3d") {
+    heModel = DynamicInteractionProcess<StackType>{sibyll};
+  } else if (modelStr == "QGSJet-II.04") {
+    heModel = DynamicInteractionProcess<StackType>{
+        std::make_shared<corsika::qgsjetII::Interaction>()};
+  } else if (modelStr == "EPOS-LHC") {
+    heModel = DynamicInteractionProcess<StackType>{
+        std::make_shared<corsika::epos::Interaction>(corsika::setup::C7trackedParticles)};
+  } else if (modelStr == "Pythia8") {
+    heModel = DynamicInteractionProcess<StackType>{
+        std::make_shared<corsika::pythia8::Interaction>(
+            corsika::setup::C7trackedParticles)};
+  } else {
+    CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr);
+    return EXIT_FAILURE;
+  }
+
+  InteractionCounter heCounted{heModel};
+
+  corsika::pythia8::Decay decayPythia;
+  // tau decay via TAUOLA (hard coded to left handed)
+  corsika::tauola::Decay decayTauola(corsika::tauola::Helicity::LeftHanded);
+
+  struct IsTauSwitch {
+    bool operator()(const Particle& p) const {
+      return (p.getPID() == Code::TauMinus || p.getPID() == Code::TauPlus);
+    }
+  };
+
+  auto decaySequence = make_select(IsTauSwitch(), decayTauola, decayPythia);
+
+  // neutrino interactions with pythia (options are: NC, CC)
+  bool NC = false;
+  bool CC = false;
+  if (auto const nuIntStr = app["--neutrino-interaction-type"]->as<std::string>();
+      nuIntStr == "neutral" || nuIntStr == "NC") {
+    NC = true;
+    CC = false;
+  } else if (nuIntStr == "charged" || nuIntStr == "CC") {
+    NC = false;
+    CC = true;
+  } else if (nuIntStr == "both") {
+    NC = true;
+    CC = true;
+  }
+  corsika::pythia8::NeutrinoInteraction neutrinoPrimaryPythia(
+      corsika::setup::C7trackedParticles, NC, CC);
+
+  // hadronic photon interactions in resonance region
+  corsika::sophia::InteractionModel sophia;
+
+  HEPEnergyType const emcut = 1_GeV * app["--emcut"]->as<double>();
+  HEPEnergyType const hadcut = 1_GeV * app["--hadcut"]->as<double>();
+  HEPEnergyType const mucut = 1_GeV * app["--mucut"]->as<double>();
+  HEPEnergyType const taucut = 1_GeV * app["--taucut"]->as<double>();
+  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut, taucut,
+                                             !track_neutrinos, dEdX);
+
+  // tell proposal that we are interested in all energy losses above the particle cut
+  auto const prod_threshold = std::min({emcut, hadcut, mucut, taucut});
+  set_energy_production_threshold(Code::Electron, prod_threshold);
+  set_energy_production_threshold(Code::Positron, prod_threshold);
+  set_energy_production_threshold(Code::Photon, prod_threshold);
+  set_energy_production_threshold(Code::MuMinus, prod_threshold);
+  set_energy_production_threshold(Code::MuPlus, prod_threshold);
+  set_energy_production_threshold(Code::TauMinus, prod_threshold);
+  set_energy_production_threshold(Code::TauPlus, prod_threshold);
+
+  // energy threshold for high energy hadronic model. Affects LE/HE switch for
+  // hadron interactions and the hadronic photon model in proposal
+  HEPEnergyType const heHadronModelThreshold =
+      1_GeV * app["--hadronModelTransitionEnergy"]->as<double>();
+
+  corsika::proposal::Interaction emCascade(
+      env, sophia, sibyll->getHadronInteractionModel(), heHadronModelThreshold);
+
+  // use BetheBlochPDG for hadronic continuous losses, and proposal otherwise
+  corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuousProposal(
+      env, dEdX);
+  BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuousBethe{dEdX};
+  struct EMHadronSwitch {
+    EMHadronSwitch() = default;
+    bool operator()(const Particle& p) const { return is_hadron(p.getPID()); }
+  };
+  auto emContinuous =
+      make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
+
+  LongitudinalWriter profile{showerAxis, dX};
+  output.add("profile", profile);
+  LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
+
+// for ICRC2023
+#ifdef WITH_FLUKA
+  corsika::fluka::Interaction leIntModel{all_elements};
+#else
+  corsika::urqmd::UrQMD leIntModel{};
+#endif
+  InteractionCounter leIntCounted{leIntModel};
+
+  // assemble all processes into an ordered process list
+  struct EnergySwitch {
+    HEPEnergyType cutE_;
+    EnergySwitch(HEPEnergyType cutE)
+        : cutE_(cutE) {}
+    bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
+  };
+  auto hadronSequence =
+      make_select(EnergySwitch(heHadronModelThreshold), leIntCounted, heCounted);
+
+  // observation plane
+  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+  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<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
+  output.add("primary", primaryWriter);
+
+  int ring_number{app["--ring"]->as<int>()};
+  auto const radius_{ring_number * 25_m};
+  const int rr_ = static_cast<int>(radius_ / 1_m);
+
+  // Radio observers and relevant information
+  // the observer time variables
+  const TimeType duration_{4e-7_s};
+  const InverseTimeType sampleRate_{1e+9_Hz};
+
+  // the observer collection for CoREAS and ZHS
+  ObserverCollection<TimeDomainObserver> detectorCoREAS;
+  ObserverCollection<TimeDomainObserver> detectorZHS;
+
+  auto const showerCoreX_{showerCore.getCoordinates().getX()};
+  auto const showerCoreY_{showerCore.getCoordinates().getY()};
+  auto const injectionPosX_{injectionPos.getCoordinates().getX()};
+  auto const injectionPosY_{injectionPos.getCoordinates().getY()};
+  auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
+  auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
+
+  if (ring_number != 0) {
+    // setup CoREAS observers - 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 const point_1{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_1),
+                               showerCoreY_ + radius_ * sin(phiRad_1),
+                               constants::EarthRadius::Mean)};
+      std::cout << "Observer point CoREAS: " << point_1 << std::endl;
+      auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
+      std::string name_1 = "CoREAS_R=" + std::to_string(rr_) +
+                           "_m--Phi=" + std::to_string(phi_1) + "degrees";
+      TimeDomainObserver observer_1(name_1, point_1, rootCS, triggertime_1, duration_,
+                                    sampleRate_, triggertime_1);
+      detectorCoREAS.addObserver(observer_1);
+    }
+
+    // setup ZHS observers - use the for loop for star shape pattern
+    for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
+      auto phiRad_ = phi_ / 180. * M_PI;
+      auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
+                              showerCoreY_ + radius_ * sin(phiRad_),
+                              constants::EarthRadius::Mean)};
+      std::cout << "Observer point ZHS: " << point_ << std::endl;
+      auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
+      std::string name_ =
+          "ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
+      TimeDomainObserver observer_2(name_, point_, rootCS, triggertime_, duration_,
+                                    sampleRate_, triggertime_);
+      detectorZHS.addObserver(observer_2);
+    }
+  }
+  LengthType const step = 1_m;
+  auto TP =
+      make_tabulated_flat_atmosphere_radio_propagator(env, injectionPos, surface_, step);
+
+  // initiate CoREAS
+  RadioProcess<decltype(detectorCoREAS), CoREAS<decltype(detectorCoREAS), decltype(TP)>,
+               decltype(TP)>
+      coreas(detectorCoREAS, TP);
+
+  // register CoREAS with the output manager
+  output.add("CoREAS", coreas);
+
+  // initiate ZHS
+  RadioProcess<decltype(detectorZHS), ZHS<decltype(detectorZHS), decltype(TP)>,
+               decltype(TP)>
+      zhs(detectorZHS, TP);
+
+  // register ZHS with the output manager
+  output.add("ZHS", zhs);
+
+  // make and register the first interaction writer
+  InteractionWriter<setup::Tracking, ParticleWriterParquet> inter_writer(
+      showerAxis, observationLevel);
+  output.add("interactions", inter_writer);
+
+  /* === END: SETUP PROCESS LIST === */
+
+  // trigger the output manager to open the library for writing
+  output.startOfLibrary();
+
+  // loop over each shower
+  for (int i_shower = 1; i_shower < nevent + 1; i_shower++) {
+
+    CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, nevent);
+
+    // randomize the primary energy
+    double const eSlope = app["--eslope"]->as<double>();
+    PowerLawDistribution<HEPEnergyType> powerLawRng(eSlope, eMin, eMax);
+    HEPEnergyType const primaryTotalEnergy =
+        (eMax == eMin) ? eMin
+                       : powerLawRng(RNGManager<>::getInstance().getRandomStream(
+                             "primary_particle"));
+
+    auto const eKin = primaryTotalEnergy - get_mass(beamCode);
+
+    // set up thinning based on primary parameters
+    double const emthinfrac = app["--emthin"]->as<double>();
+    double const maxWeight = std::invoke([&]() {
+      if (auto const wm = app["--max-weight"]->as<double>(); wm > 0)
+        return wm;
+      else
+        return 0.5 * emthinfrac * primaryTotalEnergy / 1_GeV;
+    });
+    EMThinning thinning{emthinfrac * primaryTotalEnergy, maxWeight, !multithin};
+
+    // set up the stack inspector
+    StackInspector<StackType> stackInspect(10000, false, primaryTotalEnergy);
+
+    // assemble the final process sequence
+    auto sequence =
+        make_sequence(stackInspect, neutrinoPrimaryPythia, hadronSequence, decaySequence,
+                      emCascade, emContinuous, coreas, zhs, longprof, observationLevel,
+                      inter_writer, thinning, cut);
+
+    // create the cascade object using the default stack and tracking
+    // implementation
+    TrackingType tracking(app["--max-deflection-angle"]->as<double>());
+    StackType stack;
+    Cascade EAS(env, tracking, sequence, output, stack);
+
+    // setup particle stack, and add primary particle
+    stack.clear();
+
+    // print our primary parameters all in one place
+    CORSIKA_LOG_INFO("Primary name:         {}", beamCode);
+    if (app["--pdg"]->count() > 0) {
+      CORSIKA_LOG_INFO("Primary PDG ID:       {}", app["--pdg"]->as<int>());
+    } else {
+      CORSIKA_LOG_INFO("Primary Z/A:          {}/{}", Z, A);
+    }
+    CORSIKA_LOG_INFO("Primary Total Energy: {}", primaryTotalEnergy);
+    CORSIKA_LOG_INFO("Primary Momentum:     {}",
+                     calculate_momentum(primaryTotalEnergy, get_mass(beamCode)));
+    CORSIKA_LOG_INFO("Primary Direction:    {}", propDir.getNorm());
+    CORSIKA_LOG_INFO("Point of Injection:   {}", injectionPos.getCoordinates());
+    CORSIKA_LOG_INFO("Shower Axis Length:   {}",
+                     (showerCore - injectionPos).getNorm() * 1.2);
+
+    // add the desired particle to the stack
+    auto const primaryProperties =
+        std::make_tuple(beamCode, eKin, propDir.normalized(), injectionPos, 0_ns);
+    stack.addParticle(primaryProperties);
+
+    // if we want to fix the first location of the shower
+    if (force_interaction) {
+      CORSIKA_LOG_INFO("Fixing first interaction at injection point.");
+      EAS.forceInteraction();
+    }
+
+    if (force_decay) {
+      CORSIKA_LOG_INFO("Forcing the primary to decay");
+      EAS.forceDecay();
+    }
+
+    primaryWriter.recordPrimary(primaryProperties);
+
+    // run the shower
+    EAS.run();
+
+    HEPEnergyType const Efinal =
+        dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+
+    CORSIKA_LOG_INFO(
+        "total energy budget (GeV): {} (dEdX={} ground={}), "
+        "relative difference (%): {}",
+        Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV,
+        observationLevel.getEnergyGround() / 1_GeV,
+        (Efinal / primaryTotalEnergy - 1) * 100);
+
+    if (!disable_interaction_hists) {
+      CORSIKA_LOG_INFO("Saving interaction histograms");
+      auto const hists = heCounted.getHistogram() + leIntCounted.getHistogram();
+
+      // directory for output of interaction histograms
+      string const outdir(app["--filename"]->as<std::string>() + "/interaction_hist");
+      boost::filesystem::create_directories(outdir);
+
+      string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
+      string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
+      save_hist(hists.labHist(), labHist_file, true);
+      save_hist(hists.CMSHist(), cMSHist_file, true);
+    }
+  }
+
+  // and finalize the output on disk
+  output.endOfLibrary();
+
+  return EXIT_SUCCESS;
+}

CMakeModules/FindCONEX.cmake → cmake/FindCONEX.cmake

@@ -3,9 +3,8 @@
 #
 # 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.
+# This software is distributed under the terms of the 3-clause BSD license.
+# See file LICENSE for a full version of the license.
 #
 
 # - find Conex