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
Showing
with 0 additions and 1789 deletions
/**
* (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.
*/
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
#include <corsika/geometry/CoordinateSystem.h>
#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
using namespace corsika::geometry;
using namespace corsika::units::si;
double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
.isApprox(EigenTransform::Identity()));
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point p1(rootCS, coordinates);
QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla, 0. * tesla};
Vector<magnetic_flux_density_d> v1(rootCS, components);
REQUIRE((p1.GetCoordinates() - coordinates).norm().magnitude() ==
Approx(0).margin(absMargin));
REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() ==
Approx(0).margin(absMargin));
/*
SECTION("unconnected CoordinateSystems") {
CoordinateSystem rootCS2 = CoordinateSystem::CreateRootCS();
REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2));
}*/
SECTION("translations") {
QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
CoordinateSystem translatedCS = rootCS.translate(translationVector);
REQUIRE(translatedCS.GetReference() == &rootCS);
REQUIRE((p1.GetCoordinates(translatedCS) + translationVector).norm().magnitude() ==
Approx(0).margin(absMargin));
// Vectors are not subject to translations
REQUIRE(
(v1.GetComponents(rootCS) - v1.GetComponents(translatedCS)).norm().magnitude() ==
Approx(0).margin(absMargin));
Point p2(translatedCS, {0_m, 0_m, 0_m});
REQUIRE(((p2 - p1).GetComponents() - translationVector).norm().magnitude() ==
Approx(0).margin(absMargin));
}
SECTION("multiple translations") {
QuantityVector<length_d> const tv1{0_m, 5_m, 0_m};
CoordinateSystem cs2 = rootCS.translate(tv1);
QuantityVector<length_d> const tv2{3_m, 0_m, 0_m};
CoordinateSystem cs3 = rootCS.translate(tv2);
QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
CoordinateSystem cs4 = cs3.translate(tv3);
REQUIRE(cs4.GetReference()->GetReference() == &rootCS);
REQUIRE(CoordinateSystem::GetTransformation(cs3, cs2).isApprox(
rootCS.translate({3_m, -5_m, 0_m}).GetTransform()));
REQUIRE(CoordinateSystem::GetTransformation(cs2, cs3).isApprox(
rootCS.translate({-3_m, +5_m, 0_m}).GetTransform()));
}
SECTION("rotations") {
QuantityVector<length_d> const axis{0_m, 0_m, 1_km};
double const angle = 90. / 180. * M_PI;
CoordinateSystem rotatedCS = rootCS.rotate(axis, angle);
REQUIRE(rotatedCS.GetReference() == &rootCS);
REQUIRE(v1.GetComponents(rotatedCS)[1].magnitude() ==
Approx((-1. * tesla).magnitude()));
// vector norm invariant under rotation
REQUIRE(v1.GetComponents(rotatedCS).norm().magnitude() ==
Approx(v1.GetComponents(rootCS).norm().magnitude()));
}
SECTION("multiple rotations") {
QuantityVector<length_d> const zAxis{0_m, 0_m, 1_km};
QuantityVector<length_d> const yAxis{0_m, 7_nm, 0_m};
QuantityVector<length_d> const xAxis{2_m, 0_nm, 0_m};
QuantityVector<magnetic_flux_density_d> components{1. * tesla, 2. * tesla,
3. * tesla};
Vector<magnetic_flux_density_d> v1(rootCS, components);
double const angle = 90. / 180. * M_PI;
CoordinateSystem rotated1 = rootCS.rotate(zAxis, angle);
CoordinateSystem rotated2 = rotated1.rotate(yAxis, angle);
CoordinateSystem rotated3 = rotated2.rotate(zAxis, -angle);
CoordinateSystem combined = rootCS.rotate(xAxis, -angle);
auto comp1 = v1.GetComponents(rotated3);
auto comp3 = v1.GetComponents(combined);
REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
}
}
TEST_CASE("Sphere") {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
SECTION("isInside") {
REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
}
}
TEST_CASE("Trajectories") {
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
Point r0(rootCS, {0_m, 0_m, 0_m});
SECTION("Line") {
Vector<SpeedType::dimension_type> v0(rootCS,
{1_m / second, 0_m / second, 0_m / second});
Line const line(r0, v0);
CHECK(
(line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(2_m, 0_m, 0_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
Trajectory<Line> base(line, 1_s);
CHECK(line.GetPosition(2_s).GetCoordinates() ==
base.GetPosition(2_s).GetCoordinates());
CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
}
SECTION("Helix") {
Vector<SpeedType::dimension_type> const vPar(
rootCS, {0_m / second, 0_m / second, 4_m / second});
Vector<SpeedType::dimension_type> const vPerp(
rootCS, {3_m / second, 0_m / second, 0_m / second});
auto const omegaC = 2 * M_PI / 1_s;
Helix const helix(r0, omegaC, vPar, vPerp);
CHECK((helix.GetPosition(1_s).GetCoordinates() -
QuantityVector<length_d>(0_m, 0_m, 4_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
CHECK((helix.GetPosition(0.25_s).GetCoordinates() -
QuantityVector<length_d>(-3_m / (2 * M_PI), -3_m / (2 * M_PI), 1_m))
.norm()
.magnitude() == Approx(0).margin(absMargin));
Trajectory<Helix> const base(helix, 1_s);
CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
base.GetPosition(1234_s).GetCoordinates());
CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(5));
}
}
/**
* (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_BufferedSink_h_
#define _include_BufferedSink_h_
namespace corsika::logging {
namespace sink {
/**
Output buffer template. NoBuffer does nothingk.
*/
/*
struct NoBuffer {
inline bool Test(const std::string&) const { return false; }
inline std::string GetString() const { return std::string(""); }
inline void Clear() {}
inline void Add(const std::string&) {}
};
*/
/**
Output buffer template. StdBuffer records fSize characters in
local memeory before passing it on to further output stages.
*/
struct StdBuffer {
StdBuffer(const int size)
: fSize(size) {}
inline bool Test(const std::string& s) {
return int(fBuffer.tellp()) + s.length() < fSize;
}
inline std::string GetString() const { return fBuffer.str(); }
inline void Clear() { fBuffer.str(""); }
inline void Add(const std::string& s) { fBuffer << s; }
private:
unsigned int fSize;
std::ostringstream fBuffer;
};
/**
Definition of Sink for log output.
*/
template <typename TStream, typename TBuffer = StdBuffer>
class BufferedSink {
public:
BufferedSink(TStream& out, TBuffer buffer = {})
: fOutput(out)
, fBuffer(std::move(buffer)) {}
void operator<<(const std::string& msg) {
if (!fBuffer.Test(msg)) {
fOutput << fBuffer.GetString();
fBuffer.Clear();
}
if (!fBuffer.Test(msg))
fOutput << msg;
else
fBuffer.Add(msg);
}
void Close() { fOutput << fBuffer.GetString(); }
private:
TStream& fOutput;
TBuffer fBuffer;
};
typedef BufferedSink<std::ostream, StdBuffer> BufferedSinkStream;
} // namespace sink
} // namespace corsika::logging
#endif
add_library (CORSIKAlogging INTERFACE)
# namespace of library -> location of header files
set (
CORSIKAlogging_NAMESPACE
corsika/logging
)
# header files of this library
set (
CORSIKAlogging_HEADERS
Logger.h
Sink.h
MessageOn.h
MessageOff.h
NoSink.h
Sink.h
BufferedSink.h
)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAlogging ${CORSIKAlogging_NAMESPACE} ${CORSIKAlogging_HEADERS})
# include directive for upstream code
target_include_directories (
CORSIKAlogging
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/>
)
# install library
install (
FILES ${CORSIKAlogging_HEADERS}
DESTINATION include/${CORSIKAlogging_NAMESPACE}
)
# ----------------
# code unit testing
add_executable (
testLogging
testLogging.cc
)
target_link_libraries (
testLogging
CORSIKAlogging
CORSIKAthirdparty # for catch2
)
add_test (
NAME testLogging
COMMAND testLogging
)
/**
* (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 Logger.h
Everything around logfile generation and text output.
*/
#ifndef _include_logger_h_
#define _include_logger_h_
#include <iosfwd>
#include <sstream>
#include <string>
#include <typeinfo>
#include <boost/format.hpp>
#include <corsika/logging/BufferedSink.h>
#include <corsika/logging/MessageOff.h>
#include <corsika/logging/MessageOn.h>
#include <corsika/logging/NoSink.h>
#include <corsika/logging/Sink.h>
using namespace std;
using namespace boost;
namespace corsika::logging {
/**
@class Logger
Defines one stream to accept messages, and to wrote those into
TSink. The helper class MessageOn will convert input at
compile-time into message strings. The helper class MessageOff,
will just do nothing and will be optimized away at compile time.
*/
template <typename MSG = MessageOn, typename TSink = sink::NoSink>
class Logger : private MSG {
using MSG::Message;
public:
// Logger() : fName("") {}
Logger(const std::string color, const std::string name, TSink& sink)
: fSink(sink)
, fName(color + "[" + name + "]\033[39m ") {}
~Logger() { fSink.Close(); }
// Logger(const Logger&) = delete;
/**
Function to add string-concatenation of all inputs to output
sink.
*/
template <typename... Strings>
void Log(const Strings&... inputs) {
fSink << MSG::Message(inputs...);
}
const std::string& GetName() const { return fName; }
private:
TSink& fSink;
std::string fName;
};
} // namespace corsika::logging
/**
* @def LOG(...)
*
* This is the main interface to the logging facilities. If Logger
* object are defined (e.g. log1) use as
* @example LOG(log1, "var1=", variable1int, "var2=", variabl2double)
* for arbitrary long sequence
* of arguments. This may also include boost::format objects the
* output is concatenated, if log1 is switched off at compile time,
* the whole LOG command is optimized away by the compiler.
*/
#define LOG(__LOGGER, ...) \
__LOGGER.Log(__LOGGER.GetName(), __FILE__, ":", __LINE__, " (", __func__, ") -> ", \
##__VA_ARGS__);
#endif
/**
* (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_MessageOff_h_
#define _include_MessageOff_h_
namespace corsika::logging {
/**
Helper class to ignore all arguments to MessagesOn::Message and
always return empty string "".
*/
class MessageOff {
protected:
template <typename First, typename... Strings>
std::string Message(const First&, const Strings&...) {
return "";
}
};
} // namespace corsika::logging
#endif
/**
* (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_MessageOn_h_
#define _include_MessageOn_h_
namespace corsika::logging {
/**
Helper class to convert all input arguments of MessageOn::Message
into string-concatenated version and return this as string.
*/
class MessageOn {
protected:
std::string Message() { return "\n"; }
template <typename First, typename... Strings>
std::string Message(const First& arg, const Strings&... rest) {
std::ostringstream ss;
ss << arg << Message(rest...);
return ss.str();
}
template <typename... Strings>
std::string Message(const int& arg, const Strings&... rest) {
return std::to_string(arg) + Message(rest...);
}
template <typename... Strings>
std::string Message(const double& arg, const Strings&... rest) {
return std::to_string(arg) + Message(rest...);
}
template <typename... Strings>
std::string Message(char const* arg, const Strings&... rest) {
return std::string(arg) + Message(rest...);
}
template <typename... Strings>
std::string Message(const std::string& arg, const Strings&... rest) {
return arg + Message(rest...);
}
// ----------------------
// boost format
template <typename... Strings>
std::string Message(const boost::format& fmt, const Strings&... rest) {
boost::format FMT(fmt);
return bformat(FMT, rest...);
}
template <typename Arg, typename... Strings>
std::string bformat(boost::format& fmt, const Arg& arg, const Strings&... rest) {
fmt % arg;
return bformat(fmt, rest...);
}
std::string bformat(boost::format& fmt) { return fmt.str() + "\n"; }
};
} // namespace corsika::logging
#endif
/**
* (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_NoSink_h_
#define _include_NoSink_h_
namespace corsika::logging {
namespace sink {
struct NoSink {
inline void operator<<(const std::string&) {}
inline void Close() {}
};
} // namespace sink
} // namespace corsika::logging
#endif
/**
* (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_Sink_h_
#define _include_Sink_h_
namespace corsika::logging {
/**
a sink for the logger must implement the two functions
operator<<(const std::string&)
and
Close()
See example: NoSink
*/
namespace sink {
/**
Definition of Sink for log output.
*/
template <typename TStream>
class Sink {
public:
Sink(TStream& out)
: fOutput(out) {}
void operator<<(const std::string& msg) { fOutput << msg; }
void Close() {}
private:
TStream& fOutput;
};
typedef Sink<std::ostream> SinkStream;
} // namespace sink
} // namespace corsika::logging
#endif
/**
* (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>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
TEST_CASE("Logging", "[Logging]") {
SECTION("sectionOne") {}
SECTION("sectionTwo") {}
SECTION("sectionThree") {}
}
add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl
COMMAND ${PROJECT_SOURCE_DIR}/Framework/Particles/pdxml_reader.py
${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleData.xml
${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleClassNames.xml
DEPENDS pdxml_reader.py
ParticleData.xml
ParticleClassNames.xml
WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Framework/Particles/
COMMENT "Read PYTHIA8 particle data and produce C++ source code GeneratedParticleProperties.inc"
VERBATIM
)
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
add_executable (
testParticles
testParticles.cc
${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
)
target_link_libraries (
testParticles
CORSIKAparticles
CORSIKAunits
CORSIKAthirdparty # for catch2
)
add_test (
NAME
testParticles
COMMAND
testParticles
)
/**
* (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>
namespace corsika::particles::io {
std::ostream& operator<<(std::ostream& stream, Code const p) {
return stream << GetName(p);
}
} // namespace corsika::particles::io
/**
* (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 <iostream>
#include <type_traits>
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
/**
* @namespace particle
*
* The properties of all elementary particles is stored here. The data
* is taken from the Pythia ParticleData.xml file.
*
*/
namespace corsika::particles {
using corsika::units::si::second;
enum class Code : int16_t;
using PDGCodeType = int16_t;
using CodeIntType = std::underlying_type<Code>::type;
// forward declarations to be used in GeneratedParticleProperties
int16_t constexpr GetElectricChargeNumber(Code const);
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const);
corsika::units::si::MassType constexpr GetMass(Code const);
PDGCodeType constexpr GetPDG(Code const);
constexpr std::string const& GetName(Code const);
corsika::units::si::TimeType constexpr GetLifetime(Code const);
#include <corsika/particles/GeneratedParticleProperties.inc>
/*!
* returns mass of particle
*/
corsika::units::si::MassType constexpr GetMass(Code const p) {
return masses[static_cast<CodeIntType const>(p)];
}
PDGCodeType constexpr GetPDG(Code const p) {
return pdg_codes[static_cast<CodeIntType const>(p)];
}
/*!
* returns electric charge of particle / (e/3).
*/
int16_t constexpr GetElectricChargeNumber(Code const p) {
return electric_charges[static_cast<CodeIntType const>(p)];
}
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
return GetElectricChargeNumber(p) * (corsika::units::si::constants::e / 3.);
}
constexpr std::string const& GetName(Code const p) {
return names[static_cast<CodeIntType const>(p)];
}
corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
return lifetime[static_cast<CodeIntType const>(p)];
}
namespace io {
std::ostream& operator<<(std::ostream& stream, Code const p);
} // namespace io
} // namespace corsika::particles
// to inject the operator<< into the root namespace
using namespace corsika::particles::io;
#endif
#!/usr/bin/env python3
import sys, math, itertools, re, csv, pprint
import xml.etree.ElementTree as ET
from collections import OrderedDict
import pickle
##############################################################
#
# reading xml input data, return line by line particle data
#
def parse(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"]
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(0)
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)
##############################################################
#
# returns dict with particle codes and class names
#
def class_names(filename):
tree = ET.parse(filename)
root = tree.getroot()
map = {}
for particle in root.iter("particle"):
name = particle.attrib["classname"]
pdg_id = int(particle.attrib["pdgID"])
map[pdg_id] = name
return map
##############################################################
#
# Automatically produce a string qualifying as C++ class name
#
# This function produces names of type "DELTA_PLUS_PLUS"
#
def c_identifier(name):
orig = name
name = name.upper()
for c in "() /":
name = name.replace(c, "_")
name = name.replace("BAR", "_BAR")
name = name.replace("0", "_0")
name = name.replace("*", "_STAR")
name = name.replace("'", "_PRIME")
name = name.replace("+", "_PLUS")
name = name.replace("-", "_MINUS")
while True:
tmp = name.replace("__", "_")
if tmp == name:
break
else:
name = tmp
pattern = re.compile(r'^[A-Z_][A-Z_0-9]*$')
if pattern.match(name):
return name.strip("_")
else:
raise Exception("could not generate C identifier for '{:s}'".format(orig))
##############################################################
#
# Automatically produce a string qualifying as C++ class name
#
# This function produces names of type "DeltaPlusPlus"
#
def c_identifier_camel(name):
orig = name
name = name[0].upper() + name[1:].lower() # all lower case
for c in "() /": # replace funny characters
name = name.replace(c, "_")
name = name.replace("bar", "Bar")
name = name.replace("*", "Star")
name = name.replace("'", "Prime")
name = name.replace("+", "Plus")
name = name.replace("-", "Minus")
# move "Bar" to end of name
ibar = name.find('Bar')
if ibar > 0 and ibar < len(name)-3:
name = name[:ibar] + name[ibar+3:] + 'Bar'
# cleanup "_"s
while True:
tmp = name.replace("__", "_")
if tmp == name:
break
else:
name = tmp
name.strip("_")
# remove all "_", if this does not by accident concatenate two number
istart = 0
while True:
i = name.find('_', istart)
if i < 1 or i > len(name)-1:
break
istart = i
if name[i-1].isdigit() and name[i+1].isdigit():
# there is a number on both sides
break
name = name[:i] + name[i+1:]
# and last, for example: make NuE out of Nue
if name[i-1].islower() and name[i].islower():
if i < len(name)-1:
name = name[:i] + name[i].upper() + name[i+1:]
else:
name = name[:i] + name[i].upper()
# check if name is valid C++ identifier
pattern = re.compile(r'^[a-zA-Z_][a-zA-Z_0-9]*$')
if pattern.match(name):
return name
else:
raise Exception("could not generate C identifier for '{:s}': result '{:s}'".format(orig, name))
##########################################################
#
# returns dict containing all data from pythia-xml input
#
def build_pythia_db(filename, classnames):
particle_db = OrderedDict()
counter = itertools.count(0)
for (pdg, name, mass, electric_charge, antiName, lifetime) in parse(filename):
c_id = "Unknown"
if pdg in classnames:
c_id = classnames[pdg]
else:
c_id = c_identifier_camel(name) # the camel case names
particle_db[c_id] = {
"name" : name,
"antiName" : antiName,
"pdg" : pdg,
"mass" : mass, # in GeV
"electric_charge" : electric_charge, # in e/3
"lifetime" : lifetime,
"ngc_code" : next(counter)
}
return particle_db
###############################################################
#
# return string with enum of all internal particle codes
#
def gen_internal_enum(pythia_db):
string = ("enum class Code : int16_t {\n"
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
for k in filter(lambda k: "ngc_code" in pythia_db[k], pythia_db):
last_ngc_id = pythia_db[k]['ngc_code']
string += " {key:s} = {code:d},\n".format(key = k, code = last_ngc_id)
string += (" LastParticle = {:d},\n" # identifier for eventual loops...
"}};").format(last_ngc_id + 1)
if last_ngc_id > 0x7fff: # does not fit into int16_t
raise Exception("Integer overflow in internal particle code definition prevented!")
return string
###############################################################
#
# return string with all data arrays
#
def gen_properties(pythia_db):
# number of particles, size of tables
string = "static constexpr std::size_t size = {size:d};\n".format(size = len(pythia_db))
string += "\n"
# particle masses table
string += "static constexpr std::array<corsika::units::si::MassType const, size> masses = {\n"
for p in pythia_db.values():
string += " {mass:f} * (1e9 * corsika::units::si::constants::eV / corsika::units::si::constants::cSquared), // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n"
# PDG code table
string += "static constexpr std::array<PDGCodeType const, size> pdg_codes = {\n"
for p in pythia_db.values():
string += " {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name'])
string += "};\n"
# name string table
string += "static const std::array<std::string const, size> names = {\n"
for p in pythia_db.values():
string += " \"{name:s}\",\n".format(name = p['name'])
string += "};\n"
# electric charges table
string += "static constexpr std::array<int16_t, size> electric_charges = {\n"
for p in pythia_db.values():
string += " {charge:d},\n".format(charge = p['electric_charge'])
string += "};\n"
# anti-particle table
# string += "static constexpr std::array<size, size> anti_particle = {\n"
# for p in pythia_db.values():
# string += " {anti:d},\n".format(charge = p['anti_particle'])
# string += "};\n"
# lifetime
string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n"
for p in pythia_db.values():
if p['lifetime'] == float("Inf") :
string += " std::numeric_limits<double>::infinity() * corsika::units::si::second, \n"
else :
string += " {tau:f} * corsika::units::si::second, \n".format(tau = p['lifetime'])
string += "};\n"
return string
###############################################################
#
# return string with a list of classes for all particles
#
def gen_classes(pythia_db):
string = "// list of C++ classes to access particle properties"
for cname in pythia_db:
antiP = 'Unknown'
for cname_anti in pythia_db:
if (pythia_db[cname_anti]['name'] == pythia_db[cname]['antiName']):
antiP = cname_anti
break
string += "\n";
string += "/** @class " + cname + "\n\n"
string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\n"
string += " * - pdg=" + str(pythia_db[cname]['pdg']) +"\n"
string += " * - mass=" + str(pythia_db[cname]['mass']) + " GeV/c² \n"
string += " * - charge= " + str(pythia_db[cname]['electric_charge']/3) + " \n"
string += " * - name=" + str(cname) + "\n"
string += " * - anti=" + str(antiP) + "\n"
string += "*/\n\n"
string += "class " + cname + " {\n"
string += " public:\n"
string += " static constexpr Code GetCode() { return Type; }\n"
string += " static constexpr corsika::units::si::MassType GetMass() { return corsika::particles::GetMass(Type); }\n"
string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetElectricCharge(Type); }\n"
string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n"
string += " static std::string const& GetName() { return corsika::particles::GetName(Type); }\n"
string += " static constexpr Code GetAntiParticle() { return AntiType; }\n"
string += " static constexpr Code Type = Code::" + cname + ";\n"
string += " static constexpr Code AntiType = Code::" + antiP + ";\n"
string += " private:\n"
string += " static constexpr CodeIntType TypeIndex = static_cast<CodeIntType const>(Type);\n"
string += "};\n"
return string
###############################################################
#
#
def inc_start():
string = ('// generated by pdxml_reader.py\n'
'// MANUAL EDITS ON OWN RISK\n')
return string
###############################################################
#
#
def inc_end():
string = ""
return string
###################################################################
#
# Serialize pythia_db into 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)
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))
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
print(gen_internal_enum(pythia_db), file=f)
print(gen_properties(pythia_db), file=f)
print(gen_classes(pythia_db), file=f)
print(inc_end(), file=f)
with open("pythia_db.pkl", "wb") as f:
serialize_pythia_db(pythia_db, f)
/**
* (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>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika::units::si;
using namespace corsika::particles;
TEST_CASE("ParticleProperties", "[Particles]") {
SECTION("Types") {
REQUIRE(Electron::GetCode() == Code::Electron);
REQUIRE(Positron::GetCode() == Code::Positron);
REQUIRE(Proton::GetCode() == Code::Proton);
REQUIRE(Neutron::GetCode() == Code::Neutron);
REQUIRE(Gamma::GetCode() == Code::Gamma);
REQUIRE(PiPlus::GetCode() == Code::PiPlus);
}
SECTION("Masses") {
REQUIRE(Electron::GetMass() / (511_keV / constants::cSquared) == Approx(1));
REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
}
SECTION("Charges") {
REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
REQUIRE(GetElectricCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
}
SECTION("Names") {
REQUIRE(Electron::GetName() == "e-");
REQUIRE(PiMinus::GetName() == "pi-");
}
SECTION("PDG") {
REQUIRE(GetPDG(Code::PiPlus) == 211);
REQUIRE(GetPDG(Code::DPlus) == 411);
REQUIRE(GetPDG(Code::NuMu) == 14);
REQUIRE(GetPDG(Code::NuE) == 12);
REQUIRE(GetPDG(Code::MuMinus) == 13);
}
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)));
}
}
/**
* (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
namespace corsika::process {
/**
\class BaseProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type BaseProcess<T>
*/
template <typename derived>
struct BaseProcess {
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
};
template <typename T>
struct is_base {
static const bool value = false;
};
template <typename T>
struct is_base<BaseProcess<T>> {
static const bool value = true;
};
} // namespace corsika::process
#endif
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
ContinuousProcess.h
DiscreteProcess.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}
)
#-- -- -- -- -- -- -- --
#code unit testing
add_executable (
testProcessSequence
testProcessSequence.cc
)
target_link_libraries (
testProcessSequence
CORSIKAsetup
CORSIKAgeometry
CORSIKAprocesssequence
CORSIKAthirdparty # for catch2
)
add_test (
NAME testProcessSequence
COMMAND testProcessSequence
)
/**
* (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
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); }
};
} // namespace corsika::process
#endif
/**
* (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_discreteprocess_h_
#define _include_corsika_discreteprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
#include <iostream> // debug
namespace corsika::process {
/**
\class DiscreteProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type DiscreteProcess<T>
*/
template <typename derived>
struct DiscreteProcess {
// DiscreteProcess() {
// static_assert(mustProvide<derived>::mustProvide, "");
//}
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 DoDiscrete...
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle&, Stack&) const; // {}
// private:
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T&, S&) const {
std::cout << "yeah" << std::endl;
return EProcessReturn::eOk;
} // find out how to make this FINAL
// void DoContinuous;
};
} // namespace corsika::process
#endif
/**
* (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 {
eOk = 1,
eParticleAbsorbed = 2,
};
} // namespace corsika::process
#endif
/**
* (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/ContinuousProcess.h>
#include <corsika/process/DiscreteProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/setup/SetupTrajectory.h>
#include <variant>
//#include <type_traits> // still needed ?
using corsika::setup::Trajectory;
namespace corsika::process {
/* namespace detail { */
/* /\* template<typename TT1, typename TT2, typename Type = void> *\/ */
/* /\* struct CallHello { *\/ */
/* /\* static void Call(const TT1&, const TT2&) { *\/ */
/* /\* std::cout << "normal" << std::endl; *\/ */
/* /\* } *\/ */
/* /\* }; *\/ */
/* /\* template<typename TT1, typename TT2> *\/ */
/* /\* struct CallHello<TT1, TT2, typename
* std::enable_if<std::is_base_of<ContinuousProcess<TT2>, TT2>::value>::type> *\/ */
/* /\* { *\/ */
/* /\* static void Call(const TT1&, const TT2&) { *\/ */
/* /\* std::cout << "special" << std::endl; *\/ */
/* /\* } *\/ */
/* /\* }; *\/ */
/* template<typename T1, typename T2, typename Particle, typename Trajectory, typename
* Stack> //, typename Type = void> */
/* struct DoContinuous { */
/* static EProcessReturn Call(const T1& A, const T2& B, Particle& p, Trajectory& t,
* Stack& s) { */
/* EProcessReturn ret = EProcessReturn::eOk; */
/* if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) { */
/* A.DoContinuous(p, t, s); */
/* } */
/* if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) { */
/* B.DoContinuous(p, t, s); */
/* } */
/* return ret; */
/* } */
/* }; */
/* /\* */
/* template<typename T1, typename T2, typename Particle, typename Trajectory, typename
* Stack> */
/* struct DoContinuous<T1,T2,Particle,Trajectory,Stack, typename
* std::enable_if<std::is_base_of<DiscreteProcess<T1>, T1>::value>::type> { */
/* static EProcessReturn Call(const T1& A, const T2& B, Particle& p, Trajectory& t,
* Stack& s) { */
/* EProcessReturn ret = EProcessReturn::eOk; */
/* A.DoContinuous(p, t, s); */
/* B.DoContinuous(p, t, s); */
/* return ret; */
/* } */
/* }; */
/* template<typename T1, typename T2, typename Particle, typename Trajectory,
* typename Stack> */
/* struct DoContinuous<T1,T2,Particle,Trajectory,Stack, typename
* std::enable_if<std::is_base_of<DiscreteProcess<T2>, T2>::value>::type> { */
/* static EProcessReturn Call(const T1& A, const T2&, Particle& p, Trajectory& t,
* Stack& s) { */
/* EProcessReturn ret = EProcessReturn::eOk; */
/* A.DoContinuous(p, t, s); */
/* B.DoContinuous(p, t, s); */
/* return ret; */
/* } */
/* }; */
/* *\/ */
/* template<typename T1, typename T2, typename Particle, typename Stack>//, typename
* Type = void> */
/* struct DoDiscrete { */
/* static EProcessReturn Call(const T1& A, const T2& B, Particle& p, Stack& s) { */
/* if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) { */
/* A.DoDiscrete(p, s); */
/* } */
/* if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) { */
/* B.DoDiscrete(p, s); */
/* } */
/* return EProcessReturn::eOk; */
/* } */
/* }; */
/* /\* */
/* template<typename T1, typename T2, typename Particle, typename Stack> */
/* struct DoDiscrete<T1,T2,Particle,Stack, typename
* std::enable_if<std::is_base_of<ContinuousProcess<T1>, T1>::value>::type> { */
/* static EProcessReturn Call(const T1&, const T2& B, Particle& p, Stack& s) { */
/* // A.DoDiscrete(p, s); */
/* B.DoDiscrete(p, s); */
/* return EProcessReturn::eOk; */
/* } */
/* }; */
/* template<typename T1, typename T2, typename Particle, typename Stack> */
/* struct DoDiscrete<T1,T2,Particle,Stack, typename
* std::enable_if<std::is_base_of<ContinuousProcess<T2>, T2>::value>::type> { */
/* static EProcessReturn Call(const T1& A, const T2&, Particle& p, Stack& s) { */
/* A.DoDiscrete(p, s); */
/* //B.DoDiscrete(p, s); */
/* return EProcessReturn::eOk; */
/* } */
/* }; */
/* *\/ */
/* } // end namespace detail */
/**
\class ProcessSequence
A compile time static list of processes. The compiler will
generate a new type based on template logic containing all the
elements.
\comment Using CRTP pattern,
https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
template <typename T1, typename T2>
class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
public:
const T1& A;
const T2& B;
ProcessSequence(const T1& in_A, const T2& in_B)
: A(in_A)
, B(in_B) {}
// example for a trait-based call:
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
template <typename Particle, typename Stack>
inline EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (!std::is_base_of<DiscreteProcess<T1>, T1>::value) {
A.DoContinuous(p, t, s);
}
if constexpr (!std::is_base_of<DiscreteProcess<T2>, T2>::value) {
B.DoContinuous(p, t, s);
}
return ret;
}
template <typename Particle>
inline void MinStepLength(Particle& p, Trajectory& step) const {
A.MinStepLength(p, step);
B.MinStepLength(p, step);
}
/*
template <typename Particle, typename Trajectory>
inline Trajectory Transport(Particle& p, double& length) const {
A.Transport(p, length); // todo: maybe check (?) if there is more than one Transport
// process implemented??
return B.Transport(
p, length); // need to do this also to decide which Trajectory to return!!!!
}
*/
template <typename Particle, typename Stack>
inline EProcessReturn DoDiscrete(Particle& p, Stack& s) const {
if constexpr (!std::is_base_of<ContinuousProcess<T1>, T1>::value) {
A.DoDiscrete(p, s);
}
if constexpr (!std::is_base_of<ContinuousProcess<T2>, T2>::value) {
B.DoDiscrete(p, s);
}
return EProcessReturn::eOk;
}
/// TODO the const_cast is not nice, think about the constness here
inline void Init() const {
const_cast<T1*>(&A)->Init();
const_cast<T2*>(&B)->Init();
}
};
/// the +operator assembles many BaseProcess, ContinuousProcess, and
/// DiscreteProcess objects into a ProcessSequence, all combinatorics
/// must be allowed, this is why we define a macro to define all
/// combinations here:
#define OPSEQ(C1, C2) \
template <typename T1, typename T2> \
inline const ProcessSequence<T1, T2> operator+(const C1<T1>& A, const C2<T2>& B) { \
return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef()); \
}
OPSEQ(BaseProcess, BaseProcess)
OPSEQ(BaseProcess, DiscreteProcess)
OPSEQ(BaseProcess, ContinuousProcess)
OPSEQ(ContinuousProcess, BaseProcess)
OPSEQ(ContinuousProcess, DiscreteProcess)
OPSEQ(ContinuousProcess, ContinuousProcess)
OPSEQ(DiscreteProcess, BaseProcess)
OPSEQ(DiscreteProcess, DiscreteProcess)
OPSEQ(DiscreteProcess, ContinuousProcess)
/*
template <typename T1>
struct depth_lhs
{
static const int num = 0;
};
// terminating condition
template <typename T1, typename T2>
struct depth_lhs< Sequence<T1,T2> >
{
// try to expand the left node (T1) which might be a Sequence type
static const int num = 1 + depth_lhs<T1>::num;
};
*/
/*
template <typename T1>
struct mat_ptrs
{
static const int num = 0;
inline static void
get_ptrs(const Process** ptrs, const T1& X)
{
ptrs[0] = reinterpret_cast<const Process*>(&X);
}
};
template <typename T1, typename T2>
struct mat_ptrs< Sequence<T1,T2> >
{
static const int num = 1 + mat_ptrs<T1>::num;
inline static void
get_ptrs(const Process** in_ptrs, const Sequence<T1,T2>& X)
{
// traverse the left node
mat_ptrs<T1>::get_ptrs(in_ptrs, X.A);
// get address of the matrix on the right node
in_ptrs[num] = reinterpret_cast<const Process*>(&X.B);
}
};
*/
/*
template<typename T1, typename T2>
const Process&
Process::operator=(const Sequence<T1,T2>& X)
{
int N = 1 + depth_lhs< Sequence<T1,T2> >::num;
const Process* ptrs[N];
mat_ptrs< Sequence<T1,T2> >::get_ptrs(ptrs, X);
int r = ptrs[0]->rows;
int c = ptrs[0]->cols;
// ... check that all matrices have the same size ...
set_size(r, c);
for(int j=0; j<r*c; ++j)
{
double sum = ptrs[0]->data[j];
for(int i=1; i<N; ++i)
{
sum += ptrs[i]->data[j];
}
data[j] = sum;
}
return *this;
}
*/
} // namespace corsika::process
#endif