Framework/Particles/ParticleProperties.h

-/**
-   @file Particles.h
-
-   Interface to particle properties
- */
-
-#ifndef _include_Particle_h_
-#define _include_Particle_h_
-
-#include <array>
-#include <cstdint>
-#include <iostream>
-
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/particles/GeneratedParticleProperties.inc>
-
-/**
- * @namespace particle
- *
- * The properties of all elementary particles is stored here. The data
- * is taken from the Pythia ParticleData.xml file.
- *
- */
-
-namespace corsika::particles {
-
-  /**
-   * @function GetMass
-   *
-   * return mass of particle
-   */
-  auto constexpr GetMass(Code const p) { return masses[static_cast<uint8_t const>(p)]; }
-
-  auto constexpr GetPDG(Code const p) { return pdg_codes[static_cast<uint8_t const>(p)]; }
-
-  auto constexpr GetElectricChargeNumber(Code const p) {
-    return electric_charge[static_cast<uint8_t const>(p)] / 3;
-  }
-
-  auto constexpr GetElectricCharge(Code const p) {
-    return GetElectricChargeNumber(p) * (corsika::units::constants::e);
-  }
-
-  auto const GetName(Code const p) { return names[static_cast<uint8_t const>(p)]; }
-
-  namespace io {
-
-    std::ostream& operator<<(std::ostream& stream, Code const p) {
-      stream << GetName(p);
-      return stream;
-    }
-
-  } // namespace io
-
-} // namespace corsika::particles
-
-// to inject the operator<< into the root namespace
-using namespace corsika::particles::io;
-
-#endif

Framework/Particles/pdxml_reader.py

-#!/usr/bin/env python3
-
-import sys, math, itertools, re, csv, pprint
-import xml.etree.ElementTree as ET
-from collections import OrderedDict
-
-
-
-##############################################################
-# 
-# reading xml input data, return line by line particle data
-# 
-
-def parse(filename):
-    tree = ET.parse(filename)
-    root = tree.getroot()
-        
-    for particle in root.iter("particle"):
-        name = particle.attrib["name"]
-        antiName = "Unknown"
-        # print (str(particle.attrib))
-        if ("antiName" in particle.attrib):
-            antiName = particle.attrib["antiName"]            
-#            print ("found anti: " + name + " " + antiName + "\n" )
-        pdg_id = int(particle.attrib["id"])
-        mass = float(particle.attrib["m0"]) # GeV
-        electric_charge = int(particle.attrib["chargeType"]) # in units of e/3
-        
-        decay_width = float(particle.attrib.get("mWidth", 0)) # GeV
-        lifetime = float(particle.attrib.get("tau0", math.inf)) # mm / c
-        
-        yield (pdg_id, name, mass, electric_charge, antiName)
-                
-        # TODO: read decay channels from child elements
-        
-        if "antiName" in particle.attrib:
-            yield (-pdg_id, antiName, mass, -electric_charge, name)
-
-
-            
-##############################################################
-# 
-# returns dict with particle codes and class names
-# 
-
-def class_names(filename):
-    tree = ET.parse(filename)
-    root = tree.getroot()
-
-    map = {}
-    
-    for particle in root.iter("particle"):
-        name = particle.attrib["classname"]
-        pdg_id = int(particle.attrib["id"])
-        map[pdg_id] = name
-        
-    return map 
-
-
-
-##############################################################
-# 
-# Automatically produce a string qualifying as C++ class name
-# 
-# This function produces names of type "DELTA_PLUS_PLUS"
-# 
-def c_identifier(name):
-    orig = name
-    name = name.upper()
-    for c in "() /":
-        name = name.replace(c, "_")
-    
-    name = name.replace("BAR", "_BAR")
-    name = name.replace("0", "_0")
-    name = name.replace("*", "_STAR")
-    name = name.replace("'", "_PRIME")
-    name = name.replace("+", "_PLUS")
-    name = name.replace("-", "_MINUS")
-    
-    while True:
-        tmp = name.replace("__", "_")
-        if tmp == name:
-            break
-        else:
-            name = tmp    
-
-    pattern = re.compile(r'^[A-Z_][A-Z_0-9]*$')
-    if pattern.match(name):
-        return name.strip("_")
-    else:
-        raise Exception("could not generate C identifier for '{:s}'".format(orig))
-
-
-##############################################################
-# 
-# Automatically produce a string qualifying as C++ class name
-# 
-# This function produces names of type "DeltaPlusPlus"
-# 
-def c_identifier_cammel(name):
-    orig = name
-    name = name[0].upper() + name[1:].lower() # all lower case
-    for c in "() /": # replace funny characters
-        name = name.replace(c, "_")
-    
-    name = name.replace("bar", "Bar")
-    name = name.replace("*", "Star")
-    name = name.replace("'", "Prime")
-    name = name.replace("+", "Plus")
-    name = name.replace("-", "Minus")
-
-    # move "Bar" to end of name
-    ibar = name.find('Bar')
-    if (ibar>0 and ibar<len(name)-3) :
-        name = name[:ibar] + name[ibar+3:] + str('Bar')
-            
-    # cleanup "_"s
-    while True:
-        tmp = name.replace("__", "_")
-        if tmp == name:
-            break
-        else:
-            name = tmp
-    name.strip("_")
-
-    # remove all "_", if this does not by accident concatenate two number
-    istart = 0
-    while True:
-        i = name.find('_', istart)
-        if (i<1 or i>len(name)-1):
-            break
-        istart = i
-        if (name[i-1].isdigit() and name[i+1].isdigit()):
-            # there is a number on both sides
-            break
-        name = name[:i] + name[i+1:]
-        # and last, for example: make NuE out of Nue
-        if (name[i-1].islower() and name[i].islower()) :
-            if (i<len(name)-1) :
-                name = name[:i] + name[i].upper() + name[i+1:]
-            else :
-                name = name[:i] + name[i].upper()
-
-    # check if name is valid C++ identifier
-    pattern = re.compile(r'^[a-zA-Z_][a-zA-Z_0-9]*$')
-    if pattern.match(name):
-        return name
-    else:
-        raise Exception("could not generate C identifier for '{:s}' name='{:s}'".format(orig, name))
-
-
-    
-# ########################################################
-# 
-# returns dict containing all data from pythia-xml input
-# 
-    
-def build_pythia_db(filename, classnames):    
-    particle_db = OrderedDict()    
-    for (pdg, name, mass, electric_charge, antiName) in parse(filename):
-
-        c_id = "unknown"
-        if (pdg in classnames):
-            c_id = classnames[pdg]
-        else:
-            c_id = c_identifier_cammel(name) # the cammel case names
-        
-        #~ print(name, c_id, sep='\t', file=sys.stderr)
-        #~ enums += "{:s} = {:d}, ".format(c_id, corsika_id)
-        particle_db[c_id] = {
-            "name" : name,
-            "antiName" : antiName,
-            "pdg" : pdg,
-            "mass" : mass, # in GeV
-            "electric_charge" : electric_charge # in e/3
-        }
-    
-    return particle_db
-    
-
-
-
-###############################################################
-# 
-# return string with enum of all internal particle codes
-# 
-
-def gen_internal_enum(pythia_db):
-    string = "enum class Code : uint8_t {\n"
-
-    string += "  FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\"  \n" # identifier for eventual loops...
-    last_ngc_id = 0 
-    for k in filter(lambda k: "ngc_code" in pythia_db[k], pythia_db):
-        last_ngc_id = pythia_db[k]['ngc_code']
-        string += "  {key:s} = {code:d},\n".format(key = k, code = last_ngc_id)
-
-    string += "  LastParticle = " + str(last_ngc_id+1) + ",\n" # identifier for eventual loops...
-    
-    string += "};"
-    return string
-
-
-
-
-###############################################################
-# 
-# return string with all data arrays 
-# 
-
-def gen_properties(pythia_db):
-
-    # number of particles, size of tables
-    string = "static constexpr std::size_t size = {size:d};\n".format(size = len(pythia_db))
-    string += "\n"
-    
-    # particle masses table
-    string += "static constexpr std::array<quantity<energy_d> const, size> masses = {\n"    
-    for p in pythia_db.values():
-        string += "  {mass:f}_GeV, // {name:s}\n".format(mass = p['mass'], name = p['name'])              
-    string += "};\n\n"
-                   
-    # PDG code table
-    string += "static constexpr std::array<PDGCode const, size> pdg_codes = {\n"
-    for p in pythia_db.values():
-        string += "  {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name'])    
-    string += "};\n"
-    
-    # name string table
-    string += "static const std::array<std::string const, size> names = {\n"
-    for p in pythia_db.values():
-        string += "  \"{name:s}\",\n".format(name = p['name'])            
-    string += "};\n"
-    
-    # electric charges table
-    string += "static constexpr std::array<int16_t, size> electric_charge = {\n"
-    for p in pythia_db.values():
-        string += "  {charge:d},\n".format(charge = p['electric_charge'])
-    string += "};\n"
-
-    # anti-particle table
-#    string += "static constexpr std::array<size, size> anti_particle = {\n"
-#    for p in pythia_db.values():
-#        string += "  {anti:d},\n".format(charge = p['anti_particle'])
-#    string += "};\n"
-
-    return string
-
-
-
-
-###############################################################
-# 
-# return string with a list of classes for all particles
-# 
-
-def gen_classes(pythia_db):
-
-    string = "// list of C++ classes to access particle properties"
-    
-    for cname in pythia_db:
-
-        antiP = 'unknown'
-        for cname_anti in pythia_db:
-            if (pythia_db[cname_anti]['name'] == pythia_db[cname]['antiName']):
-                antiP = cname_anti
-                break
-        
-        string += "\n";
-        string += "/** @class " + cname + "\n\n"
-        string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\n"
-        string += " *  - pdg=" + str(pythia_db[cname]['pdg']) +"\n"
-        string += " *  - mass=" + str(pythia_db[cname]['mass']) + " GeV \n"
-        string += " *  - charge= " + str(pythia_db[cname]['electric_charge']/3) + " \n"
-        string += " *  - name=" + str(cname) + "\n"
-        string += " *  - anti=" + str(antiP) + "\n"
-        string += "*/\n\n"
-        string += "class " + cname + "{\n"
-        string += "  public:\n"
-        string += "   static Code GetCode() { return Type; }\n"
-        string += "   static quantity<energy_d> GetMass() { return masses[TypeIndex]; }\n"
-        string += "   static quantity<electric_charge_d> GetCharge() { return corsika::units::constants::e*electric_charge[TypeIndex]/3; }\n"
-        string += "   static int GetChargeNumber() { return electric_charge[TypeIndex]/3; }\n"
-        string += "   static std::string GetName() { return names[TypeIndex]; }\n"
-        string += "   static Code GetAntiParticle() { return AntiType; }\n"
-        string += "   static const Code Type = Code::" + cname + ";\n"
-        string += "   static const Code AntiType = Code::" + antiP + ";\n"
-        string += " private:\n"
-        string += "   static const uint8_t TypeIndex = static_cast<uint8_t const>(Type);\n"
-        string += "};\n"
-
-    return string
-
-
-###############################################################
-# 
-# 
-
-def inc_start():
-    string = ""
-    string += "#ifndef _include_GeneratedParticleDataTable_h_\n"
-    string += "#define _include_GeneratedParticleDataTable_h_\n\n"
-    string += "#include <corsika/units/PhysicalUnits.h>\n"
-    string += "#include <corsika/units/PhysicalConstants.h>\n"
-    string += "#include <array>\n"
-    string += "#include <cstdint>\n"
-#    string += "#include <iostream>\n\n"
-    string += "namespace corsika { \n\n"
-#    string += "using namespace literals; \n"
-    string += "namespace particles { \n\n"
-    string += "using corsika::units::energy_d;\n"
-    string += "using corsika::units::electric_charge_d;\n"
-    string += "using corsika::units::quantity;\n"
-    string += "using corsika::units::operator\"\"_GeV;\n"
-    string += "typedef int16_t PDGCode;\n\n"
-    return string
-
-
-###############################################################
-# 
-# 
-
-def inc_end():
-    string = ""
-    string += "\n}\n\n"
-    string += "\n}\n\n"
-    string += "#endif\n"
-    return string
-
-
-###################################################################
-# 
-# Main function
-# 
-
-if __name__ == "__main__":
-
-    if (len(sys.argv)!=3) :
-        print ("pdxml_reader.py Pythia8.xml ClassNames.xml")
-        sys.exit(0)
-        
-    names = class_names(sys.argv[2])
-    pythia_db = build_pythia_db(sys.argv[1], names)
-
-    print ("\n       pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n")
-    
-    counter = itertools.count(0)
-    
-    not_modeled = []
-    for p in pythia_db:
-        pythia_db[p]['ngc_code'] = next(counter)               
-    
-    with open("GeneratedParticleProperties.inc", "w") as f:
-        print(inc_start(), file=f) 
-        print(gen_internal_enum(pythia_db), file=f)
-        print(gen_properties(pythia_db), file=f)
-        print(gen_classes(pythia_db), file=f)
-        print(inc_end(), file=f) 
-    
-    #~ print(pdg_id_table, mass_table, name_table, enums, sep='\n\n')
-
-
-    

Framework/Particles/testParticles.cc

-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-using namespace corsika::units;
-using namespace corsika::particles;
-
-TEST_CASE("Particles", "[Particles]") {
-
-  SECTION("Types") { REQUIRE(Electron::GetCode() == Code::Electron); }
-
-  SECTION("Data") {
-    REQUIRE(Electron::GetMass() / 0.511_MeV == Approx(1));
-    REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
-    REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
-    REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
-    REQUIRE(GetElectricCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
-    REQUIRE(Electron::GetName() == "e-");
-  }
-}

Framework/ProcessSequence/CMakeLists.txt

-
-add_library (CORSIKAprocesssequence INTERFACE)
-
-# namespace of library -> location of header files
-set (
-  CORSIKAprocesssequence_NAMESPACE
-  corsika/process
-  )
-
-# header files of this library
-set (
-  CORSIKAprocesssequence_HEADERS
-  ProcessSequence.h
-  )
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAprocesssequence ${CORSIKAprocesssequence_NAMESPACE} ${CORSIKAprocesssequence_HEADERS})
-
-# include directive for upstream code
-target_include_directories (
-  CORSIKAprocesssequence
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-# install library
-install (
-  FILES ${CORSIKAprocesssequence_HEADERS}
-  DESTINATION include/${CORSIKAprocesssequence_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-#add_executable (
-#  testLogging
-#  testLogging.cc
-#  )
-
-#target_link_libraries (
-#  testLogging
-#  CORSIKAprocesssequence
-#  CORSIKAthirdparty # for catch2
-#  )
-
-#add_test (
-#  NAME testLogging
-#  COMMAND testLogging
-#  )
-
-
-  

Framework/ProcessSequence/ProcessSequence.h

-#ifndef _include_ProcessSequence_h_
-#define _include_ProcessSequence_h_
-
-#include <iostream>
-#include <typeinfo>
-
-namespace corsika::process {
-
-  /**
-     /class BaseProcess
-
-     The structural base type of a process object in a
-     ProcessSequence. Both, the ProcessSequence and all its elements
-     are of type BaseProcess<T>
-
-   */
-
-  template <typename derived>
-  struct BaseProcess {
-    const derived& GetRef() const { return static_cast<const derived&>(*this); }
-  };
-
-  /**
-     \class ProcessSequence
-
-     A compile time static list of processes. The compiler will
-     generate a new type based on template logic containing all the
-     elements.
-
-     \comment Using CRTP pattern,
-     https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
-   */
-
-  template <typename T1, typename T2>
-  class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
-  public:
-    const T1& A;
-    const T2& B;
-
-    ProcessSequence(const T1& in_A, const T2& in_B)
-        : A(in_A)
-        , B(in_B) {}
-
-    template <typename D>
-    inline void DoContinuous(D& d) const {
-      A.DoContinuous(d);
-      B.DoContinuous(d);
-    } // add trajectory
-
-    template <typename D>
-    inline double MinStepLength(D& d) const {
-      return min(A.MinStepLength(d), B.MinStepLength(d));
-    }
-
-    // template<typename D>
-    // inline Trajectory Transport(D& d, double& length) const { A.Transport(d, length);
-    // B.Transport(d, length); }
-
-    template <typename D>
-    inline void DoDiscrete(D& d) const {
-      A.DoDiscrete(d);
-      B.DoDiscrete(d);
-    }
-  };
-
-  template <typename T1, typename T2>
-  inline const ProcessSequence<T1, T2> operator+(const BaseProcess<T1>& A,
-                                                 const BaseProcess<T2>& B) {
-    return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef());
-  }
-
-  /*
-    template <typename T1>
-    struct depth_lhs
-    {
-    static const int num = 0;
-    };
-
-
-
-    // terminating condition
-    template <typename T1, typename T2>
-    struct depth_lhs< Sequence<T1,T2> >
-    {
-    // try to expand the left node (T1) which might be a Sequence type
-    static const int num = 1 + depth_lhs<T1>::num;
-    };
-  */
-
-  /*
-    template <typename T1>
-    struct mat_ptrs
-    {
-    static const int num = 0;
-
-    inline static void
-    get_ptrs(const Process** ptrs, const T1& X)
-    {
-    ptrs[0] = reinterpret_cast<const Process*>(&X);
-    }
-    };
-
-
-    template <typename T1, typename T2>
-    struct mat_ptrs< Sequence<T1,T2> >
-    {
-    static const int num = 1 + mat_ptrs<T1>::num;
-
-    inline static void
-    get_ptrs(const Process** in_ptrs, const Sequence<T1,T2>& X)
-    {
-    // traverse the left node
-    mat_ptrs<T1>::get_ptrs(in_ptrs, X.A);
-    // get address of the matrix on the right node
-    in_ptrs[num] = reinterpret_cast<const Process*>(&X.B);
-    }
-    };
-  */
-
-  /*
-    template<typename T1, typename T2>
-    const Process&
-    Process::operator=(const Sequence<T1,T2>& X)
-    {
-    int N = 1 + depth_lhs< Sequence<T1,T2> >::num;
-    const Process* ptrs[N];
-    mat_ptrs< Sequence<T1,T2> >::get_ptrs(ptrs, X);
-    int r = ptrs[0]->rows;
-    int c = ptrs[0]->cols;
-    // ... check that all matrices have the same size ...
-    set_size(r, c);
-    for(int j=0; j<r*c; ++j)
-    {
-    double sum = ptrs[0]->data[j];
-    for(int i=1; i<N; ++i)
-    {
-    sum += ptrs[i]->data[j];
-    }
-    data[j] = sum;
-    }
-    return *this;
-    }
-  */
-
-} // namespace corsika::process
-
-#endif

Framework/StackInterface/CMakeList.txt

-
-add_library (CORSIKAstack INTERFACE)
-
-target_include_directories (CORSIKAstack INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
-                                                   $<INSTALL_INTERFACE:include/Framework>
-                                                   )
-
-install (FILES Stack.h StackIterator.h
-         DESTINATION include/Stack)
-

Framework/StackInterface/CMakeLists.txt

-set (
-  CORSIKAstackinterface_HEADERS
-  Stack.h
-  StackIterator.h
-  )
-
-set (
-  CORSIKAstackinterface_NAMESPACE
-  corsika/stack
-  )
-
-add_library (
-  CORSIKAstackinterface
-  INTERFACE
-  )
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAstackinterface ${CORSIKAstackinterface_NAMESPACE} ${CORSIKAstackinterface_HEADERS})
-
-target_include_directories (
-  CORSIKAstackinterface
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (
-  FILES
-  ${CORSIKAstackinterface_HEADERS}
-  DESTINATION
-  include/${CORSIKAstackinterface_NAMESPACE}
-  )

Framework/StackInterface/Stack.h

-#ifndef _include_Stack_h__
-#define _include_Stack_h__
-
-#include <corsika/stack/StackIterator.h> // include here, to help application programmres
-
-/**
-   All classes around management of particles on a stack.
- */
-
-namespace corsika::stack {
-
-  /**
-     Interface definition of a Stack object. The Stack implements the
-     std-type begin/end function to allow integration in normal for
-     loops etc.
-   */
-
-  template <typename DataImpl, typename Particle>
-  class Stack : public DataImpl {
-
-  public:
-    using DataImpl::GetCapacity;
-    using DataImpl::GetSize;
-
-    using DataImpl::Clear;
-    using DataImpl::Copy;
-
-    using DataImpl::DecrementSize;
-    using DataImpl::IncrementSize;
-
-  public:
-    typedef Particle iterator;
-    typedef const Particle const_iterator;
-
-    /// these are functions required by std containers and std loops
-    iterator begin() { return iterator(*this, 0); }
-    iterator end() { return iterator(*this, GetSize()); }
-    iterator last() { return iterator(*this, GetSize() - 1); }
-
-    /// these are functions required by std containers and std loops
-    const_iterator cbegin() const { return const_iterator(*this, 0); }
-    const_iterator cend() const { return const_iterator(*this, GetSize()); }
-    const_iterator clast() const { return const_iterator(*this, GetSize() - 1); }
-
-    /// increase stack size, create new particle at end of stack
-    iterator NewParticle() {
-      IncrementSize();
-      return iterator(*this, GetSize() - 1);
-    }
-    /// delete last particle on stack by decrementing stack size
-    void DeleteLast() { DecrementSize(); }
-  };
-
-} // namespace corsika::stack
-
-#endif

Framework/StackInterface/StackIterator.h

-#ifndef _include_StackIterator_h__
-#define _include_StackIterator_h__
-
-#include <iomanip>
-#include <iostream>
-
-namespace corsika::stack {
-
-  // forward decl.
-  template <class Stack, class Particle>
-  class StackIteratorInfo;
-
-  /**
-     @class StackIterator
-
-     The StackIterator is the main interface to iterator over
-     particles on a stack. At the same time StackIterator is a
-     Particle object by itself, thus there is no difference between
-     type and ref_type for convenience of the physicist.
-
-     This allows to write code like
-     \verbatim
-     for (auto& p : theStack) { p.SetEnergy(newEnergy); }
-     \endverbatim
-
-     The template argument Stack determines the type of Stack object
-     the data is stored in. A pointer to the Stack object is part of
-     the StackIterator. In addition to Stack the iterator only knows
-     the index fIndex in the Stack data.
-
-     The template argument Particles acts as a policy to provide
-     readout function of Particle data from the stack. The Particle
-     class must know how to retrieve information from the Stack data
-     for a particle entry at any index fIndex.
-  */
-
-  template <typename Stack, typename Particle>
-  class StackIterator : public Particle {
-    friend Stack;
-    friend Particle;
-    friend StackIteratorInfo<Stack, Particle>;
-
-  private:
-    int fIndex;
-
-    //#warning stacks should not be copied because of this:
-    Stack* fData;
-
-  public:
-    StackIterator()
-        : fData(0)
-        , fIndex(0) {}
-    StackIterator(Stack& data, const int index)
-        : fData(&data)
-        , fIndex(index) {}
-    StackIterator(const StackIterator& mit)
-        : fData(mit.fData)
-        , fIndex(mit.fIndex) {}
-
-    StackIterator& operator++() {
-      ++fIndex;
-      return *this;
-    }
-    StackIterator operator++(int) {
-      StackIterator tmp(*this);
-      ++fIndex;
-      return tmp;
-    }
-    bool operator==(const StackIterator& rhs) { return fIndex == rhs.fIndex; }
-    bool operator!=(const StackIterator& rhs) { return fIndex != rhs.fIndex; }
-
-    StackIterator& operator*() { return *this; }
-    const StackIterator& operator*() const { return *this; }
-
-  protected:
-    int GetIndex() const { return fIndex; }
-    Stack& GetStack() { return *fData; }
-    const Stack& GetStack() const { return *fData; }
-
-    // this is probably not needed rigth now:
-    // inline StackIterator<Stack,Particle>& BaseRef() { return
-    // static_cast<StackIterator<Stack, Particle>&>(*this); } inline const
-    // StackIterator<Stack,Particle>& BaseRef() const { return static_cast<const
-    // StackIterator<Stack, Particle>&>(*this); }
-  };
-
-  /**
-     @class StackIteratorInfo
-
-     This is the class where custom ...
-     Internal helper class for StackIterator. Document better...
-   */
-
-  template <typename _Stack, typename Particle>
-  class StackIteratorInfo {
-
-    friend Particle;
-
-  private:
-    StackIteratorInfo() {}
-
-  protected:
-    inline _Stack& GetStack() {
-      return static_cast<StackIterator<_Stack, Particle>*>(this)->GetStack();
-    }
-    inline int GetIndex() const {
-      return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetIndex();
-    }
-    inline const _Stack& GetStack() const {
-      return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetStack();
-    }
-  };
-
-} // namespace corsika::stack
-
-#endif

Framework/Units/CMakeLists.txt

-
-add_library (CORSIKAunits INTERFACE)
-
-set (CORSIKAunits_NAMESPACE corsika/units)
-set (
-  CORSIKAunits_HEADERS
-  PhysicalUnits.h
-  PhysicalConstants.h
-  )
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAunits ${CORSIKAunits_NAMESPACE} ${CORSIKAunits_HEADERS})
-
-target_include_directories (CORSIKAunits
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/ThirdParty>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (FILES ${CORSIKAunits_HEADERS} DESTINATION include/${CORSIKAunits_NAMESPACE})
-
-# code testing
-add_executable (testUnits testUnits.cc)
-target_link_libraries (testUnits CORSIKAunits CORSIKAthirdparty) # for catch2
-add_test(NAME testUnits COMMAND testUnits)
-

Framework/Units/PhysicalConstants.h

-/**
- * \file PhysicalConstants
- *
- * \brief   Several physical constants.
- * \author  Michael S. Kenniston, Martin Moene
- * \date    7 September 2013
- * \since   0.4
- *
- * Copyright 2013 Universiteit Leiden. All rights reserved.
- *
- * Copyright (c) 2001 by Michael S. Kenniston.  For the most
- * recent version check www.xnet.com/~msk/quantity.  Permission is granted
- * to use this code without restriction so long as this copyright
- * notice appears in all source files.
- *
- * This code is provided as-is, with no warrantee of correctness.
- *
- * Distributed under the Boost Software License, Version 1.0. (See accompanying
- * file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- *
- *
- *
- */
-
-#ifndef INCLUDE_PHYSICAL_CONSTANTS_H
-#define INCLUDE_PHYSICAL_CONSTANTS_H
-
-#include <phys/units/quantity.hpp>
-
-namespace corsika::units::constants {
-
-  using namespace phys::units;
-
-  // acceleration of free-fall, standard
-  constexpr phys::units::quantity<phys::units::acceleration_d> g_sub_n{
-      phys::units::Rep(9.80665L) * phys::units::meter /
-      phys::units::square(phys::units::second)};
-
-  // Avogadro constant
-  constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) / mole};
-  // electronvolt
-  constexpr quantity<energy_d> eV{Rep(1.60217733e-19L) * joule};
-
-  // elementary charge
-  constexpr quantity<electric_charge_d> e{Rep(1.602176462e-19L) * coulomb};
-
-  // Planck constant
-  constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second};
-
-  // speed of light in a vacuum
-  constexpr quantity<speed_d> c{Rep(299792458L) * meter / second};
-
-  // unified atomic mass unit
-  constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
-
-  // etc.
-
-} // namespace corsika::units::constants
-
-#endif // PHYS_UNITS_PHYSICAL_CONSTANTS_HPP_INCLUDED

Framework/Units/PhysicalUnits.h

-#ifndef _include_PhysicalUnits_h_
-#define _include_PhysicalUnits_h_
-
-#include <corsika/units/PhysicalConstants.h>
-
-#include <phys/units/io.hpp>
-#include <phys/units/quantity.hpp>
-
-/**
- * @file PhysicalUnits
- *
- * Define _XeV literals, alowing 10_GeV in the code.
- */
-
-namespace phys {
-  namespace units {
-    namespace literals {
-      QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d,
-                                       magnitude(corsika::units::constants::eV))
-    }
-  } // namespace units
-} // namespace phys
-
-namespace corsika::units {
-
-  using namespace phys::units;
-  using namespace phys::units::literals;
-  // namespace literals = phys::units::literals;
-
-  using LengthType = phys::units::quantity<phys::units::length_d, double>;
-  using TimeType = phys::units::quantity<phys::units::time_interval_d, double>;
-  using SpeedType = phys::units::quantity<phys::units::speed_d, double>;
-  using FrequencyType = phys::units::quantity<phys::units::frequency_d, double>;
-  using ElectricChargeType =
-      phys::units::quantity<phys::units::electric_charge_d, double>;
-  using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
-
-} // end namespace corsika::units
-
-// we want to call the operator<< without namespace... I think
-using namespace phys::units::io;
-
-#endif

Framework/Units/testUnits.cc

-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <array>
-
-using namespace corsika::units;
-
-TEST_CASE("PhysicalUnits", "[Units]") {
-
-  SECTION("Consistency") {
-    REQUIRE(1_m / 1_m == Approx(1));
-    // REQUIRE_FALSE( 1_m/1_s == 1 ); // static assert
-  }
-
-  SECTION("Constructors") {
-    auto E1 = 10_GeV;
-    REQUIRE(E1 == 10_GeV);
-
-    LengthType l1 = 10_nm;
-
-    LengthType arr0[5];
-    arr0[0] = 5_m;
-
-    LengthType arr1[2] = {{1_mm}, {2_cm}};
-
-    std::array<EnergyType, 4> arr2; // empty array
-
-    std::array<EnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
-  }
-
-  SECTION("Powers in literal units") {
-    REQUIRE(1_s / 1_ds == Approx(1e1));
-    REQUIRE(1_m / 1_cm == Approx(1e2));
-    REQUIRE(1_m / 1_mm == Approx(1e3));
-    REQUIRE(1_V / 1_uV == Approx(1e6));
-    REQUIRE(1_s / 1_ns == Approx(1e9));
-    REQUIRE(1_eV / 1_peV == Approx(1e12));
-    REQUIRE(1_A / 1_fA == Approx(1e15));
-    REQUIRE(1_mol / 1_amol == Approx(1e18));
-    REQUIRE(1_K / 1_zK == Approx(1e21));
-    REQUIRE(1_K / 1_yK == Approx(1e24));
-
-    REQUIRE(1_A / 1_hA == Approx(1e-2));
-    REQUIRE(1_m / 1_km == Approx(1e-3));
-    REQUIRE(1_m / 1_Mm == Approx(1e-6));
-    REQUIRE(1_V / 1_GV == Approx(1e-9));
-    REQUIRE(1_s / 1_Ts == Approx(1e-12));
-    REQUIRE(1_eV / 1_PeV == Approx(1e-15));
-    REQUIRE(1_A / 1_EA == Approx(1e-18));
-    REQUIRE(1_K / 1_ZK == Approx(1e-21));
-    REQUIRE(1_mol / 1_Ymol == Approx(1e-24));
-  }
-
-  SECTION("Powers and units") {
-    REQUIRE(1 * ampere / 1_A == Approx(1e0));
-    REQUIRE(mega * bar / bar == Approx(1e6));
-  }
-
-  SECTION("Formulas") {
-    const EnergyType E2 = 20_GeV * 2;
-    REQUIRE(E2 == 40_GeV);
-
-    const double lgE = log10(E2 / 1_GeV);
-    REQUIRE(lgE == Approx(log10(40.)));
-
-    const auto E3 = E2 + 100_GeV + pow(10, lgE) * 1_GeV;
-    REQUIRE(E3 == 180_GeV);
-  }
-}

Main/CMakeLists.txt

-
-

Main/shower.cc

-#include <utl/Stack.h>
-
-int main(int argc, char** argv) { return 0; }

Processes/CMakeLists.txt

-add_subdirectory (NullModel)
-add_subdirectory (Sibyll)

Processes/NullModel/CMakeLists.txt

-
-set (
-  MODEL_SOURCES
-  NullModel.cc
-  )
-
-set (
-  MODEL_HEADERS
-  NullModel.h
-  )
-
-set (
-  MODEL_NAMESPACE
-  corsika/process/null_model
-  )
-
-add_library (ProcessNullModel STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessNullModel ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-set_target_properties (
-  ProcessNullModel
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-#  PUBLIC_HEADER "${MODEL_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  ProcessNullModel
-  CORSIKAunits
-  )
-
-target_include_directories (
-  ProcessNullModel 
-  INTERFACE 
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS ProcessNullModel
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-#  PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-add_executable (testNullModel testNullModel.cc)
-
-target_link_libraries (
-  testNullModel
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (NAME testNullModel COMMAND testNullModel)
-