Merge branch 'master' into 89-consistent-units-in-cascade-step

AUTHORS

 Maximilian Reininghaus, maximilian.reininghaus@kit.edu
 Ralf Ulrich, ralf.ulrich@kit.edu
+Felix Riehn, friehn@lip.pt

Documentation/Examples/cascade_example.cc

@@ -55,7 +55,7 @@ static int fEmCount;
 static EnergyType fInvEnergy;
 static int fInvCount;
 
-class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
+class ProcessEMCut : public corsika::process::ContinuousProcess<ProcessEMCut> {
 public:
   ProcessEMCut() {}
   template <typename Particle>
@@ -131,7 +131,24 @@ public:
   }
 
   template <typename Particle, typename Stack>
-  EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+  EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
+    cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
+    const Code pid = p.GetPID();
+    if (isEmParticle(pid)) {
+      cout << "removing em. particle..." << endl;
+      fEmEnergy += p.GetEnergy();
+      fEmCount += 1;
+      p.Delete();
+    } else if (isInvisible(pid)) {
+      cout << "removing inv. particle..." << endl;
+      fInvEnergy += p.GetEnergy();
+      fInvCount += 1;
+      p.Delete();
+    } else if (isBelowEnergyCut(p)) {
+      cout << "removing low en. particle..." << endl;
+      fEnergy += p.GetEnergy();
+      p.Delete();
+    }
     // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
     // cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
     // cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
@@ -153,27 +170,6 @@ public:
     return EProcessReturn::eOk;
   }
 
-  template <typename Particle, typename Stack>
-  void DoDiscrete(Particle& p, Stack&) const {
-    cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
-    const Code pid = p.GetPID();
-    if (isEmParticle(pid)) {
-      cout << "removing em. particle..." << endl;
-      fEmEnergy += p.GetEnergy();
-      fEmCount += 1;
-      p.Delete();
-    } else if (isInvisible(pid)) {
-      cout << "removing inv. particle..." << endl;
-      fInvEnergy += p.GetEnergy();
-      fInvCount += 1;
-      p.Delete();
-    } else if (isBelowEnergyCut(p)) {
-      cout << "removing low en. particle..." << endl;
-      fEnergy += p.GetEnergy();
-      p.Delete();
-    }
-  }
-
   void Init() {
     fEmEnergy = 0. * 1_GeV;
     fEmCount = 0;
@@ -223,15 +219,22 @@ int main() {
 
   universe.AddChild(std::move(theMedium));
 
-  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
   tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
 
+<<<<<<< HEAD
   corsika::process::sibyll::Interaction /*<setup::Stack>,setup::Trajectory>*/ p1;
   corsika::process::sibyll::Decay p2;
   ProcessEMCut p3;
   const auto sequence = /*p0 +*/ p1 + p2 + p3;
+=======
+  corsika::process::sibyll::Interaction sibyll;
+  corsika::process::sibyll::Decay decay;
+  ProcessEMCut cut;
+  const auto sequence = /*p0 +*/ sibyll + decay + cut;
+>>>>>>> master
   setup::Stack stack;
 
   corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
@@ -256,7 +259,7 @@ int main() {
        << endl;
   cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
        << std::endl;
-  p3.ShowResults();
+  cut.ShowResults();
   cout << "total energy (GeV): "
-       << (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;
+       << (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy()) / 1_GeV << endl;
 }

Documentation/Examples/geometry_example.cc

@@ -26,7 +26,7 @@ using namespace corsika::units::si;
 int main() {
   // define the root coordinate system
   geometry::CoordinateSystem& root =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
   // another CS defined by a translation relative to the root CS
   CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});

Documentation/Examples/helix_example.cc

@@ -24,7 +24,7 @@ using namespace corsika::units::si;
 
 int main() {
   geometry::CoordinateSystem& root =
-      geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
   Point const r0(root, {0_m, 0_m, 0_m});
   auto const omegaC = 2 * M_PI * 1_Hz;

Environment/CMakeLists.txt

@@ -44,6 +44,7 @@ add_executable (testEnvironment testEnvironment.cc)
 
 target_link_libraries (
   testEnvironment
+  CORSIKAsetup
   CORSIKAenvironment
   CORSIKAthirdparty # for catch2
   )

Environment/Environment.h

@@ -35,9 +35,8 @@ namespace corsika::environment {
   class Environment {
   public:
     Environment()
-        : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
-                                .GetRootCS()}
-        , fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+        : fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem()},        
+        fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
               std::make_unique<Universe>(fCoordinateSystem))) {}
 
     using IEnvironmentModel = corsika::setup::IEnvironmentModel;

Framework/Geometry/RootCoordinateSystem.h

@@ -21,8 +21,8 @@ namespace corsika::geometry {
     RootCoordinateSystem() {}
 
   public:
-    corsika::geometry::CoordinateSystem& GetRootCS() { return fRootCS; }
-    const corsika::geometry::CoordinateSystem& GetRootCS() const { return fRootCS; }
+    corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() { return fRootCS; }
+    const corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() const { return fRootCS; }
 
   private:
     corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT

Framework/Geometry/testGeometry.cc

@@ -29,7 +29,7 @@ using namespace corsika::units::si;
 double constexpr absMargin = 1.0e-8;
 
 TEST_CASE("transformations between CoordinateSystems") {
-  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
   REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
               .isApprox(EigenTransform::Identity()));
@@ -128,7 +128,7 @@ TEST_CASE("transformations between CoordinateSystems") {
 }
 
 TEST_CASE("Sphere") {
-  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
   Point center(rootCS, {0_m, 3_m, 4_m});
   Sphere sphere(center, 5_m);
 
@@ -147,7 +147,7 @@ TEST_CASE("Sphere") {
 }
 
 TEST_CASE("Trajectories") {
-  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
   Point r0(rootCS, {0_m, 0_m, 0_m});
 
   SECTION("Line") {

Framework/Particles/CMakeLists.txt

 
 add_custom_command (
   OUTPUT  ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
-          ${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl
+          ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
   COMMAND ${PROJECT_SOURCE_DIR}/Framework/Particles/pdxml_reader.py 
           ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleData.xml
+          ${PROJECT_SOURCE_DIR}/Framework/Particles/NuclearData.xml
           ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleClassNames.xml
   DEPENDS pdxml_reader.py
           ParticleData.xml
+          NuclearData.xml
           ParticleClassNames.xml
   WORKING_DIRECTORY
           ${PROJECT_BINARY_DIR}/Framework/Particles/

Framework/Particles/NuclearData.xml

+<chapter name="Nuclear Data">
+
+<particle id="1000010010" name="hydrogen" A="1" Z="1" >
+</particle> 
+
+<particle id="1000010020" name="deuterium" A="2" Z="1" >
+</particle> 
+
+<particle id="1000010030" name="tritium" A="3" Z="1" >
+</particle> 
+
+<particle id="1000020040" name="helium" A="4" Z="2" >
+</particle> 
+
+<particle id="1000060120" name="carbon" A="12" Z="6" >
+</particle> 
+
+<particle id="1000070140" name="nitrogen" A="14" Z="7" >
+</particle> 
+
+<particle id="1000080160" name="oxygen" A="16" Z="8" >
+</particle> 
+
+<particle id="1000180400" name="argon" A="40" Z="18" >
+</particle> 
+
+</chapter>

Framework/Particles/ParticleProperties.h

@@ -40,7 +40,7 @@ namespace corsika::particles {
 
   enum class Code : int16_t;
 
-  using PDGCodeType = int16_t;
+  using PDGCodeType = int32_t;
   using CodeIntType = std::underlying_type<Code>::type;
 
   // forward declarations to be used in GeneratedParticleProperties
@@ -51,24 +51,28 @@ namespace corsika::particles {
   constexpr std::string const& GetName(Code const);
   corsika::units::si::TimeType constexpr GetLifetime(Code const);
 
+  bool constexpr IsNucleus(Code const);
+  int constexpr GetNucleusA(Code const);
+  int constexpr GetNucleusZ(Code const);
+  
 #include <corsika/particles/GeneratedParticleProperties.inc>
 
   /*!
    * returns mass of particle
    */
   corsika::units::hep::MassType constexpr GetMass(Code const p) {
-    return masses[static_cast<CodeIntType const>(p)];
+    return detail::masses[static_cast<CodeIntType const>(p)];
   }
 
   PDGCodeType constexpr GetPDG(Code const p) {
-    return pdg_codes[static_cast<CodeIntType const>(p)];
+    return detail::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)];
+    return detail::electric_charges[static_cast<CodeIntType const>(p)];
   }
 
   corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
@@ -76,13 +80,26 @@ namespace corsika::particles {
   }
 
   constexpr std::string const& GetName(Code const p) {
-    return names[static_cast<CodeIntType const>(p)];
+    return detail::names[static_cast<CodeIntType const>(p)];
   }
 
   corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
-    return lifetime[static_cast<CodeIntType const>(p)];
+    return detail::lifetime[static_cast<CodeIntType const>(p)] * corsika::units::si::second;
   }
 
+  bool constexpr IsNucleus(Code const p) {
+    return detail::isNucleus[static_cast<CodeIntType const>(p)];
+  }
+
+  int constexpr GetNucleusA(Code const p) {
+      return detail::nucleusA[static_cast<CodeIntType const>(p)];
+  }
+
+  int constexpr GetNucleusZ(Code const p) {
+    return detail::nucleusZ[static_cast<CodeIntType const>(p)];
+  }
+
+
   namespace io {
 
     std::ostream& operator<<(std::ostream& stream, Code const p);

Framework/Particles/pdxml_reader.py

@@ -5,21 +5,23 @@ import xml.etree.ElementTree as ET
 from collections import OrderedDict
 import pickle
 
+GeVfm = 0.19732696312541853
+c_speed_of_light = 29.9792458e10  # mm / s
+# for nuclear masses
+mneutron = 0.9395654133 # GeV
+mproton = 0.9382720813 # GeV
+
 
 ##############################################################
 # 
 # reading xml input data, return line by line particle data
 # 
-
-def parse(filename):
+def parsePythia(filename):
     tree = ET.parse(filename)
     root = tree.getroot()
-
-    GeVfm = 0.19732696312541853
-    c_speed_of_light = 29.9792458e10  # mm / s
     
     for particle in root.iter("particle"):
-        name = particle.attrib["name"]
+        name = particle.attrib["name"]        
         antiName = "Unknown"
         if ("antiName" in particle.attrib):
             antiName = particle.attrib["antiName"]
@@ -47,12 +49,39 @@ def parse(filename):
             yield (-pdg_id, antiName, mass, -electric_charge, name, ctau/c_speed_of_light)
 
 
-            
 ##############################################################
 # 
-# returns dict with particle codes and class names
+# reading xml input data, return line by line particle data
 # 
+def parseNuclei(filename):
+    tree = ET.parse(filename)
+    root = tree.getroot()
+        
+    for particle in root.iter("particle"):
+        name = particle.attrib["name"]        
+        antiName = "Unknown"
+        if ("antiName" in particle.attrib):
+            antiName = particle.attrib["antiName"]
+        pdg_id = int(particle.attrib["id"])
+        A = int(particle.attrib["A"])
+        Z = int(particle.attrib["Z"])
+        # mass in GeV
+        if ("mass" in particle.attrib):
+            mass = particle.attrib["mass"] 
+        else:
+            mass = (A-Z)*mneutron + Z*mproton
+
+        electric_charge = Z*3  # in units of e/3        
+        ctau = float('Inf')
+        
+        yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light, A, Z)
+                
 
+        
+##############################################################
+# 
+# returns dict with particle codes and class names
+# 
 def class_names(filename):
     tree = ET.parse(filename)
     root = tree.getroot()
@@ -164,14 +193,12 @@ def c_identifier_camel(name):
 ##########################################################
 # 
 # returns dict containing all data from pythia-xml input
-# 
+#     
+def read_pythia_db(filename, particle_db, classnames):    
     
-def build_pythia_db(filename, classnames):    
-    particle_db = OrderedDict()
+    counter = itertools.count(len(particle_db))
     
-    counter = itertools.count(0)
-    
-    for (pdg, name, mass, electric_charge, antiName, lifetime) in parse(filename):
+    for (pdg, name, mass, electric_charge, antiName, lifetime) in parsePythia(filename):
 
         c_id = "Unknown"
         if pdg in classnames:
@@ -186,7 +213,41 @@ def build_pythia_db(filename, classnames):
             "mass" : mass, # in GeV
             "electric_charge" : electric_charge, # in e/3
             "lifetime" : lifetime,
-            "ngc_code" : next(counter)
+            "ngc_code" : next(counter),
+            "isNucleus" : False,
+        }
+    
+    return particle_db
+    
+
+
+##########################################################
+# 
+# returns dict containing all data from pythia-xml input
+#     
+def read_nuclei_db(filename, particle_db, classnames):    
+    
+    counter = itertools.count(len(particle_db))
+
+    for (pdg, name, mass, electric_charge, antiName, lifetime, A, Z) in parseNuclei(filename):
+        
+        c_id = "Unknown"
+        if pdg in classnames:
+            c_id = classnames[pdg]
+        else:
+            c_id = c_identifier_camel(name)
+
+        particle_db[c_id] = {
+            "name" : name,
+            "antiName" : antiName,
+            "pdg" : pdg,
+            "mass" : mass, # in GeV
+            "electric_charge" : electric_charge, # in e/3
+            "lifetime" : lifetime,
+            "ngc_code" : next(counter),
+            "A" : A,
+            "Z" : Z,
+            "isNucleus" : True,
         }
     
     return particle_db
@@ -198,14 +259,13 @@ def build_pythia_db(filename, classnames):
 # 
 # return string with enum of all internal particle codes
 # 
-
-def gen_internal_enum(pythia_db):
+def gen_internal_enum(particle_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']
+    for k in filter(lambda k: "ngc_code" in particle_db[k], particle_db):
+        last_ngc_id = particle_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...
@@ -223,52 +283,83 @@ def gen_internal_enum(pythia_db):
 # 
 # return string with all data arrays 
 # 
-
-def gen_properties(pythia_db):
+def gen_properties(particle_db):
 
     # number of particles, size of tables
-    string = "static constexpr std::size_t size = {size:d};\n".format(size = len(pythia_db))
+    string = "static constexpr std::size_t size = {size:d};\n".format(size = len(particle_db))
     string += "\n"
     
     # particle masses table
     string += "static constexpr std::array<corsika::units::hep::MassType const, size> masses = {\n"    
-    for p in pythia_db.values():
-        string += "  {mass:f} * (1e9 * corsika::units::si::constants::eV), // {name:s}\n".format(mass = p['mass'], name = p['name'])              
+    for p in particle_db.values():
+        string += "  {mass:e} * (1e9 * corsika::units::si::constants::eV), // {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():
+    for p in particle_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():
+    for p in particle_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():
+    for p in particle_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"
+    #    string += "static constexpr std::array<size, size> anti_particle = {\n"
+    #    for p in particle_db.values():
+    #        string += "  {anti:d},\n".format(charge = p['anti_particle'])
+    #    string += "};\n"
 
     # lifetime
-    string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n"
-    for p in pythia_db.values():
+    #string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n"
+    string += "static constexpr std::array<double const, size> lifetime = {\n"
+    for p in particle_db.values():
         if p['lifetime'] == float("Inf") :
-            string += "  std::numeric_limits<double>::infinity() * corsika::units::si::second, \n"
+            string += "  std::numeric_limits<double>::infinity(), \n" # * corsika::units::si::second, \n"
         else :
-            string += "  {tau:f} * corsika::units::si::second, \n".format(tau = p['lifetime'])
+            string += "  {tau:e}, \n".format(tau = p['lifetime'])
+            #string += "  {tau:e} * corsika::units::si::second, \n".format(tau = p['lifetime'])
+    string += "};\n"
+
+    
+    ### nuclear data ###
+    
+    # is nucleus flag
+    string += "static constexpr std::array<bool, size> isNucleus = {\n"
+    for p in particle_db.values():
+        value = 'false'
+        if p['isNucleus']:
+            value = 'true'
+        string += "  {val},\n".format(val = value)
+    string += "};\n"
+
+    # nucleus mass number A
+    string += "static constexpr std::array<int, size> nucleusA = {\n"
+    for p in particle_db.values():
+        A = 0
+        if p['isNucleus']:
+            A = p['A']
+        string += "  {val},\n".format(val = A)
     string += "};\n"
     
+    # nucleus charge number Z
+    string += "static constexpr std::array<int, size> nucleusZ = {\n"
+    for p in particle_db.values():
+        Z = 0
+        if p['isNucleus']:
+            Z = p['Z']
+        string += "  {val},\n".format(val = Z)
+    string += "};\n"
+
     return string
 
 
@@ -278,27 +369,29 @@ def gen_properties(pythia_db):
 # 
 # return string with a list of classes for all particles
 # 
-
-def gen_classes(pythia_db):
+def gen_classes(particle_db):
 
     string = "// list of C++ classes to access particle properties"
     
-    for cname in pythia_db:
+    for cname in particle_db:
 
         antiP = 'Unknown'
-        for cname_anti in pythia_db:
-            if (pythia_db[cname_anti]['name'] == pythia_db[cname]['antiName']):
+        for cname_anti in particle_db:
+            if (particle_db[cname_anti]['name'] == particle_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 += " *  - pdg=" + str(particle_db[cname]['pdg']) +"\n"
+        string += " *  - mass=" + str(particle_db[cname]['mass']) + " GeV \n"
+        string += " *  - charge= " + str(particle_db[cname]['electric_charge']/3) + " \n"
         string += " *  - name=" + str(cname) + "\n"
         string += " *  - anti=" + str(antiP) + "\n"
+        if (particle_db[cname]['isNucleus']):
+            string += " *  - nuclear A=" + str(particle_db[cname]['A']) + "\n"
+            string += " *  - nuclear Z=" + str(particle_db[cname]['Z']) + "\n"        
         string += "*/\n\n"
         string += "class " + cname + " {\n"
         string += "  public:\n"
@@ -308,6 +401,9 @@ def gen_classes(pythia_db):
         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 bool IsNucleus() { return corsika::particles::IsNucleus(Type); }\n"
+        string += "   static constexpr int GetNucleusA() { return corsika::particles::GetNucleusA(Type); }\n"
+        string += "   static constexpr int GetNucleusZ() { return corsika::particles::GetNucleusZ(Type); }\n"
         string += "   static constexpr Code Type = Code::" + cname + ";\n"
         string += "   static constexpr Code AntiType = Code::" + antiP + ";\n"
         string += " private:\n"
@@ -320,17 +416,30 @@ def gen_classes(pythia_db):
 ###############################################################
 # 
 # 
-
 def inc_start():
     string = ('// generated by pdxml_reader.py\n'
-              '// MANUAL EDITS ON OWN RISK\n')
+              '// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. \n')
     return string
 
 
 ###############################################################
 # 
 # 
+def detail_start():
+    string = ('namespace detail {\n\n')
+    return string
+
+
+###############################################################
+# 
+# 
+def detail_end():
+    string = "\n}//end namespace detail\n"
+    return string
 
+###############################################################
+# 
+# 
 def inc_end():
     string = ""
     return string
@@ -338,36 +447,38 @@ def inc_end():
 
 ###################################################################
 # 
-# Serialize pythia_db into file 
+# Serialize particle_db into file 
 # 
+def serialize_particle_db(particle_db, file):
+    pickle.dump(particle_db, file)
 
-def serialize_pythia_db(pythia_db, file):
-    pickle.dump(pythia_db, file)
-
+    
 ###################################################################
 # 
 # Main function
 # 
-
 if __name__ == "__main__":
-
-    if len(sys.argv) != 3:
-        print("usage: {:s} <Pythia8.xml> <ClassNames.xml>".format(sys.argv[0]), file=sys.stderr)
+    
+    if len(sys.argv) != 4:
+        print("usage: {:s} <Pythia8.xml> <Nuclei.xml> <ClassNames.xml>".format(sys.argv[0]), file=sys.stderr)
         sys.exit(1)
         
-    print("\n       pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n")
-
-    names = class_names(sys.argv[2])
-    pythia_db = build_pythia_db(sys.argv[1], names)
-
-    print (str(pythia_db))
+    print("\n       pdxml_reader.py: Automatically produce particle-properties from input files\n")
+    
+    names = class_names(sys.argv[3])
+    particle_db = OrderedDict()
+    read_pythia_db(sys.argv[1], particle_db, names)
+    read_nuclei_db(sys.argv[2], particle_db, names)
     
     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_start(), file=f)
+        print(gen_internal_enum(particle_db), file=f)
+        print(detail_start(), file=f)
+        print(gen_properties(particle_db), file=f)
+        print(detail_end(), file=f) 
+        print(gen_classes(particle_db), file=f)
         print(inc_end(), file=f) 
     
-    with open("pythia_db.pkl", "wb") as f:
-        serialize_pythia_db(pythia_db, f)
+    with open("particle_db.pkl", "wb") as f:
+        serialize_particle_db(particle_db, f)
+    

Framework/Particles/testParticles.cc

@@ -58,11 +58,28 @@ TEST_CASE("ParticleProperties", "[Particles]") {
     REQUIRE(GetLifetime(Code::Electron) ==
             std::numeric_limits<double>::infinity() * corsika::units::si::second);
     REQUIRE(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma));
-    // REQUIRE(GetLifetime(Code::RhoPlus)/corsika::units::si::second ==
-    // (Approx(4.414566727909413e-24).epsilon(1e-3)));
-    // REQUIRE(GetLifetime(Code::SigmaMinusBar)/corsika::units::si::second ==
-    // (Approx(8.018880848563575e-11).epsilon(1e-5)));
-    // REQUIRE(GetLifetime(Code::MuPlus)/corsika::units::si::second ==
-    // (Approx(2.1970332555864364e-06).epsilon(1e-5)));
+    REQUIRE(GetLifetime(Code::RhoPlus)/corsika::units::si::second ==
+	    (Approx(4.414566727909413e-24).epsilon(1e-3)));
+    REQUIRE(GetLifetime(Code::SigmaMinusBar)/corsika::units::si::second ==
+	    (Approx(8.018880848563575e-11).epsilon(1e-5)));
+    REQUIRE(GetLifetime(Code::MuPlus)/corsika::units::si::second ==
+	    (Approx(2.1970332555864364e-06).epsilon(1e-5)));
   }
+
+  SECTION("Nuclei") {
+    REQUIRE(IsNucleus(Code::Gamma) == false);
+    REQUIRE(IsNucleus(Code::Argon) == true);
+    REQUIRE(IsNucleus(Code::Proton) == false);
+    REQUIRE(IsNucleus(Code::Hydrogen) == true);
+    REQUIRE(Argon::IsNucleus() == true);
+    REQUIRE(EtaC::IsNucleus() == false);
+
+    REQUIRE(GetNucleusA(Code::Hydrogen) == 1);
+    REQUIRE(GetNucleusA(Code::Tritium) == 3);
+    REQUIRE(Hydrogen::GetNucleusZ() == 1);
+    REQUIRE(Tritium::GetNucleusA() == 3);
+
+
+  }
+
 }

Processes/NullModel/CMakeLists.txt

@@ -29,6 +29,8 @@ set_target_properties (
 target_link_libraries (
   ProcessNullModel
   CORSIKAunits
+  CORSIKAgeometry
+  CORSIKAsetup
   )
 
 target_include_directories (
@@ -52,6 +54,8 @@ add_executable (testNullModel testNullModel.cc)
 
 target_link_libraries (
   testNullModel
+  ProcessNullModel
+  CORSIKAsetup
   CORSIKAgeometry
   CORSIKAunits
   CORSIKAthirdparty # for catch2

Processes/NullModel/NullModel.cc

@@ -11,14 +11,40 @@
 
 #include <corsika/process/null_model/NullModel.h>
 
+
+#include <corsika/logging/Logger.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <iostream>
+#include <limits>
+using namespace std;
+
+using namespace corsika;
+using namespace corsika::units::si;
 using namespace corsika::process::null_model;
 
-NullModel::NullModel() {}
+template <typename Stack>
+NullModel<Stack>::NullModel() {}
+
+template <typename Stack>
+NullModel<Stack>::~NullModel() {}
+
+template <typename Stack>
+process::EProcessReturn NullModel<Stack>::DoContinuous(Particle&, setup::Trajectory&,
+						       Stack& ) const {
+  return EProcessReturn::eOk;
+}
+
+template <typename Stack>
+double NullModel<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
+  return std::numeric_limits<double>::infinity();
+}
 
-NullModel::~NullModel() {}
+template <typename Stack>
+void NullModel<Stack>::Init() {
+}
 
-void NullModel::init() {}
+#include <corsika/setup/SetupStack.h>
 
-void NullModel::run() {}
+template class process::null_model::NullModel<setup::Stack>;
 
-double NullModel::GetStepLength() { return 0; }

Processes/NullModel/NullModel.h

@@ -12,19 +12,26 @@
 #ifndef _Physics_NullModel_NullModel_h_
 #define _Physics_NullModel_NullModel_h_
 
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupTrajectory.h>
+
 namespace corsika::process {
 
   namespace null_model {
 
+    template <typename Stack>
     class NullModel {
 
+      typedef typename Stack::ParticleType Particle;
+
     public:
       NullModel();
       ~NullModel();
 
-      void init();
-      void run();
-      double GetStepLength();
+      void Init();
+      EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
+      double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
+
     };
 
   } // namespace null_model

Processes/NullModel/testNullModel.cc

@@ -13,11 +13,41 @@
                           // cpp file
 #include <catch2/catch.hpp>
 
+#include <corsika/process/null_model/NullModel.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using namespace corsika::units::si;
+using namespace corsika::process::null_model;
+using namespace corsika;
+
 TEST_CASE("NullModel", "[processes]") {
 
-  SECTION("bla") {}
+  auto const& cs = geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();    
+  geometry::Point const origin(cs, {0_m, 0_m, 0_m});    
+  geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+								    0_m / second, 1_m / second);
+  geometry::Line line(origin, v);
+  geometry::Trajectory<geometry::Line> track(line, 10_s);
+
+  setup::Stack stack;
+  auto particle = stack.NewParticle();
+  
+  SECTION("interface") {
+
+    NullModel<setup::Stack> model;
+
+    model.Init();
+    [[maybe_unused]] const process::EProcessReturn ret = model.DoContinuous(particle, track, stack);
+    [[maybe_unused]] const double length = model.MaxStepLength(particle, track);
 
-  SECTION("blubb") {}
+    
+  }
 }

Processes/Sibyll/CMakeLists.txt

 add_custom_command (
   OUTPUT  ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
   COMMAND ${PROJECT_SOURCE_DIR}/Processes/Sibyll/code_generator.py 
-          ${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl
+          ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
           ${PROJECT_SOURCE_DIR}/Processes/Sibyll/sibyll_codes.dat
   DEPENDS code_generator.py
           sibyll_codes.dat          
-          ${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl
+          ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
   WORKING_DIRECTORY
           ${PROJECT_BINARY_DIR}/Processes/Sibyll/
   COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA"
@@ -94,6 +94,9 @@ add_executable (testSibyll
 
 target_link_libraries (
   testSibyll
+  ProcessSibyll
+  CORSIKAsetup
+  CORSIKArandom
   CORSIKAgeometry
   CORSIKAunits
   CORSIKAthirdparty # for catch2

Processes/Sibyll/Decay.h

@@ -5,6 +5,7 @@
 #include <corsika/process/sibyll/ParticleConversion.h>
 #include <corsika/process/sibyll/SibStack.h>
 #include <corsika/setup/SetupTrajectory.h>
+#include <corsika/setup/SetupStack.h>
 
 #include <corsika/particles/ParticleProperties.h>
 

Processes/Sibyll/Interaction.h

@@ -20,14 +20,7 @@ using namespace corsika::units::si;
 
 namespace corsika::process::sibyll {
 
-  // template <typename Stack, typename Track>
-  // template <typename Stack>
-  class Interaction
-      : public corsika::process::InteractionProcess<Interaction> { // <Stack,Track>> {
-
-    // typedef typename Stack::ParticleType Particle;
-    // typedef typename corsika::setup::Stack::ParticleType Particle;
-    // typedef corsika::setup::Trajectory Track;
+  class Interaction : public corsika::process::InteractionProcess<Interaction> { 
 
   public:
     Interaction() {}
@@ -51,9 +44,8 @@ namespace corsika::process::sibyll {
 
     template <typename Particle, typename Track>
     corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) const {
-
       // coordinate system, get global frame of reference
-      CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+      CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
       const particles::Code corsikaBeamId = p.GetPID();