added and tested COAST interface

CMakeLists.txt

@@ -45,6 +45,13 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O3 -g") # -O2 would not trade speed for size, neit
 # clang produces a lot of unecessary warnings without this:
 add_compile_options("$<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>>:-Wno-nonportable-include-path>")
 
+# COAST - interface
+if (WITH_COAST)
+  message(STATUS "Compiling CORSIKA8 for the use with COAST/corsika7.")
+  add_compile_options("-fPIC")
+endif()
+
+
 # unit testing coverage, does not work yet
 #include (CodeCoverage)
 ##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')
@@ -79,3 +86,6 @@ add_subdirectory (Processes)
 add_subdirectory (Documentation)
 add_subdirectory (Main)
 add_subdirectory (Tools)
+if (WITH_COAST)
+  add_subdirectory (COAST)
+endif()

COAST/CMakeLists.txt

+
+set (
+  COAST_HEADERS
+  COASTProcess.h
+  COASTStack.h
+  ParticleConversion.h
+  )
+
+set (
+  COAST_SOURCES
+  COASTUserLib.cc
+  COASTProcess.cc
+  ParticleConversion.cc
+  )
+
+set (
+  COAST_NAMESPACE
+  corsika/coast
+  )
+
+add_library (COAST SHARED ${COAST_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (COAST ${COAST_NAMESPACE} ${COAST_HEADERS})
+
+set_target_properties (
+  COAST
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+#  PUBLIC_HEADER "${MODEL_HEADERS}"
+  )
+
+target_link_libraries (
+  COAST
+  PUBLIC
+  CORSIKAgeometry
+  CORSIKAunits
+  CORSIKAparticles
+  CORSIKAgeometry
+  CORSIKAsetup
+  # SuperStupidStack
+  )
+
+target_include_directories (
+  COAST
+  PUBLIC
+  $ENV{COAST_DIR}/include
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS COAST
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+  PUBLIC_HEADER DESTINATION include/${COAST_NAMESPACE}
+  )
+
+#install (
+#  FILES ${COAST_HEADERS} 
+#  DESTINATION include/${COAST_NAMESPACE}
+#  )

COAST/COASTProcess.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/coast/COASTProcess.h>
+#include <corsika/coast/COASTStack.h>
+
+#include <iostream>
+
+using namespace std;
+
+namespace corsika::coast {
+
+  /**
+     the init function is called during the start of corsika.
+   */
+  void COASTProcess::Init() {
+    cout << "************* Greetings from CORSIKA8 *****************" << endl;
+  }
+
+  /**
+     the docontinous function is called for each tracking step in
+     corsika. Take care: you cannot modify the particle, the track or
+     the stack from here (docontinuous) inside corisika7. In corsika8
+     you will be able to do that.
+   */
+  corsika::process::EProcessReturn COASTProcess::DoContinuous(Particle& p, Track& t,
+                                                              Stack&) {
+    using namespace corsika::units::si;
+    auto const start = t.GetPosition(0).GetCoordinates();
+    auto const delta = t.GetPosition(1).GetCoordinates() - start;
+    auto const name = corsika::particles::GetName(p.GetPID());
+    cout << "CORSIKA8: particle=" << name << ", pos=" << start
+         << " track-l=" << delta.norm() << ", track-t=" << t.GetDuration() << endl;
+    return corsika::process::EProcessReturn::eOk;
+  }
+
+} // namespace corsika::coast

COAST/COASTProcess.h

+#ifndef _include_corsika_coast_coastprocess_h_
+#define _include_corsika_coast_coastprocess_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/coast/COASTStack.h>
+
+#include <limits>
+
+typedef corsika::coast::COASTStack Stack;
+typedef corsika::coast::COASTStack::ParticleType Particle;
+typedef corsika::geometry::Trajectory<corsika::geometry::Line> Track;
+
+namespace corsika::coast {
+
+  class COASTProcess : public corsika::process::ContinuousProcess<COASTProcess> {
+
+  public:
+    void Init();
+
+    corsika::process::EProcessReturn DoContinuous(Particle&, Track&, Stack&);
+
+    corsika::units::si::LengthType MaxStepLength(Particle&, Track&) {
+      return corsika::units::si::meter * std::numeric_limits<double>::infinity();
+    }
+
+  private:
+  };
+
+} // namespace corsika::coast
+
+#endif

COAST/COASTStack.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_coaststack_h_
+#define _include_coaststack_h_
+
+#include <corsika/coast/ParticleConversion.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/stack/ParticleBase.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h> // remove
+#include <corsika/geometry/Vector.h>
+
+#include <crs/CParticle.h>
+#include <crs/CorsikaTypes.h>
+
+#include <algorithm>
+#include <vector>
+
+namespace corsika::coast {
+
+  typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
+
+  /**
+   * Example of a particle object on the stack.
+   */
+
+  template <typename StackIteratorInterface>
+  class ParticleInterface : public corsika::stack::ParticleBase<StackIteratorInterface> {
+
+    using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+    using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+  public:
+    corsika::particles::Code GetPID() const { return GetStackData().GetPID(GetIndex()); }
+
+    corsika::units::si::HEPEnergyType GetEnergy() const {
+      return GetStackData().GetEnergy(GetIndex());
+    }
+
+    MomentumVector GetMomentum(const corsika::geometry::CoordinateSystem& cs) const {
+      using namespace corsika::units::si;
+      const HEPEnergyType mass = corsika::particles::GetMass(GetPID());
+      auto P = sqrt((GetEnergy() - mass) * (GetEnergy() + mass));
+      auto p = GetStackData().GetDirection(cs, GetIndex());
+      return p * P;
+    }
+
+    corsika::geometry::Point GetPosition(
+        const corsika::geometry::CoordinateSystem& cs) const {
+      return GetStackData().GetPosition(cs, GetIndex());
+    }
+
+    corsika::geometry::Vector<corsika::units::si::speed_d> GetVelocity(
+        const corsika::geometry::CoordinateSystem& cs) const {
+      return GetStackData().GetVelocity(cs, GetIndex());
+    }
+
+    corsika::units::si::TimeType GetTime() const {
+      return GetStackData().GetTime(GetIndex());
+    }
+
+    corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection(
+        const corsika::geometry::CoordinateSystem& cs) const {
+      return GetStackData().GetDirection(cs, GetIndex());
+    }
+
+    corsika::units::si::TimeType GetTimeInterval() const {
+      return GetStackData().GetTimeInterval();
+    }
+  };
+
+  /**
+   *
+   * Memory implementation of the most simple (stupid) particle stack object.
+   */
+
+  class COASTStackImpl {
+
+    const crs::CParticle* fParticle1 = 0;
+    const crs::CParticle* fParticle2 = 0;
+
+  public:
+    void Init() {}
+    void SetParticle(const crs::CParticle* v1, const crs::CParticle* v2) {
+      fParticle1 = v1;
+      fParticle2 = v2;
+    }
+    void Clear() {}
+
+    // there is one particle only
+    int GetSize() const { return 1; }
+    int GetCapacity() const { return 1; }
+
+    // you cannot modify the particle:
+    // there are no Set... function defined
+
+    // readout particle data, there is just one particle!
+    /*
+      double x;
+      double y;
+      double z;
+      double depth;
+      double time;
+      double energy;
+      double weight;
+      int    particleId;
+      int    hadronicGeneration;
+    */
+    corsika::particles::Code GetPID(const int) const {
+      return ConvertFromCoast(static_cast<CoastCode>(fParticle1->particleId));
+    }
+    corsika::units::si::HEPEnergyType GetEnergy(const int) const {
+      using namespace corsika::units::si;
+      return fParticle1->energy * 1_GeV;
+    }
+    corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection(
+        const corsika::geometry::CoordinateSystem& cs, const int) const {
+      using namespace corsika::units::si;
+      corsika::geometry::Point p1(
+          cs, {fParticle1->x * 1_cm, fParticle1->y * 1_cm, fParticle1->y * 1_cm});
+      corsika::geometry::Point p2(
+          cs, {fParticle2->x * 1_cm, fParticle2->y * 1_cm, fParticle2->y * 1_cm});
+      const corsika::geometry::Vector D = p2 - p1;
+      const auto magD = D.norm();
+      const corsika::geometry::Vector dir = D / magD;
+      return dir;
+    }
+    corsika::geometry::Vector<corsika::units::si::speed_d> GetVelocity(
+        const corsika::geometry::CoordinateSystem& cs, const int) const {
+      using namespace corsika::units::si;
+      corsika::geometry::Vector<corsika::units::si::dimensionless_d> dir =
+          GetDirection(cs, 0);
+      corsika::geometry::Point p1(
+          cs, {fParticle1->x * 1_cm, fParticle1->y * 1_cm, fParticle1->y * 1_cm});
+      corsika::geometry::Point p2(
+          cs, {fParticle2->x * 1_cm, fParticle2->y * 1_cm, fParticle2->y * 1_cm});
+      const corsika::geometry::Vector D = p2 - p1;
+      const LengthType magD = D.norm();
+      const TimeType deltaT = (fParticle2->time - fParticle1->time) * 1_ns;
+      return dir * magD / deltaT;
+    }
+    corsika::geometry::Point GetPosition(const corsika::geometry::CoordinateSystem& cs,
+                                         const int) const {
+      using namespace corsika::units::si;
+      return corsika::geometry::Point(
+          cs, {fParticle1->x * 1_cm, fParticle1->y * 1_cm, fParticle1->z * 1_cm});
+    }
+    corsika::units::si::TimeType GetTime(const int) const {
+      using namespace corsika::units::si;
+      return fParticle1->time * 1_ns;
+    }
+
+    corsika::units::si::TimeType GetTimeInterval() const {
+      using namespace corsika::units::si;
+      return (fParticle2->time - fParticle1->time) * 1_ns;
+    }
+
+    /**
+     *   We do not allow copying!
+     */
+    void Copy(const int, const int) {}
+
+    /**
+     *   We do not allow swapping particles, there is just one...
+     */
+    void Swap(const int, const int) {}
+
+  protected:
+    // size cannot be increased
+    void IncrementSize() {}
+
+    // size cannot be decremented
+    void DecrementSize() {}
+
+  }; // end class COASTStackImpl
+
+  typedef corsika::stack::Stack<COASTStackImpl, ParticleInterface> COASTStack;
+
+} // namespace corsika::coast
+
+#endif

COAST/COASTUserLib.cc

+#include <interface/CorsikaInterface.h>
+
+#include <corsika/coast/COASTProcess.h>
+#include <corsika/coast/COASTStack.h>
+#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <crs/CInteraction.h>
+#include <crs/CParticle.h>
+#include <crs/CorsikaTypes.h>
+#include <crs/TSubBlock.h>
+
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+using namespace corsika;
+using namespace corsika::units::si;
+
+coast::COASTStack gCOASTStack;
+corsika::coast::COASTProcess gCorsikaProcess;
+
+/*
+  Data is one CORSIKA data-block constining of 21 SubBlocks.
+  A SubBlock can be:
+  - thinned mode:     39 (Particles) * 8 (ENTRIES) * 4 (BYTES)
+  - not-thinned mode: 39 (Particles) * 7 (ENTRIES) * 4 (BYTES)
+*/
+extern "C" void wrida_([[maybe_unused]] const CREAL* Data) {
+  // crs::CParticleFortranPtr p;
+  // const bool isF = prminfo_(p);
+}
+
+extern "C" void inida_([[maybe_unused]] const char* filename,
+                       [[maybe_unused]] const int& thinning,
+                       [[maybe_unused]] const int& /*curved*/,
+                       [[maybe_unused]] const int& /*slant*/,
+                       [[maybe_unused]] const int& /*stackinput*/,
+                       [[maybe_unused]] const int& /*preshower*/,
+                       [[maybe_unused]] int str_length) {
+  gCorsikaProcess.Init();
+}
+
+extern "C" void cloda_() {
+  // crs::CParticleFortranPtr pptr;
+  // const bool isF = prminfo_(pptr);
+  // gCorsikaProcess.Close();
+}
+
+void interaction_([[maybe_unused]] const crs::CInteraction& interaction) {
+  /*
+    all interactions in the shower are available in this function !
+    the information availabel in the CInteraction class are:
+    double x;
+    double y;
+    double z;
+    double etot;      // lab energy
+    double sigma;     // cross-section of process
+    double kela;      // elasticity
+    int    projId;    // projectile
+    int    targetId;  // target
+    double time;
+  */
+}
+
+extern "C" void track_([[maybe_unused]] const crs::CParticle& pre,
+                       [[maybe_unused]] const crs::CParticle& post) {
+  /*
+    all particles in the shower are available in this function !
+    The pre and post objecte are the two endpoints for one single track
+    in the shower, where the information available in CParticle is:
+    double x;
+    double y;
+    double z;
+    double depth;
+    double time;
+    double energy;
+    double weight;
+    int    particleId;
+    int    hadronicGeneration;
+  */
+  gCOASTStack.SetParticle(&pre, &post);
+  auto particle = gCOASTStack.GetNextParticle();
+
+  geometry::CoordinateSystem& rootCS =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+  geometry::Line const line(particle.GetPosition(rootCS), particle.GetVelocity(rootCS));
+  const TimeType time = particle.GetTimeInterval();
+  geometry::Trajectory<geometry::Line> track(line, time);
+  gCorsikaProcess.DoContinuous(particle, track, gCOASTStack);
+}
+
+extern "C" void tabularizedatmosphere_([[maybe_unused]] const int& nPoints,
+                                       [[maybe_unused]] const double* height,
+                                       [[maybe_unused]] const double* refractiveIndex) {
+  // for special use only but should be defined because it is delcared in CORSIKA.F
+}

COAST/ParticleConversion.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/coast/ParticleConversion.h>
+#include <corsika/particles/ParticleProperties.h>
+
+#include <exception>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+namespace corsika::coast {
+
+  corsika::particles::Code ConvertFromCoast(CoastCode pCode) {
+    if (coast2corsika.count(pCode) == 0) {
+      ostringstream err;
+      err << "corsika::coast::ConvertFromCoast CoastCode does not exists="
+          << static_cast<CoastCodeIntType>(pCode) << endl;
+      cout << err.str() << endl;
+      throw std::runtime_error(err.str());
+    }
+    return coast2corsika.find(pCode)->second;
+  }
+  //   process::sibyll::Sibyll2Corsika = {
+  //        {PID::E_MINUS, InternalParticleCode::Electron},
+  //};
+} // namespace corsika::coast

COAST/ParticleConversion.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_processes_sibyll_particles_h_
+#define _include_processes_sibyll_particles_h_
+
+#include <corsika/particles/ParticleProperties.h>
+
+#include <bitset2/bitset2.hpp>
+
+#include <map>
+
+namespace corsika::coast {
+
+  enum class CoastCode : int16_t {
+    Gamma = 1,
+    Positron = 2,
+    Electron = 3,
+    MuonBar = 5,
+    Muon = 6,
+    Pi0 = 7,
+    PiP = 8,
+    PiM = 9,
+    Klong = 10,
+    KP = 11,
+    KM = 12,
+    Neutron = 13,
+    Proton = 14,
+    ProtonBar = 15,
+    Kshort = 16,
+    Eta = 17,
+    Lambda = 18,
+    SigmaPlus = 19,
+    Sigma0 = 20,
+    SigmaMinus = 21,
+    Xi0 = 22,
+    XiMinus = 23,
+    OmegaMinus = 24,
+    NeutronBar = 25,
+    LambdaBar = 26,
+    SigmaMinusBar = 27,
+    Sigma0Bar = 28,
+    SigmaPlusBar = 29,
+    Xi0Bar = 30,
+    XiPlusBar = 31,
+    OmegaPlusBar = 32,
+
+    EtaPrime = 48,
+    Phi = 49,
+    omega = 50,
+    Rho0 = 51,
+    RhoPlus = 52,
+    RhoMinus = 53,
+    DeltaPlusPlus = 54,
+    DeltaPlus = 55,
+    Delta0 = 56,
+    DeltaMinus = 57,
+    DeltaMinusMinusBar = 58,
+    DeltaMinusBar = 59,
+    Delta0Bar = 60,
+    DeltaPlusBar = 61,
+    KStar0 = 62,
+    KStarPlus = 63,
+    KStarMinus = 64,
+    KStar0Bar = 65,
+
+    NeutrinoE = 66,
+    NeutrinoEBar = 67,
+    NeutrinoMu = 68,
+    NeutrinoMuBar = 69,
+
+    Code71 = 71,
+    Code72 = 72,
+    Code73 = 73,
+    Code74 = 74,
+    Code75 = 75,
+    Code76 = 76,
+
+    Code85 = 85,
+    Code86 = 86,
+
+    Code95 = 95,
+    Code96 = 96,
+
+    Helium = 402,
+    Carbon = 612,
+    Nitrogen = 714,
+    Oxygen = 816,
+    Iron = 5628,
+  };
+  using CoastCodeIntType = std::underlying_type<CoastCode>::type;
+
+  const std::map<corsika::coast::CoastCode, corsika::particles::Code> coast2corsika = {
+      {CoastCode::Gamma, corsika::particles::Code::Gamma},
+      {CoastCode::Positron, corsika::particles::Code::Positron},
+      {CoastCode::Electron, corsika::particles::Code::Electron},
+      //{CoastCode::    ,corsika::particles::Code::Unknown}, // 4
+      {CoastCode::MuonBar, corsika::particles::Code::MuPlus},
+      {CoastCode::Muon, corsika::particles::Code::MuMinus},
+      {CoastCode::Pi0, corsika::particles::Code::Pi0},
+      {CoastCode::PiP, corsika::particles::Code::PiPlus},
+      {CoastCode::PiM, corsika::particles::Code::PiMinus},
+      {CoastCode::Klong, corsika::particles::Code::K0Long}, // 10
+      {CoastCode::KP, corsika::particles::Code::KPlus},
+      {CoastCode::KM, corsika::particles::Code::KMinus},
+      {CoastCode::Neutron, corsika::particles::Code::Neutron},
+      {CoastCode::Proton, corsika::particles::Code::Proton}, // 14
+      {CoastCode::ProtonBar, corsika::particles::Code::AntiProton},
+      {CoastCode::Kshort, corsika::particles::Code::K0Short},
+      {CoastCode::Eta, corsika::particles::Code::Eta}, // 17
+      {CoastCode::Lambda, corsika::particles::Code::Lambda0},
+      {CoastCode::SigmaPlus, corsika::particles::Code::SigmaPlus},
+      {CoastCode::Sigma0, corsika::particles::Code::Sigma0}, // 20
+      {CoastCode::SigmaMinus, corsika::particles::Code::SigmaMinus},
+      {CoastCode::Xi0, corsika::particles::Code::Xi0},
+      {CoastCode::XiMinus, corsika::particles::Code::XiMinus},
+      {CoastCode::OmegaMinus, corsika::particles::Code::OmegaMinus},
+      {CoastCode::NeutronBar, corsika::particles::Code::AntiNeutron}, // 25
+      {CoastCode::LambdaBar, corsika::particles::Code::Lambda0Bar},
+      {CoastCode::SigmaMinusBar, corsika::particles::Code::SigmaMinusBar},
+      {CoastCode::Sigma0Bar, corsika::particles::Code::Sigma0Bar},
+      {CoastCode::SigmaPlusBar, corsika::particles::Code::SigmaPlusBar},
+      {CoastCode::Xi0Bar, corsika::particles::Code::Xi0Bar},
+      {CoastCode::XiPlusBar, corsika::particles::Code::XiPlusBar},
+      {CoastCode::OmegaPlusBar, corsika::particles::Code::OmegaPlusBar}, // 32
+      //{CoastCode::    ,corsika::particles::Code::Unknown}, // eta-prime
+      //{CoastCode::    ,corsika::particles::Code::Unknown}, // PHI
+      //{CoastCode::    ,corsika::particles::Code::Unknown}, // omega
+      {CoastCode::Rho0, corsika::particles::Code::Rho0}, // 51
+      {CoastCode::RhoPlus, corsika::particles::Code::RhoPlus},
+      {CoastCode::RhoMinus, corsika::particles::Code::RhoMinus},
+      {CoastCode::DeltaPlusPlus, corsika::particles::Code::DeltaPlusPlus},
+      {CoastCode::DeltaPlus, corsika::particles::Code::DeltaPlus},
+      {CoastCode::Delta0, corsika::particles::Code::Delta0}, // 56
+      //{CoastCode::    ,corsika::particles::Code::Unknown}, // DeltaMinus},
+      {CoastCode::DeltaMinusMinusBar, corsika::particles::Code::DeltaMinusMinusBar},
+      {CoastCode::DeltaMinusBar, corsika::particles::Code::DeltaMinusBar},
+      {CoastCode::Delta0Bar, corsika::particles::Code::Delta0Bar},
+      //{CoastCode::    ,corsika::particles::Code::Unknown}, // DeltaPlusBar
+      {CoastCode::KStar0, corsika::particles::Code::KStar0}, // 62
+      {CoastCode::KStarPlus, corsika::particles::Code::KStarPlus},
+      {CoastCode::KStarMinus, corsika::particles::Code::KStarMinus},
+      {CoastCode::KStar0Bar, corsika::particles::Code::KStar0Bar},
+      {CoastCode::NeutrinoE, corsika::particles::Code::NuE},
+      {CoastCode::NeutrinoEBar, corsika::particles::Code::NuEBar},
+      {CoastCode::NeutrinoMu, corsika::particles::Code::NuMu},
+      {CoastCode::NeutrinoMuBar, corsika::particles::Code::NuMuBar}, // 69
+      {CoastCode::Code71, corsika::particles::Code::Unknown},
+      {CoastCode::Code72, corsika::particles::Code::Unknown},
+      {CoastCode::Code73, corsika::particles::Code::Unknown},
+      {CoastCode::Code74, corsika::particles::Code::Unknown},
+      {CoastCode::Code75, corsika::particles::Code::Unknown},
+      {CoastCode::Code76, corsika::particles::Code::Unknown},
+      {CoastCode::Code85, corsika::particles::Code::Unknown},
+      {CoastCode::Code86, corsika::particles::Code::Unknown},
+      {CoastCode::Code95, corsika::particles::Code::Unknown},
+      {CoastCode::Code96, corsika::particles::Code::Unknown},
+
+      {CoastCode::Helium, corsika::particles::Code::Helium}, // 402
+      {CoastCode::Carbon, corsika::particles::Code::Carbon},
+      {CoastCode::Nitrogen, corsika::particles::Code::Nitrogen},
+      {CoastCode::Oxygen, corsika::particles::Code::Oxygen},
+      {CoastCode::Iron, corsika::particles::Code::Iron}, // 5628
+
+      //{CoastCode::, corsika::particles::Code::}
+
+  };
+
+  corsika::particles::Code ConvertFromCoast(CoastCode pCode);
+
+} // namespace corsika::coast
+
+#endif

COAST/README.md

+This is based on corsika7/trunk/coast/CoastOptions/example with an
+additional interface to CORSIKA8. 
+
+This is an example for a "COAST user library" using CORSIKA8
+technology. It explains the steps, how to use the COAST_USER_LIB
+option of CORSIKA together with CORSIKA8/COAST
+
+
+Step 1:
+-------
+
+Configure CORSIKA8 with 'cmake -DWITH_COAST=1' and other options you
+prefer. You have to define COAST_DIR environment variable to the
+location of your existing CORSIKA7 installation. Check that 'ls
+$COAST_DIR/include/interface' contains the file
+'CorsikaInterface.h'. If this is not the case, this is not a valid or
+proper installation of CORSIKA. First compile CORSIKA, e.g. with
+ROOTOUT option, or CoReas to get COAST installed!
+
+
+Step 2:
+-------
+
+Compile CORSIKA8, edit COAST/COASTProcess.cc to modify/add your
+physics module. 
+
+There should be a libCOAST.so as a result! This is what you need.
+
+
+Step 3:
+-------
+
+Create COAST_USER_LIB environment variable to point at your current
+directory, where now the 'libCOAST.so' library is located.
+Add the path in $COAST_USER_LIB to your LD_LIBRARY_PATH environment
+variable.
+
+
+Step 4:
+-------
+
+Go back to your CORSIKA directory and re-start 'coconut'. The option
+COAST_USER_LIB will now be visible. Please select it, copile CORSIKA
+and start the executable. In the generated console output you will the
+statements from your COAST library during Init and Close of the
+simulation.  
+Add any kind of code now working on the CParticle or
+CInteraction class to start using the full power of COAST.
+
+Note: the default COASTProcess just prints out tons of ASCII. This is
+not very useful, don't run a full MC with this...
+
+
+Step 5:
+-------
+
+enjoy...
+

Documentation/Examples/cascade_example.cc

@@ -251,7 +251,8 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-  const HEPEnergyType E0 = 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
+  const HEPEnergyType E0 =
+      100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
   double theta = 0.;
   double phi = 0.;
   {

Framework/Cascade/testCascade.cc

@@ -73,9 +73,10 @@ class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
   int fCount = 0;
   int fCalls = 0;
   HEPEnergyType fEcrit;
-  
+
 public:
-  ProcessSplit(HEPEnergyType e) : fEcrit(e) {}
+  ProcessSplit(HEPEnergyType e)
+      : fEcrit(e) {}
 
   template <typename Particle, typename T>
   LengthType MaxStepLength(Particle&, T&) const {
@@ -83,7 +84,7 @@ public:
   }
 
   template <typename Particle, typename T, typename Stack>
-  EProcessReturn DoContinuous(Particle& p, T&, Stack& s)  {
+  EProcessReturn DoContinuous(Particle& p, T&, Stack& s) {
     fCalls++;
     HEPEnergyType E = p.GetEnergy();
     if (E < fEcrit) {
@@ -102,7 +103,10 @@ public:
     return EProcessReturn::eOk;
   }
 
-  void Init() { fCount = 0; fCalls =0; }
+  void Init() {
+    fCount = 0;
+    fCalls = 0;
+  }
 
   int GetCount() const { return fCount; }
   int GetCalls() const { return fCalls; }
@@ -140,8 +144,8 @@ TEST_CASE("Cascade", "[Cascade]") {
   EAS.Init();
   EAS.Run();
 
-  CHECK( p1.GetCount() == 2048 );
-  CHECK( p1.GetCalls() == 4095 );
+  CHECK(p1.GetCount() == 2048);
+  CHECK(p1.GetCalls() == 4095);
 
   /*
   SECTION("sectionTwo") {

Framework/Geometry/testFourVector.cc

@@ -63,15 +63,15 @@ TEST_CASE("four vectors") {
 
     CHECK(p0.IsSpacelike());
     CHECK(!p0.IsTimelike());
-    //CHECK(!p0.IsPhotonlike());
+    // CHECK(!p0.IsPhotonlike());
 
     CHECK(!p1.IsSpacelike());
     CHECK(p1.IsTimelike());
-    //CHECK(!p1.IsPhotonlike());
+    // CHECK(!p1.IsPhotonlike());
 
     CHECK(!p2.IsSpacelike());
     CHECK(!p2.IsTimelike());
-    //CHECK(p2.IsPhotonlike());
+    // CHECK(p2.IsPhotonlike());
   }
 
   /*

Framework/Particles/NuclearData.xml

@@ -24,4 +24,7 @@
 <particle id="1000180400" name="argon" A="40" Z="18" >
 </particle> 
 
+<particle id="1000280560" name="iron" A="56" Z="28" >
+</particle> 
+
 </chapter>

Framework/Random/testRandom.cc

@@ -12,9 +12,9 @@
                           // cpp file
 #include <catch2/catch.hpp>
 
+#include <corsika/random/ExponentialDistribution.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/random/UniformRealDistribution.h>
-#include <corsika/random/ExponentialDistribution.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <iostream>
 #include <limits>
@@ -93,17 +93,17 @@ TEST_CASE("ExponentialDistribution") {
   auto const beta = 15_m;
 
   corsika::random::ExponentialDistribution dist(beta);
-  
+
   SECTION("mean") {
     std::remove_const<decltype(beta)>::type mean = beta * 0;
-    
+
     int constexpr N = 1'000'000;
-        
+
     for (int i{0}; i < N; ++i) {
       decltype(beta) x = dist(rng);
       mean += x / N;
     }
-    
+
     CHECK(mean / beta == Approx(1).margin(1e-2));
   }
 }

Framework/Utilities/COMBoost.cc

@@ -17,7 +17,8 @@ using namespace corsika::utl;
 using namespace corsika::units::si;
 
 template <typename FourVector>
-COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const HEPMassType massTarget)
+COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile,
+                               const HEPMassType massTarget)
     : fRotation(Eigen::Matrix3d::Identity())
     , fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
   // calculate matrix for rotating pProjectile to z-axis first
@@ -44,8 +45,7 @@ COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const HEPMassType
   }
 
   // calculate boost
-  double const beta =
-    pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
+  double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
 
   /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
   double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
@@ -83,8 +83,8 @@ FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
       (p.GetSpaceLikeComponents().GetComponents().eVector(2) * (1 / 1_GeV).magnitude());
 
   std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << " GeV, "
-            << " pcm=" << p.GetSpaceLikeComponents().GetComponents().squaredNorm() / 1_GeV << " GeV"
-            << std::endl;
+            << " pcm=" << p.GetSpaceLikeComponents().GetComponents().squaredNorm() / 1_GeV
+            << " GeV" << std::endl;
 
   auto const boostedZ = fInverseBoost * com;
   auto const E_lab = boostedZ(0) * 1_GeV;

Framework/Utilities/COMBoost.h

@@ -31,7 +31,8 @@ namespace corsika::utl {
 
   public:
     //! construct a COMBoost given four-vector of prjectile and mass of target
-    COMBoost(const FourVector& Pprojectile, const corsika::units::si::HEPEnergyType massTarget);
+    COMBoost(const FourVector& Pprojectile,
+             const corsika::units::si::HEPEnergyType massTarget);
 
     //! transforms a 4-momentum from lab frame to the center-of-mass frame
     FourVector toCoM(const FourVector& p) const;

Framework/Utilities/CorsikaFenv.h

@@ -27,11 +27,8 @@
  */
 extern "C" {
 
-  int
-  feenableexcept(int excepts);
-  int
-  fedisableexcept(int excepts);
-
+int feenableexcept(int excepts);
+int fedisableexcept(int excepts);
 }
 
-#endif //CORSIKA_CORSIKAFENV_H
+#endif // CORSIKA_CORSIKAFENV_H

Framework/Utilities/CorsikaFenvFallback.cc

@@ -24,14 +24,7 @@
 extern "C" {
 #warning No enabling/disabling of floating point exceptions - platform needs better implementation
 
-  int feenableexcept(int excepts)
-  {
-    return -1;
-  }
-
-  int fedisableexcept(int excepts)
-  {
-    return -1;
-  }
+int feenableexcept(int excepts) { return -1; }
 
+int fedisableexcept(int excepts) { return -1; }
 }

Framework/Utilities/CorsikaFenvOSX.cc

@@ -14,47 +14,42 @@
 #include <cfenv>
 
 // Implementation for OS X on intel X64_86
-// code from https://stackoverflow.com/questions/37819235/how-do-you-enable-floating-point-exceptions-for-clang-in-os-x
-// based on http://www-personal.umich.edu/~williams/archive/computation/fe-handling-example.c
+// code from
+// https://stackoverflow.com/questions/37819235/how-do-you-enable-floating-point-exceptions-for-clang-in-os-x
+// based on
+// http://www-personal.umich.edu/~williams/archive/computation/fe-handling-example.c
 
 extern "C" {
 
-  int
-  feenableexcept(int excepts) {
-    static fenv_t fenv;
-    int new_excepts = excepts & FE_ALL_EXCEPT;
-    // previous masks
-    int old_excepts;
+int feenableexcept(int excepts) {
+  static fenv_t fenv;
+  int new_excepts = excepts & FE_ALL_EXCEPT;
+  // previous masks
+  int old_excepts;
 
-    if (fegetenv(&fenv)) {
-      return -1;
-    }
-    old_excepts = fenv.__control & FE_ALL_EXCEPT;
+  if (fegetenv(&fenv)) { return -1; }
+  old_excepts = fenv.__control & FE_ALL_EXCEPT;
 
-    // unmask
-    fenv.__control &= ~new_excepts;
-    fenv.__mxcsr &= ~(new_excepts << 7);
+  // unmask
+  fenv.__control &= ~new_excepts;
+  fenv.__mxcsr &= ~(new_excepts << 7);
 
-    return fesetenv(&fenv) ? -1 : old_excepts;
-  }
-
-  int
-  fedisableexcept(int excepts) {
-    static fenv_t fenv;
-    int new_excepts = excepts & FE_ALL_EXCEPT;
-    // all previous masks
-    int old_excepts;
+  return fesetenv(&fenv) ? -1 : old_excepts;
+}
 
-    if (fegetenv(&fenv)) {
-      return -1;
-    }
-    old_excepts = fenv.__control & FE_ALL_EXCEPT;
+int fedisableexcept(int excepts) {
+  static fenv_t fenv;
+  int new_excepts = excepts & FE_ALL_EXCEPT;
+  // all previous masks
+  int old_excepts;
 
-    // mask
-    fenv.__control |= new_excepts;
-    fenv.__mxcsr |= new_excepts << 7;
+  if (fegetenv(&fenv)) { return -1; }
+  old_excepts = fenv.__control & FE_ALL_EXCEPT;
 
-    return fesetenv(&fenv) ? -1 : old_excepts;
-  }
+  // mask
+  fenv.__control |= new_excepts;
+  fenv.__mxcsr |= new_excepts << 7;
 
+  return fesetenv(&fenv) ? -1 : old_excepts;
+}
 }

Framework/Utilities/testCOMBoost.cc

@@ -36,7 +36,7 @@ TEST_CASE("boosts") {
   auto energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
     return sqrt(m * m + p.squaredNorm());
   };
-  
+
   // helper function for mandelstam-s
   auto s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
     return E * E - p.squaredNorm();
@@ -50,7 +50,7 @@ TEST_CASE("boosts") {
   /*
     General tests check the interface and basic operation
    */
-  
+
   SECTION("General tests") {
 
     // define projectile kinematics in lab frame
@@ -124,7 +124,7 @@ TEST_CASE("boosts") {
   /*
     special case: projectile with arbitrary direction
    */
-  
+
   SECTION("Test boost along tilted axis") {
 
     const HEPMomentumType P0 = 1_PeV;
@@ -161,7 +161,7 @@ TEST_CASE("boosts") {
   /*
     test the ultra-high energy behaviour: E=ZeV
    */
-  
+
   SECTION("High energy") {
     // define projectile kinematics in lab frame
     HEPMassType const projectileMass = 1_GeV;