IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 2155 deletions
Framework/Logging/testLogging.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/logging/Logger.h>
#include <catch2/catch.hpp>
/*
this does nothing yet
*/
TEST_CASE("Logging", "[Logging]") {
SECTION("sectionOne") {}
SECTION("sectionTwo") {}
SECTION("sectionThree") {}
}
Framework/Particles/CMakeLists.txt deleted 100644 → 0
View file @ 3deb039c
set(Python_ADDITIONAL_VERSIONS 3)
find_package(PythonInterp 3 REQUIRED)
add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
${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/
COMMENT "Read PYTHIA8 particle data and produce C++ source code GeneratedParticleProperties.inc"
VERBATIM
)
set (
PARTICLE_SOURCES
ParticleProperties.cc
)
# 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 STATIC ${PARTICLE_SOURCES})
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 ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc)
add_dependencies (CORSIKAparticles SourceDirLink)
# .....................................................
target_link_libraries (
CORSIKAparticles
PUBLIC
CORSIKAunits
)
set_target_properties (
CORSIKAparticles
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
PUBLIC_HEADER "${PARTICLE_HEADERS}"
)
target_include_directories (
CORSIKAparticles
PUBLIC
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
FILES
${PARTICLE_HEADERS}
DESTINATION
include/${PARTICLE_NAMESPACE}
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testParticles
SOURCES
testParticles.cc
${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
)
target_link_libraries (
testParticles
CORSIKAparticles
CORSIKAunits
CORSIKAtesting
)
Framework/Particles/ParticleProperties.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/particles/ParticleProperties.h>
#include <iostream>
namespace corsika::particles {
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p) {
return stream << corsika::particles::GetName(p);
}
Code ConvertFromPDG(PDGCode p) {
static_assert(detail::conversionArray.size() % 2 == 1);
// this will fail, for the strange case where the maxPDG is negative...
unsigned int constexpr maxPDG{(detail::conversionArray.size() - 1) >> 1};
auto k = static_cast<PDGCodeType>(p);
if ((unsigned int)abs(k) <= maxPDG) {
return detail::conversionArray[k + maxPDG];
} else {
return detail::conversionMap.at(p);
}
}
} // namespace corsika::particles
Framework/Particles/ParticleProperties.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
/**
@file Particles.h
Interface to particle properties
*/
#ifndef _include_corsika_particles_ParticleProperties_h_
#define _include_corsika_particles_ParticleProperties_h_
#include <array>
#include <cstdint>
#include <iosfwd>
#include <type_traits>
#include <corsika/units/PhysicalUnits.h>
/**
*
* The properties of all elementary particles is stored here. The data
* are taken from the Pythia ParticleData.xml file.
*
*/
namespace corsika::particles {
/**
* @enum Code
* The Code enum is the actual place to define CORSIKA 8 particle codes.
*/
enum class Code : int16_t;
enum class PDGCode : int32_t;
using CodeIntType = std::underlying_type<Code>::type;
using PDGCodeType = std::underlying_type<PDGCode>::type;
// forward declarations to be used in GeneratedParticleProperties
int16_t constexpr GetChargeNumber(Code const);
corsika::units::si::ElectricChargeType constexpr GetCharge(Code const);
corsika::units::si::HEPMassType constexpr GetMass(Code const);
PDGCode constexpr GetPDG(Code const);
constexpr std::string const& GetName(Code const);
corsika::units::si::TimeType constexpr GetLifetime(Code const);
bool constexpr IsNucleus(Code const);
bool constexpr IsHadron(Code const);
bool constexpr IsEM(Code const);
bool constexpr IsMuon(Code const);
bool constexpr IsNeutrino(Code const);
int constexpr GetNucleusA(Code const);
int constexpr GetNucleusZ(Code const);
#include <corsika/particles/GeneratedParticleProperties.inc>
/*!
* returns mass of particle in natural units
*/
corsika::units::si::HEPMassType constexpr GetMass(Code const p) {
if (p == Code::Nucleus)
throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
return detail::masses[static_cast<CodeIntType>(p)];
}
/*!
* returns PDG id
*/
PDGCode constexpr GetPDG(Code const p) {
return detail::pdg_codes[static_cast<CodeIntType>(p)];
}
/*!
* returns electric charge number of particle return 1 for a proton.
*/
int16_t constexpr GetChargeNumber(Code const p) {
if (p == Code::Nucleus)
throw std::runtime_error(
"Cannot GetChargeNumber() of particle::Nucleus -> unspecified");
// electric_charges stores charges in units of (e/3), e.g. 3 for a proton
return detail::electric_charges[static_cast<CodeIntType>(p)] / 3;
}
/*!
* returns electric charge of particle, e.g. return 1.602e-19_C for a proton.
*/
corsika::units::si::ElectricChargeType constexpr GetCharge(Code const p) {
if (p == Code::Nucleus)
throw std::runtime_error("Cannot GetCharge() of particle::Nucleus -> unspecified");
return GetChargeNumber(p) * corsika::units::constants::e;
}
constexpr std::string const& GetName(Code const p) {
return detail::names[static_cast<CodeIntType>(p)];
}
corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
return detail::lifetime[static_cast<CodeIntType>(p)] * corsika::units::si::second;
}
bool constexpr IsHadron(Code const p) {
return detail::isHadron[static_cast<CodeIntType>(p)];
}
bool constexpr IsEM(Code c) {
return c == Code::Electron || c == Code::Positron || c == Code::Gamma;
}
bool constexpr IsMuon(Code c) { return c == Code::MuPlus || c == Code::MuMinus; }
bool constexpr IsNeutrino(Code c) {
return c == Code::NuE || c == Code::NuMu || c == Code::NuTau || c == Code::NuEBar ||
c == Code::NuMuBar || c == Code::NuTauBar;
}
int constexpr GetNucleusA(Code const p) {
if (p == Code::Nucleus) {
throw std::runtime_error("GetNucleusA(Code::Nucleus) is impossible!");
}
return detail::nucleusA[static_cast<CodeIntType>(p)];
}
int constexpr GetNucleusZ(Code const p) {
if (p == Code::Nucleus) {
throw std::runtime_error("GetNucleusZ(Code::Nucleus) is impossible!");
}
return detail::nucleusZ[static_cast<CodeIntType>(p)];
}
bool constexpr IsNucleus(Code const p) {
return (p == Code::Nucleus) || (GetNucleusA(p) != 0);
}
/**
* the output operator for humand-readable particle codes
**/
std::ostream& operator<<(std::ostream&, corsika::particles::Code);
Code ConvertFromPDG(PDGCode);
/**
* Get mass of nucleus
**/
corsika::units::si::HEPMassType constexpr GetNucleusMass(const int vA, const int vZ) {
return Proton::GetMass() * vZ + (vA - vZ) * Neutron::GetMass();
}
} // namespace corsika::particles
#endif
Framework/Particles/particles.dox deleted 100644 → 0
View file @ 3deb039c
/**
@page Particles Particle properties
The properties of all particles are saved in static and flat
arrays. There is a enum corsika::particles::Code to identify each
particles, and each individual particles has its own static class,
which can be used to retrieve its physical properties.
See namespace corsika::particles for all classes.
*/
\ No newline at end of file
Framework/Particles/pdxml_reader.py deleted 100755 → 0
View file @ 3deb039c
#!/usr/bin/env python3
import sys, math, itertools, re, csv, pprint
import xml.etree.ElementTree as ET
from collections import OrderedDict
import pickle
import io
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 parsePythia(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"])
mass = float(particle.attrib["m0"]) # GeV
electric_charge = int(particle.attrib["chargeType"]) # in units of e/3
ctau = 0.
if pdg_id in (11, 12, 14, 16, 22, 2212): # these are the stable particles !
ctau = float('Inf')
elif 'tau0' in particle.attrib:
ctau = float(particle.attrib['tau0']) # mm / c
elif 'mWidth' in particle.attrib:
ctau = GeVfm / float(particle.attrib['mWidth']) * 1e-15 * 1000.0 # mm / s
elif pdg_id in (0, 423, 433, 4312, 4322, 5112, 5222): # those are certainly not stable....
ctau = 0.
else:
print ("missing lifetime: " + str(pdg_id) + " " + str(name))
sys.exit(1)
yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light)
# TODO: read decay channels from child elements
if "antiName" in particle.attrib:
yield (-pdg_id, antiName, mass, -electric_charge, name, ctau/c_speed_of_light)
##############################################################
#
# 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()
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 "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:] + '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 numbers
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 read_pythia_db(filename, particle_db, classnames):
counter = itertools.count(len(particle_db))
for (pdg, name, mass, electric_charge, antiName, lifetime) in parsePythia(filename):
c_id = "Unknown"
if pdg in classnames:
c_id = classnames[pdg]
else:
c_id = c_identifier_camel(name) # the camel case names
hadron = False
if abs(pdg) > 100:
hadron = True
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),
"isNucleus" : False,
"isHadron" : hadron,
}
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,
"isHadron" : True,
}
return particle_db
###############################################################
#
# build conversion table PDG -> ngc
#
def gen_conversion_PDG_ngc(particle_db):
# todo: find a optimum value, think about cache miss with array vs lookup time with map
P_MAX = 500 # the maximum PDG code that is still filled into the table
conversionDict = dict()
conversionTable = [None] * (2*P_MAX + 1)
for cId, p in particle_db.items():
pdg = p['pdg']
if abs(pdg) < P_MAX:
if conversionTable[pdg + P_MAX]:
raise Exception("table entry already occupied")
else:
conversionTable[pdg + P_MAX] = cId
else:
if pdg in conversionDict.keys():
raise Exception(f"map entry {pdg} already occupied")
else:
conversionDict[pdg] = cId
output = io.StringIO()
def oprint(*args, **kwargs):
print(*args, **kwargs, file=output)
oprint(f"static std::array<Code, {len(conversionTable)}> constexpr conversionArray {{")
for ngc in conversionTable:
oprint(" Code::{0},".format(ngc if ngc else "Unknown"))
oprint("};")
oprint()
oprint("static std::map<PDGCode, Code> const conversionMap {")
for ngc in conversionDict.values():
oprint(f" {{PDGCode::{ngc}, Code::{ngc}}},")
oprint("};")
oprint()
return output.getvalue()
###############################################################
#
# return string with enum of all internal particle codes
#
def gen_internal_enum(particle_db):
string = ("enum class Code : CodeIntType {\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 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...
"}};").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 enum of all PDG particle codes
#
def gen_pdg_enum(particle_db):
string = "enum class PDGCode : PDGCodeType {\n"
for cId in particle_db:
pdgCode = particle_db[cId]['pdg']
string += " {key:s} = {code:d},\n".format(key = cId, code = pdgCode)
string += " };\n"
return string
###############################################################
#
# return string with all data arrays
#
def gen_properties(particle_db):
# number of particles, size of tables
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::si::HEPMassType const, size> masses = {\n"
for p in particle_db.values():
string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n"
# PDG code table
string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
for p in particle_db.keys():
string += f" PDGCode::{p},\n"
string += "};\n"
# name string table
string += "static const std::array<std::string const, size> names = {\n"
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 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 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"
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(), \n" # * corsika::units::si::second, \n"
else :
string += " {tau:e}, \n".format(tau = p['lifetime'])
#string += " {tau:e} * corsika::units::si::second, \n".format(tau = p['lifetime'])
string += "};\n"
# is Hadron flag
string += "static constexpr std::array<bool, size> isHadron = {\n"
for p in particle_db.values():
value = 'false'
if p['isHadron']:
value = 'true'
string += " {val},\n".format(val = value)
string += "};\n"
### nuclear data ###
# nucleus mass number A
string += "static constexpr std::array<int16_t, size> nucleusA = {\n"
for p in particle_db.values():
A = p.get('A', 0)
string += " {val},\n".format(val = A)
string += "};\n"
# nucleus charge number Z
string += "static constexpr std::array<int16_t, size> nucleusZ = {\n"
for p in particle_db.values():
Z = p.get('Z', 0)
string += " {val},\n".format(val = Z)
string += "};\n"
return string
###############################################################
#
# return string with a list of classes for all particles
#
def gen_classes(particle_db):
string = "// list of C++ classes to access particle properties"
for cname in particle_db:
antiP = 'Unknown'
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(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"
string += " static constexpr Code GetCode() { return Type; }\n"
string += " static constexpr corsika::units::si::HEPMassType GetMass() { return corsika::particles::GetMass(Type); }\n"
string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetCharge(Type); }\n"
string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetChargeNumber(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 int16_t GetNucleusA() { return corsika::particles::GetNucleusA(Type); }\n"
string += " static constexpr int16_t 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"
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. 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
###################################################################
#
# Serialize particle_db into file
#
def serialize_particle_db(particle_db, file):
pickle.dump(particle_db, file)
###################################################################
#
# Main function
#
if __name__ == "__main__":
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 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(particle_db), file=f)
print(gen_pdg_enum(particle_db), file=f)
print(detail_start(), file=f)
print(gen_properties(particle_db), file=f)
print(gen_conversion_PDG_ngc(particle_db), file=f)
print(detail_end(), file=f)
print(gen_classes(particle_db), file=f)
print(inc_end(), file=f)
with open("particle_db.pkl", "wb") as f:
serialize_particle_db(particle_db, f)
Framework/Particles/testParticles.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
using namespace corsika::units;
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) == Approx(1));
REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
REQUIRE((Proton::GetMass() + Neutron::GetMass()) /
corsika::units::constants::nucleonMass ==
Approx(2));
}
SECTION("Charges") {
REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
REQUIRE(GetCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
}
SECTION("Names") {
REQUIRE(Electron::GetName() == "e-");
REQUIRE(PiMinus::GetName() == "pi-");
REQUIRE(Nucleus::GetName() == "nucleus");
REQUIRE(Iron::GetName() == "iron");
}
SECTION("PDG") {
REQUIRE(GetPDG(Code::PiPlus) == PDGCode::PiPlus);
REQUIRE(GetPDG(Code::DPlus) == PDGCode::DPlus);
REQUIRE(GetPDG(Code::NuMu) == PDGCode::NuMu);
REQUIRE(GetPDG(Code::NuE) == PDGCode::NuE);
REQUIRE(GetPDG(Code::MuMinus) == PDGCode::MuMinus);
REQUIRE(static_cast<int>(GetPDG(Code::PiPlus)) == 211);
REQUIRE(static_cast<int>(GetPDG(Code::DPlus)) == 411);
REQUIRE(static_cast<int>(GetPDG(Code::NuMu)) == 14);
REQUIRE(static_cast<int>(GetPDG(Code::NuEBar)) == -12);
REQUIRE(static_cast<int>(GetPDG(Code::MuMinus)) == 13);
}
SECTION("Conversion PDG -> internal") {
REQUIRE(ConvertFromPDG(PDGCode::KStarMinus) == Code::KStarMinus);
REQUIRE(ConvertFromPDG(PDGCode::MuPlus) == Code::MuPlus);
REQUIRE(ConvertFromPDG(PDGCode::SigmaStarCMinusBar) == Code::SigmaStarCMinusBar);
}
SECTION("Lifetimes") {
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)));
}
SECTION("Particle groups: electromagnetic") {
REQUIRE(IsEM(Code::Gamma));
REQUIRE(IsEM(Code::Electron));
REQUIRE_FALSE(IsEM(Code::MuPlus));
REQUIRE_FALSE(IsEM(Code::NuE));
REQUIRE_FALSE(IsEM(Code::Proton));
REQUIRE_FALSE(IsEM(Code::PiPlus));
REQUIRE_FALSE(IsEM(Code::Oxygen));
}
SECTION("Particle groups: hadrons") {
REQUIRE_FALSE(IsHadron(Code::Gamma));
REQUIRE_FALSE(IsHadron(Code::Electron));
REQUIRE_FALSE(IsHadron(Code::MuPlus));
REQUIRE_FALSE(IsHadron(Code::NuE));
REQUIRE(IsHadron(Code::Proton));
REQUIRE(IsHadron(Code::PiPlus));
REQUIRE(IsHadron(Code::Oxygen));
REQUIRE(IsHadron(Code::Nucleus));
}
SECTION("Particle groups: muons") {
REQUIRE_FALSE(IsMuon(Code::Gamma));
REQUIRE_FALSE(IsMuon(Code::Electron));
REQUIRE(IsMuon(Code::MuPlus));
REQUIRE_FALSE(IsMuon(Code::NuE));
REQUIRE_FALSE(IsMuon(Code::Proton));
REQUIRE_FALSE(IsMuon(Code::PiPlus));
REQUIRE_FALSE(IsMuon(Code::Oxygen));
}
SECTION("Particle groups: neutrinos") {
REQUIRE_FALSE(IsNeutrino(Code::Gamma));
REQUIRE_FALSE(IsNeutrino(Code::Electron));
REQUIRE_FALSE(IsNeutrino(Code::MuPlus));
REQUIRE(IsNeutrino(Code::NuE));
REQUIRE_FALSE(IsNeutrino(Code::Proton));
REQUIRE_FALSE(IsNeutrino(Code::PiPlus));
REQUIRE_FALSE(IsNeutrino(Code::Oxygen));
}
SECTION("Nuclei") {
REQUIRE_FALSE(IsNucleus(Code::Gamma));
REQUIRE(IsNucleus(Code::Argon));
REQUIRE_FALSE(IsNucleus(Code::Proton));
REQUIRE(IsNucleus(Code::Hydrogen));
REQUIRE(Argon::IsNucleus());
REQUIRE_FALSE(EtaC::IsNucleus());
REQUIRE(GetNucleusA(Code::Hydrogen) == 1);
REQUIRE(GetNucleusA(Code::Tritium) == 3);
REQUIRE(Hydrogen::GetNucleusZ() == 1);
REQUIRE(Tritium::GetNucleusA() == 3);
REQUIRE_THROWS(GetNucleusA(Code::Nucleus));
REQUIRE_THROWS(GetNucleusZ(Code::Nucleus));
}
}
Framework/ProcessSequence/BaseProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_baseprocess_h_
#define _include_corsika_baseprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <type_traits>
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 {
private:
BaseProcess() {}
friend Derived;
public:
Derived& GetRef() { return static_cast<Derived&>(*this); }
const Derived& GetRef() const { return static_cast<const Derived&>(*this); }
};
// overwrite the default trait class, to mark BaseProcess<T> as useful process
template <class T>
std::true_type is_process_impl(const BaseProcess<T>* impl);
} // namespace corsika::process
#endif
Framework/ProcessSequence/BoundaryCrossingProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_processes_BoundaryCrossingProcess_h_
#define _include_corsika_processes_BoundaryCrossingProcess_h_
#include <corsika/environment/Environment.h>
#include <corsika/process/ProcessReturn.h>
namespace corsika::process {
template <typename TDerived>
struct BoundaryCrossingProcess {
auto& GetRef() { return static_cast<TDerived&>(*this); }
auto const& GetRef() const { return static_cast<const TDerived&>(*this); }
/**
* This method is called when a particle crosses the boundary between the nodes
* \p from and \p to.
*/
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle&, VTNType const& from, VTNType const& to);
};
template <class T>
std::true_type is_process_impl(BoundaryCrossingProcess<T> const* impl);
} // namespace corsika::process
#endif
Framework/ProcessSequence/CMakeLists.txt deleted 100644 → 0
View file @ 3deb039c
add_library (CORSIKAprocesssequence INTERFACE)
#namespace of library->location of header files
set (
CORSIKAprocesssequence_NAMESPACE
corsika/process
)
#header files of this library
set (
CORSIKAprocesssequence_HEADERS
BaseProcess.h
BoundaryCrossingProcess.h
ContinuousProcess.h
SecondariesProcess.h
InteractionProcess.h
StackProcess.h
DecayProcess.h
ProcessSequence.h
ProcessReturn.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}
)
target_link_libraries (
CORSIKAprocesssequence
INTERFACE
CORSIKAenvironment
)
#-- -- -- -- -- -- -- --
#code unit testing
CORSIKA_ADD_TEST (testProcessSequence)
target_link_libraries (
testProcessSequence
ProcessSwitch
CORSIKAsetup
CORSIKAgeometry
CORSIKAprocesssequence
CORSIKAtesting
)
Framework/ProcessSequence/ContinuousProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_continuousprocess_h_
#define _include_corsika_continuousprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class ContinuousProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type ContinuousProcess<T>
*/
template <typename derived>
struct ContinuousProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
// here starts the interface part
// -> enforce derived to implement DoContinuous...
template <typename Particle, typename Track>
EProcessReturn DoContinuous(Particle&, Track const&) const;
// -> enforce derived to implement MaxStepLength...
template <typename Particle, typename Track>
units::si::LengthType MaxStepLength(Particle const& p, Track const& track) const;
};
// overwrite the default trait class, to mark BaseProcess<T> as useful process
template <class T>
std::true_type is_process_impl(const ContinuousProcess<T>* impl);
} // namespace corsika::process
#endif
Framework/ProcessSequence/DecayProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_decayprocess_h_
#define _include_corsika_decayprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class DecayProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type DecayProcess<T>
*/
template <typename derived>
struct DecayProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/// here starts the interface-definition part
// -> enforce derived to implement DoDecay...
template <typename Particle>
EProcessReturn DoDecay(Particle&);
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle& p);
template <typename Particle>
corsika::units::si::InverseTimeType GetInverseLifetime(Particle& vP) {
return 1. / GetRef().GetLifetime(vP);
}
};
// overwrite the default trait class, to mark DecayProcess<T> as useful process
template <class T>
std::true_type is_process_impl(const DecayProcess<T>* impl);
} // namespace corsika::process
#endif
Framework/ProcessSequence/InteractionProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_interactionprocess_h_
#define _include_corsika_interactionprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class InteractionProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type InteractionProcess<T>
*/
template <typename derived>
struct InteractionProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/// here starts the interface-definition part
// -> enforce derived to implement DoInteraction...
template <typename Particle>
EProcessReturn DoInteraction(Particle&);
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
template <typename TParticle>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
return 1. / GetRef().GetInteractionLength(p);
}
};
// overwrite the default trait class, to mark BaseProcess<T> as useful process
template <class T>
std::true_type is_process_impl(const InteractionProcess<T>* impl);
} // namespace corsika::process
#endif
Framework/ProcessSequence/ProcessReturn.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_ProcessReturn_h_
#define _include_ProcessReturn_h_
namespace corsika::process {
/**
since in a process sequence many status updates can accumulate
for a single particle, this enum should define only bit-flags
that can be accumulated easily with "|="
*/
enum class EProcessReturn : int {
eOk = (1 << 0),
eParticleAbsorbed = (1 << 2),
eInteracted = (1 << 3),
eDecayed = (1 << 4),
};
inline EProcessReturn operator|(EProcessReturn a, EProcessReturn b) {
return static_cast<EProcessReturn>(static_cast<int>(a) | static_cast<int>(b));
}
inline EProcessReturn& operator|=(EProcessReturn& a, EProcessReturn b) {
return a = a | b;
}
inline bool operator==(EProcessReturn a, EProcessReturn b) {
return (static_cast<int>(a) & static_cast<int>(b)) != 0;
}
} // namespace corsika::process
#endif
Framework/ProcessSequence/ProcessSequence.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_ProcessSequence_h_
#define _include_ProcessSequence_h_
#include <corsika/process/BaseProcess.h>
#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/StackProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
#include <limits>
#include <type_traits>
namespace corsika::process {
/**
\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
*/
// define a marker (trait class) to tag any class that qualifies as "Process" for the
// "ProcessSequence"
std::false_type is_process_impl(...);
template <class T>
using is_process = decltype(is_process_impl(std::declval<T*>()));
// this is a marker to track which BaseProcess is also a ProcessSequence
template <typename T>
struct is_process_sequence : std::false_type {};
template <typename T>
bool constexpr is_process_sequence_v = is_process_sequence<T>::value;
namespace switch_process {
template <typename A, typename B>
class SwitchProcess; // fwd-decl.
}
// to detect SwitchProcesses inside the ProcessSequence
template <typename T>
struct is_switch_process : std::false_type {};
template <typename T>
bool constexpr is_switch_process_v = is_switch_process<T>::value;
template <typename A, typename B>
struct is_switch_process<switch_process::SwitchProcess<A, B>> : std::true_type {};
/**
T1 and T2 are both references if possible (lvalue), otherwise
(rvalue) they are just classes. This allows us to handle both,
rvalue as well as lvalue Processes in the ProcessSequence.
*/
template <typename T1, typename T2>
class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2>> {
using T1type = typename std::decay<T1>::type;
using T2type = typename std::decay<T2>::type;
static bool constexpr t1ProcSeq = is_process_sequence_v<T1type>;
static bool constexpr t2ProcSeq = is_process_sequence_v<T2type>;
static bool constexpr t1SwitchProc = is_switch_process_v<T1type>;
static bool constexpr t2SwitchProc = is_switch_process_v<T2type>;
public:
T1 A; // this is a reference, if possible
T2 B; // this is a reference, if possible
ProcessSequence(T1 in_A, T2 in_B)
: A(in_A)
, B(in_B) {}
// example for a trait-based call:
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from,
VTNType const& to) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T1type>, T1type> ||
t1ProcSeq) {
ret |= A.DoBoundaryCrossing(p, from, to);
}
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T2type>, T2type> ||
t2ProcSeq) {
ret |= B.DoBoundaryCrossing(p, from, to);
}
return ret;
}
template <typename TParticle, typename TTrack>
EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
ret |= A.DoContinuous(vP, vT);
}
if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
ret |= B.DoContinuous(vP, vT);
}
return ret;
}
template <typename TSecondaries>
EProcessReturn DoSecondaries(TSecondaries& vS) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of_v<SecondariesProcess<T1type>, T1type> || t1ProcSeq) {
ret |= A.DoSecondaries(vS);
}
if constexpr (std::is_base_of_v<SecondariesProcess<T2type>, T2type> || t2ProcSeq) {
ret |= B.DoSecondaries(vS);
}
return ret;
}
/**
The processes of type StackProcess do have an internal counter,
so they can be exectuted only each N steps. Often these are
"maintenacne processes" that do not need to run after each
single step of the simulations. In the CheckStep function it is
tested if either A or B are StackProcess and if they are due
for execution.
*/
bool CheckStep() {
bool ret = false;
if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
ret |= A.CheckStep();
}
if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
ret |= B.CheckStep();
}
return ret;
}
/**
Execute the StackProcess-es in the ProcessSequence
*/
template <typename TStack>
EProcessReturn DoStack(TStack& vS) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
if (A.CheckStep()) { ret |= A.DoStack(vS); }
}
if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
if (B.CheckStep()) { ret |= B.DoStack(vS); }
}
return ret;
}
template <typename TParticle, typename TTrack>
corsika::units::si::LengthType MaxStepLength(TParticle& vP, TTrack& vTrack) {
corsika::units::si::LengthType
max_length = // if no other process in the sequence implements it
std::numeric_limits<double>::infinity() * corsika::units::si::meter;
if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack);
max_length = std::min(max_length, len);
}
if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack);
max_length = std::min(max_length, len);
}
return max_length;
}
template <typename TParticle>
corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP) {
return 1. / GetInverseInteractionLength(vP);
}
template <typename TParticle>
corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
TParticle& vP) {
return GetInverseInteractionLength(vP);
}
template <typename TParticle>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP) {
using namespace corsika::units::si;
InverseGrammageType tot = 0 * meter * meter / gram;
if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type> || t1ProcSeq ||
t1SwitchProc) {
tot += A.GetInverseInteractionLength(vP);
}
if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type> || t2ProcSeq ||
t2SwitchProc) {
tot += B.GetInverseInteractionLength(vP);
}
return tot;
}
template <typename TParticle, typename TSecondaries>
EProcessReturn SelectInteraction(
TParticle& vP, TSecondaries& vS,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if constexpr (t1ProcSeq || t1SwitchProc) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += A.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
A.DoInteraction(vS);
return EProcessReturn::eInteracted;
}
} // end branch A
if constexpr (t2ProcSeq || t2SwitchProc) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
B.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += B.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
B.DoInteraction(vS);
return EProcessReturn::eInteracted;
}
} // end branch A
return EProcessReturn::eOk;
}
template <typename TParticle>
corsika::units::si::TimeType GetTotalLifetime(TParticle& p) {
return 1. / GetInverseLifetime(p);
}
template <typename TParticle>
corsika::units::si::InverseTimeType GetTotalInverseLifetime(TParticle& p) {
return GetInverseLifetime(p);
}
template <typename TParticle>
corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& p) {
using namespace corsika::units::si;
corsika::units::si::InverseTimeType tot = 0 / second;
if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type> || t1ProcSeq) {
tot += A.GetInverseLifetime(p);
}
if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type> || t2ProcSeq) {
tot += B.GetInverseLifetime(p);
}
return tot;
}
// select decay process
template <typename TParticle, typename TSecondaries>
EProcessReturn SelectDecay(
TParticle& vP, TSecondaries& vS,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
corsika::units::si::InverseTimeType& decay_inv_count) {
if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret = A.SelectDecay(vP, vS, decay_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += A.GetInverseLifetime(vP);
// check if we should execute THIS process and then EXIT
if (decay_select < decay_inv_count) { // more pedagogical: rndm_select <
// decay_inv_count / decay_inv_tot
A.DoDecay(vS);
return EProcessReturn::eDecayed;
}
} // end branch A
if constexpr (t2ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret = B.SelectDecay(vP, vS, decay_select, decay_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_count += B.GetInverseLifetime(vP);
// check if we should execute THIS process and then EXIT
if (decay_select < decay_inv_count) {
B.DoDecay(vS);
return EProcessReturn::eDecayed;
}
} // end branch B
return EProcessReturn::eOk;
}
void Init() {
A.Init();
B.Init();
}
};
/// the << operator assembles many BaseProcess, ContinuousProcess, and
/// Interaction/DecayProcess objects into a ProcessSequence, all combinatorics
/// must be allowed, this is why we define a macro to define all
/// combinations here:
template <
typename P1, typename P2,
typename std::enable_if<is_process<typename std::decay<P1>::type>::value &&
is_process<typename std::decay<P2>::type>::value>::type...>
inline auto operator<<(P1&& vA, P2&& vB) -> ProcessSequence<P1, P2> {
return ProcessSequence<P1, P2>(vA.GetRef(), vB.GetRef());
}
/// marker to identify objectas ProcessSequence
template <typename A, typename B>
struct is_process_sequence<corsika::process::ProcessSequence<A, B>> : std::true_type {};
} // namespace corsika::process
#endif
Framework/ProcessSequence/ProcessSignature.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_process_processsignature_h_
#define _include_process_processsignature_h_
#define FORCE_SIGNATURE(nameTrait, nameMethod, signatureMethod) \
template <typename U> \
class nameTrait { \
private: \
template <typename T, T> \
struct helper; \
template <typename T> \
static std::uint8_t check(helper<signatureMethod, &nameMethod>*); \
template <typename T> \
static std::uint16_t check(...); \
\
public: \
static constexpr bool value = sizeof(check<U>(0)) == sizeof(std::uint8_t); \
}
// FORCE_SIGNATURE(thisMustBeDefined, T::thisMustBeDefined, int(*)(void));
#endif
Framework/ProcessSequence/SecondariesProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_secondariesprocess_h_
#define _include_corsika_secondariesprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class SecondariesProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type SecondariesProcess<T>
*/
template <typename derived>
struct SecondariesProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/// here starts the interface-definition part
// -> enforce derived to implement DoSecondaries...
template <typename TSecondaries>
inline EProcessReturn DoSecondaries(TSecondaries&);
};
// overwrite the default trait class, to mark BaseProcess<T> as useful process
template <class T>
std::true_type is_process_impl(const SecondariesProcess<T>* impl);
} // namespace corsika::process
#endif
Framework/ProcessSequence/StackProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_stackprocess_h_
#define _include_corsika_stackprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class StackProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type StackProcess<T>
*/
template <typename derived>
class StackProcess {
public:
StackProcess() = delete;
StackProcess(const unsigned int nStep)
: fNStep(nStep) {}
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/// here starts the interface-definition part
// -> enforce derived to implement DoStack...
template <typename TStack>
inline EProcessReturn DoStack(TStack&);
int GetStep() const { return fIStep; }
bool CheckStep() { return !((++fIStep) % fNStep); }
private:
/**
@name The number of "steps" during the cascade processing after
which this StackProcess is going to be executed. The logic is
"fIStep modulo fNStep"
@{
*/
unsigned int fNStep = 0;
unsigned long int fIStep = 0;
//! @}
};
// overwrite the default trait class, to mark BaseProcess<T> as useful process
template <class T>
std::true_type is_process_impl(const StackProcess<T>* impl);
} // namespace corsika::process
#endif
Framework/ProcessSequence/testProcessSequence.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <catch2/catch.hpp>
#include <array>
#include <iomanip>
#include <iostream>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/switch_process/SwitchProcess.h>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::process;
using namespace std;
static const int nData = 10;
int globalCount = 0;
class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
int fV = 0;
public:
ContinuousProcess1(const int v)
: fV(v) {}
void Init() {
cout << "ContinuousProcess1::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess1::DoContinuous" << endl;
for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
}
};
class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
int fV = 0;
public:
ContinuousProcess2(const int v)
: fV(v) {}
void Init() {
cout << "ContinuousProcess2::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess2::DoContinuous" << endl;
for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
return EProcessReturn::eOk;
}
};
class Process1 : public InteractionProcess<Process1> {
public:
Process1(const int v)
: fV(v) {}
void Init() {
cout << "Process1::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename D, typename S>
inline EProcessReturn DoInteraction(D& d, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
return EProcessReturn::eOk;
}
// private:
int fV;
};
class Process2 : public InteractionProcess<Process2> {
int fV = 0;
public:
Process2(const int v)
: fV(v) {}
void Init() {
cout << "Process2::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename Particle>
inline EProcessReturn DoInteraction(Particle&) const {
cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process2::GetInteractionLength" << endl;
return 3_g / (1_cm * 1_cm);
}
};
class Process3 : public InteractionProcess<Process3> {
int fV = 0;
public:
Process3(const int v)
: fV(v) {}
void Init() {
cout << "Process3::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename Particle>
inline EProcessReturn DoInteraction(Particle&) const {
cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process3::GetInteractionLength" << endl;
return 1_g / (1_cm * 1_cm);
}
};
class Process4 : public BaseProcess<Process4> {
int fV = 0;
public:
Process4(const int v)
: fV(v) {}
void Init() {
cout << "Process4::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; }
return EProcessReturn::eOk;
}
// inline double MinStepLength(D& d) {
template <typename Particle>
EProcessReturn DoInteraction(Particle&) const {
return EProcessReturn::eOk;
}
};
class Decay1 : public DecayProcess<Decay1> {
int fV = 0;
public:
Decay1(const int v)
: fV(v) {}
void Init() {
cout << "Decay1::Init" << endl;
assert(globalCount == fV);
globalCount++;
}
template <typename Particle>
TimeType GetLifetime(Particle&) const {
return 1_s;
}
template <typename Particle>
EProcessReturn DoDecay(Particle&) const {
return EProcessReturn::eOk;
}
};
class Stack1 : public StackProcess<Stack1> {
int fCount = 0;
public:
Stack1(const int n)
: StackProcess(n) {}
template <typename TStack>
EProcessReturn DoStack(TStack&) {
fCount++;
return EProcessReturn::eOk;
}
int GetCount() const { return fCount; }
};
struct DummyStack {};
struct DummyData {
double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
struct DummyTrajectory {};
TEST_CASE("Process Sequence", "[Process Sequence]") {
SECTION("Check init order") {
Process1 m1(0);
Process2 m2(1);
Process3 m3(2);
Process4 m4(3);
auto sequence = m1 << m2 << m3 << m4;
globalCount = 0;
sequence.Init();
// REQUIRE_NOTHROW( (sequence.Init()) );
// const auto sequence_wrong = m3 + m2 + m1 + m4;
// globalCount = 0;
// sequence_wrong.Init();
// REQUIRE_THROWS(sequence_wrong.Init());
}
SECTION("interaction length") {
ContinuousProcess1 cp1(0);
Process2 m2(1);
Process3 m3(2);
DummyData particle;
auto sequence2 = cp1 << m2 << m3;
GrammageType const tot = sequence2.GetTotalInteractionLength(particle);
InverseGrammageType const tot_inv =
sequence2.GetTotalInverseInteractionLength(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
SECTION("lifetime") {
ContinuousProcess1 cp1(0);
Process2 m2(1);
Process3 m3(2);
Decay1 d3(2);
DummyData particle;
auto sequence2 = cp1 << m2 << m3 << d3;
TimeType const tot = sequence2.GetTotalLifetime(particle);
InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
SECTION("sectionTwo") {
ContinuousProcess1 cp1(0);
ContinuousProcess2 cp2(3);
Process2 m2(1);
Process3 m3(2);
auto sequence2 = cp1 << m2 << m3 << cp2;
DummyData particle;
DummyTrajectory track;
cout << "-->init sequence2" << endl;
globalCount = 0;
sequence2.Init();
cout << "-->docont" << endl;
sequence2.DoContinuous(particle, track);
cout << "-->dodisc" << endl;
cout << "-->done" << endl;
const int nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(particle, track); }
for (int i = 0; i < nData; i++) {
cout << "data[" << i << "]=" << particle.p[i] << endl;
}
cout << "done" << endl;
}
SECTION("StackProcess") {
ContinuousProcess1 cp1(0);
ContinuousProcess2 cp2(3);
Process2 m2(1);
Process3 m3(2);
Stack1 s1(1);
Stack1 s2(2);
auto sequence = s1 << s2;
DummyStack stack;
const int nLoop = 20;
for (int i = 0; i < nLoop; ++i) { sequence.DoStack(stack); }
CHECK(s1.GetCount() == 20);
CHECK(s2.GetCount() == 10);
}
}
/*
Note: there is a fine-grained dedicated test-suite for SwitchProcess
in Processes/SwitchProcess/testSwtichProcess
*/
TEST_CASE("SwitchProcess") {
Process1 p1(0);
Process2 p2(0);
switch_process::SwitchProcess s(p1, p2, 10_GeV);
REQUIRE(is_switch_process_v<decltype(s)>);
}
Framework/Random/CMakeLists.txt deleted 100644 → 0
View file @ 3deb039c
set (
CORSIKArandom_SOURCES
RNGManager.cc
)
set (
CORSIKArandom_HEADERS
RNGManager.h
UniformRealDistribution.h
ExponentialDistribution.h
)
set (
CORSIKArandom_NAMESPACE
corsika/random
)
add_library (CORSIKArandom STATIC ${CORSIKArandom_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKArandom ${CORSIKArandom_NAMESPACE} ${CORSIKArandom_HEADERS})
target_link_libraries (
CORSIKArandom
INTERFACE
CORSIKAutilities
CORSIKAunits
)
target_include_directories (
CORSIKArandom
PUBLIC
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/>
)
# target dependencies on other libraries (also the header onlys)
# none
install (
TARGETS CORSIKArandom
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${CORSIKArandom_NAMESPACE}
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testRandom)
target_link_libraries (
testRandom
CORSIKArandom
CORSIKAtesting
)

data privacy - impressum