diff --git a/AUTHORS b/AUTHORS
index 7036fa76902c4d4368664f0b745b5ea621a7e62f..517adbbef2aef965bee9260037d86cb218e3def4 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -1,2 +1,3 @@
Maximilian Reininghaus, maximilian.reininghaus@kit.edu
Ralf Ulrich, ralf.ulrich@kit.edu
+Felix Riehn, friehn@lip.pt
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 888b9a31e893b01a80170556654620f9799a59c7..8607bbcbb7c16d8f8e975b8550ca41492dee26dc 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/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;
}
diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index 67f43344c55d99fadfca0483e824f862c533db4d..cfd9fd8b0af8e36de289ba7dd8b4bd600677d3bf 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/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});
diff --git a/Documentation/Examples/helix_example.cc b/Documentation/Examples/helix_example.cc
index 02ccc7c2454da62ea2fec686f9980f923ba34ceb..068d25f1fc0fe6ea320eaa6d525a9679c006d2f5 100644
--- a/Documentation/Examples/helix_example.cc
+++ b/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;
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 420fc20f087d0cb8ba3ab884f4ae0db7efc8c469..ca93c81cde5c8aec8df2f2ee3f10df371828b108 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -44,6 +44,7 @@ add_executable (testEnvironment testEnvironment.cc)
target_link_libraries (
testEnvironment
+ CORSIKAsetup
CORSIKAenvironment
CORSIKAthirdparty # for catch2
)
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 1aea20cda07f642cabb82a0fb9eb18b6a73d0e4a..ef69116febd3d4859d443be8a005ecf82e36404b 100644
--- a/Environment/Environment.h
+++ b/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;
diff --git a/Framework/Geometry/RootCoordinateSystem.h b/Framework/Geometry/RootCoordinateSystem.h
index 4c3133becd9c59e4b50a9730b893a933657c5518..a9df34edfeda18b9c7a058503ac163a0f18d4c05 100644
--- a/Framework/Geometry/RootCoordinateSystem.h
+++ b/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
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 2f1f891814cc68d7e5f585c9992dce63f244cb8a..d5ef385881ec304aff345b8a6a89a608517be521 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/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") {
diff --git a/Framework/Particles/CMakeLists.txt b/Framework/Particles/CMakeLists.txt
index 8ea8654fad658692dbbc1819dbc2ab957cb64de5..3301139d597d1bc513e8b226db4708ca4c20f6ed 100644
--- a/Framework/Particles/CMakeLists.txt
+++ b/Framework/Particles/CMakeLists.txt
@@ -1,12 +1,14 @@
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/
diff --git a/Framework/Particles/NuclearData.xml b/Framework/Particles/NuclearData.xml
new file mode 100644
index 0000000000000000000000000000000000000000..d6d2a46dee4c0ac3beb53e8ec5ca5db94bf20680
--- /dev/null
+++ b/Framework/Particles/NuclearData.xml
@@ -0,0 +1,27 @@
+<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>
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index e7d48547f87779c8da84ea2c318de55ee8b163a8..91bb74eefd485b3f8afbd67fc0eaeabdbf4f7b42 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/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);
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index aa77d0c5a4e954f713c4cac0aeb938bcf00d8bff..ca7701d40dbb467781a63ffb304adfe8bdb2a7e3 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/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)
+
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 78f34d0beb262166e51839bbedf3b83521355d6a..bf307798b57d4998115a54321de965e6e53d8bb3 100644
--- a/Framework/Particles/testParticles.cc
+++ b/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);
+
+
+ }
+
}
diff --git a/Processes/NullModel/CMakeLists.txt b/Processes/NullModel/CMakeLists.txt
index ff8b026d993b5ccf9f98c2352a571c9e4dd3e018..300e3b571c581dce68f4983537a3ac128fcdf009 100644
--- a/Processes/NullModel/CMakeLists.txt
+++ b/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
diff --git a/Processes/NullModel/NullModel.cc b/Processes/NullModel/NullModel.cc
index 3629d522a3383abcb5954a84f2fa47d1a1d45354..631891e63d11cca0fc7485c4e549496037e0b19e 100644
--- a/Processes/NullModel/NullModel.cc
+++ b/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; }
diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h
index dc3bf504bca57e2cadcc264ec50e9903dcd9a2ad..e7eac1bb0196c77f0ced559ac6145455615fd00b 100644
--- a/Processes/NullModel/NullModel.h
+++ b/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
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index c098b025b678427f1e244fc8052d0dd889f11d08..908d7370224a4237bb1332e6b6149176eace39ee 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/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") {}
+
+ }
}
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index 31d204714dba4aa72301e97ab3f5402e1abdeb73..9224d9fe5f557aa8a9d2bd8e1973a52cea4c8830 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -1,11 +1,11 @@
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
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 13524a1a03d36e27bc4d170d61d10992b2732f06..fb2f67a6786fdc3613678baac6948f3417bd6c75 100644
--- a/Processes/Sibyll/Decay.h
+++ b/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>
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 87a2a054143bbbd99d0b97b0c9f8a27ed93ad58e..839fc4461b08c61fd7f08d3242b52a8360d3a40e 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/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();
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index f79fe7da354077f7b2934afac315ba926b89ae5a..3df5f8076fb356873af57dde343f1a13faf6371c 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -24,65 +24,66 @@ namespace corsika::process::sibyll {
int GetSize() const { return s_plist_.np; }
#warning check actual capacity of sibyll stack
- int GetCapacity() const { return 8000; }
-
- void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
- void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
- void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
- auto tmp = v.GetComponents();
- for (int idx = 0; idx < 3; ++idx)
- s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
- }
-
- int GetId(const int i) const { return s_plist_.llist[i]; }
-
- EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
-
- super_stupid::MomentumVector GetMomentum(const int i) const {
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
- corsika::geometry::QuantityVector<momentum_d> components{
- s_plist_.p[0][i] * 1_GeV / constants::c,
- s_plist_.p[1][i] * 1_GeV / constants::c,
- s_plist_.p[2][i] * 1_GeV / constants::c};
- super_stupid::MomentumVector v1(rootCS, components);
- return v1;
- }
-
- void Copy(const int i1, const int i2) {
- s_plist_.llist[i1] = s_plist_.llist[i2];
- s_plist_.p[3][i1] = s_plist_.p[3][i2];
- }
-
- protected:
- void IncrementSize() { s_plist_.np++; }
- void DecrementSize() {
- if (s_plist_.np > 0) { s_plist_.np--; }
- }
- };
-
- template <typename StackIteratorInterface>
- class ParticleInterface : public ParticleBase<StackIteratorInterface> {
- using ParticleBase<StackIteratorInterface>::GetStackData;
- using ParticleBase<StackIteratorInterface>::GetIndex;
-
- public:
- void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
- EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
- void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
- corsika::process::sibyll::SibyllCode GetPID() const {
- return static_cast<corsika::process::sibyll::SibyllCode>(
- GetStackData().GetId(GetIndex()));
- }
- super_stupid::MomentumVector GetMomentum() const {
- return GetStackData().GetMomentum(GetIndex());
- }
- void SetMomentum(const super_stupid::MomentumVector& v) {
- GetStackData().SetMomentum(GetIndex(), v);
- }
- };
-
- typedef Stack<SibStackData, ParticleInterface> SibStack;
-
-} // namespace corsika::process::sibyll
+ int GetCapacity() const { return 8000; }
+
+ void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
+ void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
+ void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
+ auto tmp = v.GetComponents();
+ for (int idx = 0; idx < 3; ++idx)
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
+ }
+
+ int GetId(const int i) const { return s_plist_.llist[i]; }
+
+ EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+
+ super_stupid::MomentumVector GetMomentum(const int i) const {
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ corsika::geometry::QuantityVector<momentum_d> components{
+ s_plist_.p[0][i] * 1_GeV / constants::c,
+ s_plist_.p[1][i] * 1_GeV / constants::c,
+ s_plist_.p[2][i] * 1_GeV / constants::c};
+ super_stupid::MomentumVector v1(rootCS, components);
+ return v1;
+ }
+
+ void Copy(const int i1, const int i2) {
+ s_plist_.llist[i1] = s_plist_.llist[i2];
+ s_plist_.p[3][i1] = s_plist_.p[3][i2];
+ }
+
+protected:
+ void IncrementSize() { s_plist_.np++; }
+ void DecrementSize() {
+ if (s_plist_.np > 0) { s_plist_.np--; }
+ }
+};
+
+template <typename StackIteratorInterface>
+class ParticleInterface : public ParticleBase<StackIteratorInterface> {
+ using ParticleBase<StackIteratorInterface>::GetStackData;
+ using ParticleBase<StackIteratorInterface>::GetIndex;
+
+public:
+ void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v); }
+ EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
+ corsika::process::sibyll::SibyllCode GetPID() const {
+ return static_cast<corsika::process::sibyll::SibyllCode>(
+ GetStackData().GetId(GetIndex()));
+ }
+ super_stupid::MomentumVector GetMomentum() const {
+ return GetStackData().GetMomentum(GetIndex());
+ }
+ void SetMomentum(const super_stupid::MomentumVector& v) {
+ GetStackData().SetMomentum(GetIndex(), v);
+ }
+};
+
+typedef Stack<SibStackData, ParticleInterface> SibStack;
+
+}
+>>>>>>> master
#endif
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index 388668e1de0d5b04c8b5fa920308d20ab9982d39..e1cfcbfaf84e4d04aff4ed35b324daacf1fc66c1 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -13,16 +13,16 @@ import pickle, sys, itertools
-# loads the pickled pythia_db (which is an OrderedDict)
-def load_pythiadb(filename):
+# loads the pickled particle_db (which is an OrderedDict)
+def load_particledb(filename):
with open(filename, "rb") as f:
- pythia_db = pickle.load(f)
- return pythia_db
+ particle_db = pickle.load(f)
+ return particle_db
#
-def read_sibyll_codes(filename, pythia_db):
+def read_sibyll_codes(filename, particle_db):
with open(filename) as f:
for line in f:
line = line.strip()
@@ -30,19 +30,19 @@ def read_sibyll_codes(filename, pythia_db):
continue
identifier, sib_code, canInteractFlag, xsType = line.split()
try:
- pythia_db[identifier]["sibyll_code"] = int(sib_code)
- pythia_db[identifier]["sibyll_canInteract"] = int(canInteractFlag)
- pythia_db[identifier]["sibyll_xsType"] = int(xsType)
+ particle_db[identifier]["sibyll_code"] = int(sib_code)
+ particle_db[identifier]["sibyll_canInteract"] = int(canInteractFlag)
+ particle_db[identifier]["sibyll_xsType"] = int(xsType)
except KeyError as e:
- raise Exception("Identifier '{:s}' not found in pythia_db".format(identifier))
+ raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
# generates the enum to access sibyll particles by readable names
-def generate_sibyll_enum(pythia_db):
+def generate_sibyll_enum(particle_db):
output = "enum class SibyllCode : int8_t {\n"
- for identifier, pData in pythia_db.items():
+ for identifier, pData in particle_db.items():
if pData.get('sibyll_code') != None:
output += " {:s} = {:d},\n".format(identifier, pData['sibyll_code'])
output += "};\n"
@@ -51,9 +51,9 @@ def generate_sibyll_enum(pythia_db):
# generates the look-up table to convert corsika codes to sibyll codes
-def generate_corsika2sibyll(pythia_db):
- string = "std::array<SibyllCodeIntType, {:d}> constexpr corsika2sibyll = {{\n".format(len(pythia_db))
- for identifier, pData in pythia_db.items():
+def generate_corsika2sibyll(particle_db):
+ string = "std::array<SibyllCodeIntType, {:d}> constexpr corsika2sibyll = {{\n".format(len(particle_db))
+ for identifier, pData in particle_db.items():
sibCode = pData.get("sibyll_code", 0)
string += " {:d}, // {:s}\n".format(sibCode, identifier if sibCode else identifier + " (not implemented in SIBYLL)")
string += "};\n"
@@ -62,9 +62,9 @@ def generate_corsika2sibyll(pythia_db):
# generates the look-up table to convert corsika codes to sibyll codes
-def generate_corsika2sibyll_xsType(pythia_db):
- string = "std::array<int, {:d}> constexpr corsika2sibyllXStype = {{\n".format(len(pythia_db))
- for identifier, pData in pythia_db.items():
+def generate_corsika2sibyll_xsType(particle_db):
+ string = "std::array<int, {:d}> constexpr corsika2sibyllXStype = {{\n".format(len(particle_db))
+ for identifier, pData in particle_db.items():
sibCodeXS = pData.get("sibyll_xsType", -1)
string += " {:d}, // {:s}\n".format(sibCodeXS, identifier if sibCodeXS else identifier + " (not implemented in SIBYLL)")
string += "};\n"
@@ -72,18 +72,18 @@ def generate_corsika2sibyll_xsType(pythia_db):
# generates the look-up table to convert sibyll codes to corsika codes
-def generate_sibyll2corsika(pythia_db) :
+def generate_sibyll2corsika(particle_db) :
string = ""
minID = 0
- for identifier, pData in pythia_db.items() :
+ for identifier, pData in particle_db.items() :
if 'sibyll_code' in pData:
minID = min(minID, pData['sibyll_code'])
string += "SibyllCodeIntType constexpr minSibyll = {:d};\n\n".format(minID)
pDict = {}
- for identifier, pData in pythia_db.items() :
+ for identifier, pData in particle_db.items() :
if 'sibyll_code' in pData:
sib_code = pData['sibyll_code'] - minID
pDict[sib_code] = identifier
@@ -104,13 +104,13 @@ def generate_sibyll2corsika(pythia_db) :
# generates the bitset for the flag whether Sibyll knows the particle
-def generate_known_particle(pythia_db):
- num_particles = len(pythia_db)
+def generate_known_particle(particle_db):
+ num_particles = len(particle_db)
num_bytes = num_particles // 32 + 1
string = "Bitset2::bitset2<{:d}> constexpr isKnown{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
numeric = 0
- for identifier, pData in reversed(pythia_db.items()):
+ for identifier, pData in reversed(particle_db.items()):
handledBySibyll = int("sibyll_code" in pData) & 0x1
numeric = (numeric << 1) | handledBySibyll
@@ -125,14 +125,14 @@ def generate_known_particle(pythia_db):
# generates the bitset for the flag whether Sibyll can use particle as projectile
-def generate_interacting_particle(pythia_db):
- num_particles = len(pythia_db)
+def generate_interacting_particle(particle_db):
+ num_particles = len(particle_db)
num_bytes = num_particles // 32 + 1
string = "Bitset2::bitset2<{:d}> constexpr canInteract{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
#string = "std::array<bool, {:d}> constexpr corsika2sibyll = {{\n".format(num_particles)
numeric = 0
- for identifier, pData in reversed(pythia_db.items()):
+ for identifier, pData in reversed(particle_db.items()):
can = 0
if 'sibyll_canInteract' in pData:
if pData['sibyll_canInteract'] > 0:
@@ -151,19 +151,19 @@ def generate_interacting_particle(pythia_db):
if __name__ == "__main__":
if len(sys.argv) != 3:
- print("usage: {:s} <pythia_db.pkl> <sibyll_codes.dat>".format(sys.argv[0]), file=sys.stderr)
+ print("usage: {:s} <particle_db.pkl> <sibyll_codes.dat>".format(sys.argv[0]), file=sys.stderr)
sys.exit(1)
print("code_generator.py for SIBYLL")
- pythia_db = load_pythiadb(sys.argv[1])
- read_sibyll_codes(sys.argv[2], pythia_db)
+ particle_db = load_particledb(sys.argv[1])
+ read_sibyll_codes(sys.argv[2], particle_db)
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
- print(generate_sibyll_enum(pythia_db), file=f)
- print(generate_corsika2sibyll(pythia_db), file=f)
- print(generate_known_particle(pythia_db), file=f)
- print(generate_sibyll2corsika(pythia_db), file=f)
- print(generate_interacting_particle(pythia_db), file=f)
- print(generate_corsika2sibyll_xsType(pythia_db), file=f)
+ print(generate_sibyll_enum(particle_db), file=f)
+ print(generate_corsika2sibyll(particle_db), file=f)
+ print(generate_known_particle(particle_db), file=f)
+ print(generate_sibyll2corsika(particle_db), file=f)
+ print(generate_interacting_particle(particle_db), file=f)
+ print(generate_corsika2sibyll_xsType(particle_db), file=f)
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 67c950e5585f831f8a95ae4f0e00f59df95c8919..5d91a5b4b4f2ba06f93fdc25bad74c095c66fac6 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -9,17 +9,18 @@
* the license.
*/
-#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/particles/ParticleProperties.h>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
-#include <iostream>
-using namespace std;
using namespace corsika;
+using namespace corsika::process::sibyll;
TEST_CASE("Sibyll", "[processes]") {
@@ -57,5 +58,51 @@ TEST_CASE("Sibyll", "[processes]") {
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::K0Long) == 3);
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::SigmaPlus) == 1);
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::PiMinus) == 2);
+ }
+}
+
+#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;
+
+TEST_CASE("SibyllInterface", "[processes]") {
+
+ 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("InteractionInterface") {
+
+ Interaction model;
+
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(particle, stack);
+ [[maybe_unused]] const double length = model.GetInteractionLength(particle, track);
+
+ }
+
+ SECTION("DecayInterface") {
+
+ Decay model;
+
+ model.Init();
+ /*[[maybe_unused]] const process::EProcessReturn ret =*/ model.DoDecay(particle, stack);
+ [[maybe_unused]] const double length = model.GetLifetime(particle);
+
}
+
}
+
diff --git a/Processes/StackInspector/CMakeLists.txt b/Processes/StackInspector/CMakeLists.txt
index b0b9d379e9b5b1c7536e1619e9c9e55f322b791c..f87cc1374c6161e98383a7e755a70190d62220e4 100644
--- a/Processes/StackInspector/CMakeLists.txt
+++ b/Processes/StackInspector/CMakeLists.txt
@@ -29,6 +29,7 @@ set_target_properties (
target_link_libraries (
ProcessStackInspector
CORSIKAunits
+ CORSIKAgeometry
CORSIKAsetup
)
@@ -53,6 +54,9 @@ add_executable (testStackInspector testStackInspector.cc)
target_link_libraries (
testStackInspector
+ ProcessStackInspector
+ CORSIKAsetup
+ CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
)
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 0dce10b08f9971c3265845a858fd8b66d943e831..d84634fbad2641982bc7024e799ac3655b1722d6 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -45,7 +45,7 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
EnergyType E = iterP.GetEnergy();
Etot += E;
geometry::CoordinateSystem& rootCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCS(); // for printout
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); // for printout
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 6a7d2fc62f248813895b3bfef4e673e010b20ad1..07b987b46f8a03a2683b2d160ac86e24ac251560 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -31,11 +31,8 @@ namespace corsika::process {
~StackInspector();
void Init();
-
- // template <typename Particle, typename Trajectory, typename Stack>
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
-
- // template <typename Particle>
+ //~ template <typename Particle>
corsika::units::si::LengthType MaxStepLength(Particle&,
corsika::setup::Trajectory&) const;
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index c098b025b678427f1e244fc8052d0dd889f11d08..51f796a24a92439d2fc11a42c1fac07cc2f8ec30 100644
--- a/Processes/StackInspector/testStackInspector.cc
+++ b/Processes/StackInspector/testStackInspector.cc
@@ -13,11 +13,43 @@
// cpp file
#include <catch2/catch.hpp>
+#include <corsika/process/stack_inspector/StackInspector.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+
#include <corsika/units/PhysicalUnits.h>
-TEST_CASE("NullModel", "[processes]") {
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using namespace corsika::units::si;
+using namespace corsika::process::stack_inspector;
+using namespace corsika;
+
+
+TEST_CASE("StackInspector", "[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);
- SECTION("blubb") {}
+ setup::Stack stack;
+ auto particle = stack.NewParticle();
+
+ SECTION("interface") {
+
+ StackInspector<setup::Stack> model(true);
+
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoContinuous(particle, track, stack);
+ [[maybe_unused]] const double length = model.MaxStepLength(particle, track);
+
+ }
}
+
+
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index 83ed893c150d6d041182bbe7ab9cd30046edb640..b68564bfa284cce6048c88532e13b6f796cb6b9d 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -28,6 +28,10 @@ add_executable (testTrackingLine testTrackingLine.cc)
target_link_libraries (
testTrackingLine
+ ProcessTrackingLine
+ CORSIKAsetup
+ CORSIKAutilities
+ CORSIKAenvironment
CORSIKAunits
CORSIKAenvironment
CORSIKAgeometry
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index a889d32a652753955ba2553012350e521b1d99dc..d865c2d1e7ac1aaf5a5692ca3a32587506c81375 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -71,7 +71,7 @@ namespace corsika::stack {
Point GetPosition() const { return GetStackData().GetPosition(GetIndex()); }
TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); }
-#warning this does not really work, nor make sense:
+#warning this does not really work, nor makes sense:
Vector<SpeedType::dimension_type> GetDirection() const {
auto P = GetMomentum();
return P / P.norm() * 1e10 * (units::si::meter / units::si::second);
@@ -139,7 +139,7 @@ namespace corsika::stack {
fDataE.push_back(0 * joule);
//#TODO this here makes no sense: see issue #48
geometry::CoordinateSystem& dummyCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
fMomentum.push_back(MomentumVector(
dummyCS, {0 * newton_second, 0 * newton_second, 0 * newton_second}));
fPosition.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc
index 28cb891bdc8dfead5e3b04e8113d46f7bc127108..21f882908a39e58ff83411e82771c3379dcdbc27 100644
--- a/Stack/SuperStupidStack/testSuperStupidStack.cc
+++ b/Stack/SuperStupidStack/testSuperStupidStack.cc
@@ -35,7 +35,7 @@ TEST_CASE("SuperStupidStack", "[stack]") {
p.SetPID(particles::Code::Electron);
p.SetEnergy(1.5_GeV);
geometry::CoordinateSystem& dummyCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
p.SetMomentum(MomentumVector(
dummyCS, {1 * newton_second, 1 * newton_second, 1 * newton_second}));
p.SetPosition(Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));