Framework/Logging/NoSink.h

-#ifndef _include_NoSink_h_
-#define _include_NoSink_h_
-
-namespace corsika::logging {
-
-  namespace sink {
-
-    struct NoSink {
-      inline void operator<<(const std::string&) {}
-      inline void Close() {}
-    };
-
-  } // namespace sink
-} // namespace corsika::logging
-
-#endif

Framework/Logging/Sink.h

-#ifndef _include_Sink_h_
-#define _include_Sink_h_
-
-namespace corsika::logging {
-
-  /**
-     a sink for the logger must implement the two functions
-     operator<<(const std::string&)
-     and
-     Close()
-
-     See example: NoSink
-   */
-
-  namespace sink {
-
-    /**
-       Definition of Sink for log output.
-    */
-    template <typename TStream>
-    class Sink {
-    public:
-      Sink(TStream& out)
-          : fOutput(out) {}
-      void operator<<(const std::string& msg) { fOutput << msg; }
-      void Close() {}
-
-    private:
-      TStream& fOutput;
-    };
-
-    typedef Sink<std::ostream> SinkStream;
-
-  } // namespace sink
-} // namespace corsika::logging
-
-#endif

Framework/Particles/CMakeLists.txt

-
-add_custom_command (
-  OUTPUT  ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
-  COMMAND ${PROJECT_SOURCE_DIR}/Framework/Particles/pdxml_reader.py 
-          ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleData.xml
-          ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleClassNames.xml
-  DEPENDS pdxml_reader.py
-          ParticleData.xml
-          ParticleClassNames.xml
-  WORKING_DIRECTORY
-          ${PROJECT_BINARY_DIR}/Framework/Particles/
-  COMMENT "Read PYTHIA8 particle data and produce C++ source code GeneratedParticleProperties.inc"
-  VERBATIM
-  )
-
-# all public header files of library, includes automatic generated file(s)
-set (
-  PARTICLE_HEADERS
-  ParticleProperties.h
-  ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc # this one is auto-generated
-  )
-
-set (
-  PARTICLE_NAMESPACE
-  corsika/particles
-  )
-
-add_library (CORSIKAparticles INTERFACE)
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAparticles ${PARTICLE_NAMESPACE} ${PARTICLE_HEADERS})
-
-# ....................................................
-# since GeneratedParticleProperties.inc is an automatically produced file in the build directory,
-# create a symbolic link into the source tree, so that it can be found and edited more easily
-# this is not needed for the build to succeed! .......
-add_custom_command (
-  OUTPUT  ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc
-  COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/particles/GeneratedParticleProperties.inc ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc
-  COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc"
-  )
-add_custom_target (SourceDirLink DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc)
-add_dependencies (CORSIKAparticles SourceDirLink)
-# .....................................................
-
-
-target_include_directories (
-  CORSIKAparticles
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (
-  FILES
-  ${PARTICLE_HEADERS}
-  DESTINATION
-  include/${PARTICLE_NAMESPACE}
-  )
-
-# --------------------
-# code unit testing
-add_executable (
-  testParticles
-  testParticles.cc
-  ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
-  )
-
-target_link_libraries (
-  testParticles
-  CORSIKAparticles
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (
-  NAME
-  testParticles
-  COMMAND
-  testParticles
-  )
-

Framework/Particles/ParticleProperties.h

-/**
-   @file Particles.h
-
-   Interface to particle properties
- */
-
-#ifndef _include_ParticleProperties_h_
-#define _include_ParticleProperties_h_
-
-#include <array>
-#include <cstdint>
-#include <iostream>
-#include <type_traits>
-
-#include <corsika/units/PhysicalConstants.h>
-#include <corsika/units/PhysicalUnits.h>
-
-/**
- * @namespace particle
- *
- * The properties of all elementary particles is stored here. The data
- * is taken from the Pythia ParticleData.xml file.
- *
- */
-
-namespace corsika::particles {
-  enum class Code : int16_t;
-
-  using PDGCodeType = int16_t;
-  using CodeIntType = std::underlying_type<Code>::type;
-
-  // forward declarations to be used in GeneratedParticleProperties
-  int16_t constexpr GetElectricChargeNumber(Code const);
-  corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const);
-  corsika::units::si::MassType constexpr GetMass(Code const);
-  PDGCodeType constexpr GetPDG(Code const);
-  std::string const GetName(Code const);
-
-#include <corsika/particles/GeneratedParticleProperties.inc>
-
-  /*!
-   * returns mass of particle
-   */
-  corsika::units::si::MassType constexpr GetMass(Code const p) {
-    return masses[static_cast<CodeIntType const>(p)];
-  }
-
-  PDGCodeType constexpr GetPDG(Code const p) {
-    return pdg_codes[static_cast<CodeIntType const>(p)];
-  }
-
-  /*!
-   * returns electric charge of particle / (e/3).
-   */
-  int16_t constexpr GetElectricChargeNumber(Code const p) {
-    return electric_charges[static_cast<CodeIntType const>(p)];
-  }
-
-  corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
-    return GetElectricChargeNumber(p) * (corsika::units::si::constants::e / 3.);
-  }
-
-  std::string const GetName(Code const p) {
-    return names[static_cast<CodeIntType const>(p)];
-  }
-
-  namespace io {
-
-    std::ostream& operator<<(std::ostream& stream, Code const p) {
-      return stream << GetName(p);
-    }
-
-  } // 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["pdgID"])
-        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_camel(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}': result '{: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_camel(name) # the camel case names
-        
-        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 : int16_t {\n"
-              "  FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\"  \n") # identifier for eventual loops...
-    
-    
-    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 = {:d},\n" # identifier for eventual loops...
-               "}};").format(last_ngc_id + 1)
-               
-    if last_ngc_id > 0x7fff: # does not fit into int16_t
-        raise Exception("Integer overflow in internal particle code definition prevented!")
-    
-    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<corsika::units::si::MassType const, size> masses = {\n"    
-    for p in pythia_db.values():
-        string += "  {mass:f} * (1e9 * corsika::units::si::constants::eV / corsika::units::si::constants::cSquared), // {name:s}\n".format(mass = p['mass'], name = p['name'])              
-    string += "};\n\n"
-                   
-    # PDG code table
-    string += "static constexpr std::array<PDGCodeType 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_charges = {\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/c² \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 constexpr Code GetCode() { return Type; }\n"
-        string += "   static constexpr corsika::units::si::MassType GetMass() { return corsika::particles::GetMass(Type); }\n"
-        string += "   static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetElectricCharge(Type); }\n"
-        string += "   static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n"
-        string += "   static std::string const GetName() { return corsika::particles::GetName(Type); }\n"
-        string += "   static constexpr Code GetAntiParticle() { return AntiType; }\n"
-        string += "   static constexpr Code Type = Code::" + cname + ";\n"
-        string += "   static constexpr Code AntiType = Code::" + antiP + ";\n"
-        string += " private:\n"
-        string += "   static constexpr CodeIntType TypeIndex = static_cast<CodeIntType const>(Type);\n"
-        string += "};\n"
-
-    return string
-
-
-###############################################################
-# 
-# 
-
-def inc_start():
-    string = ('// generated by pdxml_reader.py\n'
-              '// MANUAL EDITS ON OWN RISK\n')
-    return string
-
-
-###############################################################
-# 
-# 
-
-def inc_end():
-    string = ""
-    return string
-
-
-###################################################################
-# 
-# Main function
-# 
-
-if __name__ == "__main__":
-
-    if len(sys.argv) != 3:
-        print("usage: {:s} <Pythia8.xml> <ClassNames.xml>".format(sys.argv[0]))
-        sys.exit(1)
-        
-    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::si;
-using namespace corsika::particles;
-
-TEST_CASE("ParticleProperties", "[Particles]") {
-
-  SECTION("Types") {
-    REQUIRE(Electron::GetCode() == Code::Electron);
-    REQUIRE(Positron::GetCode() == Code::Positron);
-    REQUIRE(Proton::GetCode() == Code::Proton);
-    REQUIRE(Neutron::GetCode() == Code::Neutron);
-    REQUIRE(Gamma::GetCode() == Code::Gamma);
-    REQUIRE(PiPlus::GetCode() == Code::PiPlus);
-  }
-
-  SECTION("Masses") {
-    REQUIRE(Electron::GetMass() / (511_keV / constants::cSquared) == Approx(1));
-    REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
-  }
-
-  SECTION("Charges") {
-    REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
-    REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
-    REQUIRE(GetElectricCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
-  }
-
-  SECTION("Names") {
-    REQUIRE(Electron::GetName() == "e-");
-    REQUIRE(PiMinus::GetName() == "pi-");
-  }
-
-  SECTION("PDG") {
-    REQUIRE(GetPDG(Code::PiPlus) == 211);
-    REQUIRE(GetPDG(Code::DPlus) == 411); 
-    REQUIRE(GetPDG(Code::NuMu) == 14); 
-    REQUIRE(GetPDG(Code::NuE) == 12); 
-    REQUIRE(GetPDG(Code::MuMinus) == 13); 
-  }
-  
-}

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 (
-  testProcessSequence
-  testProcessSequence.cc
-  )
-
-target_link_libraries (
-  testProcessSequence
-  CORSIKAprocesssequence
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (
-  NAME testProcessSequence
-  COMMAND testProcessSequence
-  )

Framework/ProcessSequence/ProcessSequence.h

-#ifndef _include_ProcessSequence_h_
-#define _include_ProcessSequence_h_
-
-#include <cmath>
-#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 {
-    derived& GetRef() { return static_cast<derived&>(*this); }
-    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 Particle, typename Trajectory, typename Stack>
-    inline void DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
-      A.DoContinuous(p, t, s);
-      B.DoContinuous(p, t, s);
-    } // add trajectory
-
-    template <typename D>
-    inline double MinStepLength(D& d) const {
-      return std::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 Particle, typename Stack>
-    void DoDiscrete(Particle& p, Stack& s) const {
-      A.DoDiscrete(p, s);
-      B.DoDiscrete(p, s);
-    }
-
-    /// TODO the const_cast is not nice, think about the constness here
-    inline void Init() const {
-      const_cast<T1*>(&A)->Init();
-      const_cast<T2*>(&B)->Init();
-    }
-  };
-
-  /// the + operator that assembles more BaseProcess objects into a ProcessSequence
-  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/ProcessSequence/testProcessSequence.cc

-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-#include <array>
-#include <iomanip>
-#include <iostream>
-
-#include <corsika/process/ProcessSequence.h>
-
-using namespace std;
-using namespace corsika::process;
-
-class Process1 : public BaseProcess<Process1> {
-public:
-  Process1() {}
-  void Init() {} // cout << "Process1::Init" << endl; }
-  template <typename D, typename T, typename S>
-  void DoContinuous(D& d, T& t, S& s) const {
-    for (int i = 0; i < 10; ++i) d.p[i] += 1 + i;
-  }
-};
-
-class Process2 : public BaseProcess<Process2> {
-public:
-  Process2() {}
-  void Init() {} // cout << "Process2::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline void DoContinuous(D& d, T& t, S& s) const {
-    for (int i = 0; i < 10; ++i) d.p[i] *= 0.7;
-  }
-};
-
-class Process3 : public BaseProcess<Process3> {
-public:
-  Process3() {}
-  void Init() {} // cout << "Process3::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline void DoContinuous(D& d, T& t, S& s) const {
-    for (int i = 0; i < 10; ++i) d.p[i] += 0.933;
-  }
-};
-
-class Process4 : public BaseProcess<Process4> {
-public:
-  Process4() {}
-  void Init() {} // cout << "Process4::Init" << endl; }
-  template <typename D, typename T, typename S>
-  inline void DoContinuous(D& d, T& t, S& s) const {
-    for (int i = 0; i < 10; ++i) d.p[i] /= 1.2;
-  }
-  // inline double MinStepLength(D& d) {
-  // void DoDiscrete(Particle& p, Stack& s) const {
-};
-
-struct DummyData {
-  double p[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-};
-struct DummyStack {};
-struct DummyTrajectory {};
-
-TEST_CASE("Cascade", "[Cascade]") {
-
-  SECTION("sectionTwo") {
-
-    Process1 m1;
-    Process2 m2;
-    Process3 m3;
-    Process4 m4;
-
-    const auto sequence = m1 + m2 + m3 + m4;
-
-    DummyData p;
-    DummyTrajectory t;
-    DummyStack s;
-
-    sequence.Init();
-
-    const int n = 100;
-    INFO("Running loop with n=" << n);
-    for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
-
-    for (int i = 0; i < 10; i++) { INFO("data[" << i << "]=" << p.p[i]); }
-  }
-
-  SECTION("sectionThree") {}
-}

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
-  ParticleBase.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}
-  )
-
-
-# code testing
-add_executable (testStackInterface testStackInterface.cc)
-target_link_libraries (testStackInterface CORSIKAstackinterface CORSIKAthirdparty) # for catch2
-add_test(NAME testStackInterface COMMAND testStackInterface)
-

Framework/StackInterface/ParticleBase.h

-#ifndef _include_particleBase_h_
-#define _include_particleBase_h_
-
-class StackData; // forward decl
-
-namespace corsika::stack {
-
-  //  template <typename> class PI;// : public ParticleBase<StackIteratorInterface> {
-  // template <typename, template <typename> typename> class Stack; // forward decl
-
-  /**
-   \class ParticleBase
-
-   The base class to define the readout of particle properties from a
-   particle stack. Every stack must implement this readout via the
-   ParticleBase class.
-  */
-
-  template <typename StackIterator>
-  class ParticleBase {
-
-    //    friend class Stack<StackData, PI>;  // for access to GetIterator
-  public:
-    ParticleBase() {}
-
-  private:
-    ParticleBase(ParticleBase&);
-    // ParticleBase& operation=(ParticleBase& p);
-
-  public:
-    /// delete this particle on the stack. The corresponding iterator
-    /// will be invalidated by this operation
-    void Delete() { GetIterator().GetStack().Delete(GetIterator()); }
-
-    //  protected: // todo should be proteced, but don't now how to 'friend Stack'
-    /// Function to provide CRTP access to inheriting class (type)
-    StackIterator& GetIterator() { return static_cast<StackIterator&>(*this); }
-    const StackIterator& GetIterator() const {
-      return static_cast<const StackIterator&>(*this);
-    }
-
-  protected:
-    /// access to underling stack data
-    auto& GetStackData() { return GetIterator().GetStackData(); }
-    const auto& GetStackData() const { return GetIterator().GetStackData(); }
-
-    /// return the index number of the underlying iterator object
-    int GetIndex() const { return GetIterator().GetIndex(); }
-  };
-
-}; // namespace corsika::stack
-
-#endif

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 {
-
-  template <typename>
-  class PI; // forward decl
-
-  /**
-     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 StackData, template <typename> typename PI>
-  class Stack : public StackData {
-
-  public:
-    typedef Stack<StackData, PI> StackType;
-    typedef StackIteratorInterface<StackData, PI> StackIterator;
-    typedef const StackIterator ConstStackIterator;
-    typedef typename StackIterator::ParticleInterfaceType ParticleType;
-    friend class StackIteratorInterface<StackData, PI>;
-
-  public:
-    using StackData::GetCapacity;
-    using StackData::GetSize;
-
-    using StackData::Clear;
-    using StackData::Copy;
-
-    using StackData::DecrementSize;
-    using StackData::IncrementSize;
-
-    using StackData::Init;
-
-  public:
-    /// these are functions required by std containers and std loops
-    StackIterator begin() { return StackIterator(*this, 0); }
-    StackIterator end() { return StackIterator(*this, GetSize()); }
-    StackIterator last() { return StackIterator(*this, GetSize() - 1); }
-
-    /// these are functions required by std containers and std loops
-    ConstStackIterator cbegin() const { return ConstStackIterator(*this, 0); }
-    ConstStackIterator cend() const { return ConstStackIterator(*this, GetSize()); }
-    ConstStackIterator clast() const { return ConstStackIterator(*this, GetSize() - 1); }
-
-    /// increase stack size, create new particle at end of stack
-    StackIterator NewParticle() {
-      IncrementSize();
-      return StackIterator(*this, GetSize() - 1);
-    }
-    /// delete this particle
-    void Delete(StackIterator& p) {
-      if (GetSize() == 0) { /*error*/
-      }
-      if (p.GetIndex() < GetSize() - 1) Copy(GetSize() - 1, p.GetIndex());
-      DeleteLast();
-      // p.SetInvalid();
-    }
-    void Delete(ParticleType& p) { Delete(p.GetIterator()); }
-    /// delete last particle on stack by decrementing stack size
-    void DeleteLast() { DecrementSize(); }
-    /// check if there are no further particles on stack
-    bool IsEmpty() { return GetSize() == 0; }
-    StackIterator GetNextParticle() { return last(); }
-
-  protected:
-    StackData& GetStackData() { return static_cast<StackData&>(*this); }
-    const StackData& GetStackData() const { return static_cast<const StackData&>(*this); }
-  };
-
-} // namespace corsika::stack
-
-#endif

Framework/StackInterface/StackIterator.h

-#ifndef _include_StackIterator_h__
-#define _include_StackIterator_h__
-
-#include <corsika/stack/ParticleBase.h>
-
-#include <iomanip>
-#include <iostream>
-
-class StackData; // forward decl
-
-namespace corsika::stack {
-
-  template <typename StackData, template <typename> typename ParticleInterface>
-  class Stack; // forward decl
-
-  /**
-     @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 StackData, template <typename> typename ParticleInterface>
-  class StackIteratorInterface
-      : public ParticleInterface<StackIteratorInterface<StackData, ParticleInterface> > {
-
-    typedef Stack<StackData, ParticleInterface> StackType;
-    typedef ParticleInterface<StackIteratorInterface<StackData, ParticleInterface> >
-        ParticleInterfaceType;
-
-    // friend class ParticleInterface<StackIterator<StackData>>; // to access GetStackData
-    friend class Stack<StackData, ParticleInterface>;  // for access to GetIndex
-    friend class ParticleBase<StackIteratorInterface>; // for access to GetStackData
-
-  private:
-    int fIndex = 0;
-    StackType* fData = 0; // todo is this problematic, when stacks are copied?
-
-  public:
-    // StackIterator() : fData(0), fIndex(0) { }
-    StackIteratorInterface(StackType& data, const int index)
-        : fData(&data)
-        , fIndex(index) {}
-
-  private:
-    StackIteratorInterface(const StackIteratorInterface& mit)
-        : fData(mit.fData)
-        , fIndex(mit.fIndex) {}
-    
-  public:
-    StackIteratorInterface& operator=(const StackIteratorInterface& mit) {
-      fData = mit.fData;
-      fIndex = mit.fIndex;
-      return *this;
-    }
-
-  public:
-    StackIteratorInterface& operator++() {
-      ++fIndex;
-      return *this;
-    }
-    StackIteratorInterface operator++(int) {
-      StackIteratorInterface tmp(*this);
-      ++fIndex;
-      return tmp;
-    }
-    bool operator==(const StackIteratorInterface& rhs) { return fIndex == rhs.fIndex; }
-    bool operator!=(const StackIteratorInterface& rhs) { return fIndex != rhs.fIndex; }
-
-    ParticleInterfaceType& operator*() {
-      return static_cast<ParticleInterfaceType&>(*this);
-    }
-    const ParticleInterfaceType& operator*() const {
-      return static_cast<const ParticleInterfaceType&>(*this);
-    }
-
-  protected:
-    int GetIndex() const { return fIndex; }
-    StackType& GetStack() { return *fData; }
-    const StackType& GetStack() const { return *fData; }
-    StackData& GetStackData() { return fData->GetStackData(); }
-    const StackData& GetStackData() const { return fData->GetStackData(); }
-
-  }; // end class StackIterator
-
-} // namespace corsika::stack
-
-#endif

Framework/StackInterface/testStackInterface.cc

-#include <corsika/stack/Stack.h>
-
-#include <iomanip>
-#include <iostream>
-#include <vector>
-
-#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::stack;
-using namespace std;
-
-// definition of stack-data object
-class StackOneData {
-
-public:
-  // these functions are needed for the Stack interface
-  void Init() {}
-  void Clear() { fData.clear(); }
-  int GetSize() const { return fData.size(); }
-  int GetCapacity() const { return fData.size(); }
-  void Copy(const int i1, const int i2) { fData[i2] = fData[i1]; }
-  void Swap(const int i1, const int i2) {
-    double tmp0 = fData[i1];
-    fData[i1] = fData[i2];
-    fData[i2] = tmp0;
-  }
-
-  // custom data access function
-  void SetData(const int i, const double v) { fData[i] = v; }
-  double GetData(const int i) const { return fData[i]; }
-
-protected:
-  // these functions are also needed by the Stack interface
-  void IncrementSize() { fData.push_back(0.); }
-  void DecrementSize() {
-    if (fData.size() > 0) { fData.pop_back(); }
-  }
-
-  // custom private data section
-private:
-  std::vector<double> fData;
-};
-
-// defintion of a stack-readout object, the iteractor dereference
-// operator will deliver access to these function
-template <typename StackIteratorInterface>
-class ParticleInterface : public ParticleBase<StackIteratorInterface> {
-  //  using ParticleBase<StackIteratorInterface>::Delete;
-  using ParticleBase<StackIteratorInterface>::GetStackData;
-  using ParticleBase<StackIteratorInterface>::GetIndex;
-
-public:
-  void SetData(const double v) { GetStackData().SetData(GetIndex(), v); }
-  double GetData() const { return GetStackData().GetData(GetIndex()); }
-};
-
-TEST_CASE("Stack", "[Stack]") {
-
-  SECTION("StackInterface") {
-
-    // construct a valid Stack object
-    typedef Stack<StackOneData, ParticleInterface> StackTest;
-    StackTest s;
-    s.Init();
-    s.Clear();
-    s.IncrementSize();
-    s.Copy(0, 0);
-    s.Swap(0, 0);
-    s.GetCapacity();
-    REQUIRE(s.GetSize() == 1);
-    s.DecrementSize();
-    REQUIRE(s.GetSize() == 0);
-  }
-
-  SECTION("write") {
-
-    // construct a valid Stack object
-    typedef Stack<StackOneData, ParticleInterface> StackTest;
-    StackTest s;
-  }
-
-  SECTION("read") {
-
-    typedef Stack<StackOneData, ParticleInterface> StackTest;
-    StackTest s;
-    s.NewParticle().SetData(9.9);
-    cout << "kk" << endl;
-    double v = 0;
-    for (auto& p : s) {
-      cout << typeid(p).name() << endl;
-      v += p.GetData();
-    }
-    cout << "k222k" << endl;
-
-    REQUIRE(v == 9.9);
-  }
-
-  SECTION("delete_stack") {
-
-    typedef Stack<StackOneData, ParticleInterface> StackTest;
-    StackTest s;
-    auto p = s.NewParticle();
-    p.SetData(9.9);
-    s.Delete(p);
-  }
-
-  SECTION("delete_particle") {
-
-    typedef Stack<StackOneData, ParticleInterface> StackTest;
-    StackTest s;
-    auto p = s.NewParticle();
-    p.SetData(9.9);
-    p.Delete();
-  }
-}

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::si::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.6021766208e-19L) * joule};
-
-  // elementary charge
-  constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-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};
-  constexpr auto cSquared = c * c;
-
-  // unified atomic mass unit
-  constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
-
-  // etc.
-
-} // namespace corsika::units::si::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::si::constants::eV))
-
-      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
-
-    } // namespace literals
-  }   // namespace units
-} // namespace phys
-
-namespace corsika::units::si {
-  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>;
-  using MassType = phys::units::quantity<phys::units::mass_d, double>;
-
-} // end namespace corsika::units::si
-
-// 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::si;
-
-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);
-    REQUIRE(E2 / 1_GeV == Approx(40));
-
-    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);
-  }
-}