diff --git a/MAINTAINERS.md b/MAINTAINERS.md
index 2da518cf5b6b9a8e0e3a278bf166fd246137554d..05f2e296d40db3554332511098a13b25751bc74b 100644
--- a/MAINTAINERS.md
+++ b/MAINTAINERS.md
@@ -6,7 +6,6 @@ discussions.
General and infrastructure:
- Ralf Ulrich <ralf.ulrich@kit.edu>, KIT
- Maximilian Reininghaus <maximilian.reininghaus@kit.edu>, KIT
-- Hans Dembinski <hdembins@mpi-hd.mpg.de>, Dortmund
- Antonio Augusto Alves Junior <antonio.junior@kit.edu>, KIT
High performance, GPU:
diff --git a/SCIENTIFIC_AUTHORS b/SCIENTIFIC_AUTHORS
deleted file mode 100644
index a00fe03acc791f21daa27da56ca89925febc7a9c..0000000000000000000000000000000000000000
--- a/SCIENTIFIC_AUTHORS
+++ /dev/null
@@ -1,28 +0,0 @@
-This is currently derived from
-https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/wikis/ICRC2019
-and will be revised whenever needed:
-
-- Luisa Arrabito (1)
-- Dominik Baack (2)
-- Johan Bregeon (1)
-- Hans Dembinski (3)
-- Ralph Engel (4)
-- Dieter Heck (4)
-- Tim Huege (4,5)
-- Lukas Nellen (6)
-- Tanguy Pierog (4)
-- Maximilian Reininghaus (4,7)
-- Felix Riehn (8)
-- Ralf Ulrich (4)
-- Darko Veberic (4)
-
-Affiliations:
-
-- University of Montpellier, France
-- TU Dortmund University, Germany
-- Max Planck Institute for Nuclear Physics, Heidelberg, Germany
-- Karlsruhe Institute of Technology (KIT), Germany
-- Vrije Universiteit Brussel, Belgium
-- National Autonomous University of Mexico (UNAM)
-- Instituto de TecnologÃas en Detección y AstropartÃculas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina
-- Laboratory of Instrumentation and Experimental Particles (LIP), Portugal
diff --git a/documentation/conf.py b/documentation/conf.py
index b0a15a6f281ce7008fea14840657803c94f7b6ae..d281602e9bc083acc9046491a7bfcb952303b945 100644
--- a/documentation/conf.py
+++ b/documentation/conf.py
@@ -25,13 +25,12 @@ def configureDoxyLayout(template_file, output_file, page_url, page_tile):
-read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True'
-build_version = os.environ.get('READTHEDOCS_VERSION', None)
+read_the_docs_build = os.environ.get('READTHEDOCS', None) != None
breathe_projects = {}
-doc_url = 'https://corsika-8.readthedocs.io/en/'+build_version
-
if read_the_docs_build:
+ build_version = os.environ.get('READTHEDOCS_VERSION', None)
+ doc_url = 'https://corsika-8.readthedocs.io/en/'+build_version
configureDoxyfile("Doxyfile.in", "Doxyfile", "../", "_build/workdir/doxygen")
configureDoxyLayout("DoxyLayout.in", "DoxyLayout.xml", doc_url , "CORSIKA 8 Webpage")
subprocess.call('mkdir -p _build/workdir/doxygen; doxygen Doxyfile', shell=True)
diff --git a/documentation/howto_create_module.rst b/documentation/howto_create_module.rst
new file mode 100644
index 0000000000000000000000000000000000000000..596bf723c545dc2b2af83418a99f0d1469c0c02b
--- /dev/null
+++ b/documentation/howto_create_module.rst
@@ -0,0 +1,135 @@
+Howto create new physics modules
+================================
+
+There are different types of physics modules, which you can add to the CORSIKA 8 physics process
+sequence. Modules can act on particles, secondaries or the entire stack. Modules can create new particles
+or they can modify or delete particles. They can also produce output.
+
+Types of different modules are explained in ::modules
+
+When creating new modules, we suggest to stick as close as possible to the default CORSIKA 8 coding guidelines
+and code structure. This makes code review and sharing with others not more complicated than needed.
+
+When your modules creates output, use the available CORSIKA 8 output machinery. This is not explained here. Also learn
+how to use units, and use the loggers from the very beginning. Furthermore, get aquinted with C++17.
+
+Let's consider the case of an "InteractionProcess" which will remove the projectile particle and create
+secondary particles on the stack instead. It also has a cross section in order to evaulate the probability
+with respect to other InteractionProcesses. Create a header file `SimpleProcess.hpp`, which is conceptually
+based to resemble (roughly) a Matthew-Heitler model:
+
+.. code-block:: c++
+
+ #pragma once
+
+ #include <corsika/framework/core/ParticleProperties.hpp>
+ #include <corsika/framework/core/PhysicalConstants.hpp>
+ #include <corsika/framework/core/PhysicalUnits.hpp>
+ #include <corsika/framework/core/Logging.hpp>
+ #include <corsika/framework/process/InteractionProcess.hpp>
+ #include <corsika/framework/random/RNGManager.hpp>
+
+ namespace corsika::simple_process {
+
+ class SimpleProcess : public corsika::InteractionProcess<SimpleProcess> {
+
+ public:
+ SimpleProcess();
+
+ template <typename TParticle>
+ GrammageType getInteractionLength(TParticle const&) const;
+
+ template <typename TSecondaryView>
+ void doInteraction(TSecondaryView& view) const;
+
+ private:
+ // the random number stream
+ corsika::default_prng_type& rng_ =
+ corsika::RNGManager<>::getInstance().getRandomStream("simple_process");
+
+ // the logger
+ std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_SimpleProcess");
+
+ };
+
+ }
+
+ #include <SimpleProcess.inl>
+
+And the corresponding `SimpleProcess.inl` as:
+
+.. code-block:: c++
+
+ Namespace corsika::simple_process
+ {
+
+ inline SimpleProcess::SimpleProcess() {}
+
+ template <typename TParticle>
+ inline GrammageType SimpleProcess::getInteractionLength(
+ TParticle const &particle) const
+ {
+ // this process only applies to hadrons above Ekin=100GeV
+ if (is_hadron(particle.getPID()) && particle.getKineticEnergy()>100_GeV) {
+ return 2_g/square(1_cm);
+ }
+ // otherwise its cross section is 0
+ return std::numeric_limits<double>::infinity() * 1_g/square(1_cm);
+ }
+
+ template <typename TSecondaryView>
+ inline void SimpleProcess::doInteraction(TSecondaryView &view) const
+ {
+ // the actual projectile particle, which will disappear after the call to doInteraction finished.
+ auto projectile = view.getProjectile();
+ int const nMult = 15;
+ auto pionEnergy = projectile.getEnergy() / nMult;
+
+ auto const pOrig = projectile.getPosition();
+ auto const tOrig = projectile.getTime();
+ auto const projectileMomentum = projectile.getMomentum();
+ auto const &cs = projectileMomentum.getCoordinateSystem();
+
+ for (int iMult=0; iMult<nMult; ++iMult) {
+
+ projectile.addSecondary(std::make_tuple(Code::PiPlus,
+ projectileMomentum.normalized() * sqrt(square(pionEnergy) + square(PiPlus::mass)),
+ pOrig,
+ tOrig));
+ }
+ CORSIKA_LOGGER_INFO(logger_, "Created {} new secondaries of energy {}.", nMult, pionEnergy);
+ }
+
+ }
+
+In this example, `TParticle` as used here in the SimpleModule is one
+particle on the CORSIKA8 particle stack. It has all methods you expect
+as `getPosition()`, `getEnergy()`, `getKineticEnergy()`, etc.
+
+The `TSecondaryView` provides special access to the CORSIKA8 particle
+stack as needed for a particle interaction with a projectile and
+secondaries. For example, `getProjectil()` will return the actual
+projectile particle, `addSecondary(...)` will produce a secondary of
+the projectil. The advantage of addSecondary wrt. addParticle is that
+there exists a direct reference to the projectile allowing to
+automatically copy weights, keep the cascade history,
+etc. Furthermore, you can define subsequent `SeoncariesProcesses`
+which can further process all the newly created secondaries. This
+could perform energy threshold cuts, or also calculate new weights,
+etc.
+
+The SimpleProcess is not necessarily useful for anything and its sole purpuse is to illustrate the mechanism to
+create your own processes.
+
+You can then include such a process in your programm (e.g. in vertical_EAS) with
+ - add: `#include <SimpleProcess.hpp>`
+ - register the new random stream: `RNGManager<>::getInstance().registerRandomStream("simple_process");`
+ - initialize it with `simple_process::SimpleProcess simple; `
+ - add it to the physiscs sequence i.e. via `auto extended_sequence = make_sequence(old_sequence, simple);`
+
+Please follow the style and guidelines for programming CORSIKA 8 code
+even for private project. Your code will be GPLv3, too, and thus
+should be made public to the community. Any discussion or eventual
+bug-fixing is much more complicated if there are deviations from the
+guideline.
+
diff --git a/documentation/index.rst b/documentation/index.rst
index 79440c503bd0b500d8235d5c1edd77c2bf459822..cbcb393b56f8cca704fd5379ba75ca93c342bd34 100644
--- a/documentation/index.rst
+++ b/documentation/index.rst
@@ -7,13 +7,16 @@ Welcome to the CORSIKA 8 air shower simulation framework.
:maxdepth: 2
readme_link
- modules
+ output
particles
+ particle_stack
media
units
environment
stack
utilities
+ modules
+ howto_create_module
api
diff --git a/documentation/output.rst b/documentation/output.rst
new file mode 100644
index 0000000000000000000000000000000000000000..88edc4b3e934a4bbaa0c0c47f5c136d544db02ea
--- /dev/null
+++ b/documentation/output.rst
@@ -0,0 +1,99 @@
+Output
+======
+
+The format of CORSIKA 8 is designed to allow simple and robust managmenet of large libraries, as well as high reading performance. There is a dedicated python library to help processing data.
+
+The basic structure of output is structured on the filesystem itself with a couple of subdirectories and files. Each run of CORSIKA~8 creates a library that can contain any number of showers.
+The format is equally suited for single huge showers as well as for a very high number of very low-energy showers. Each module included in the run can
+produce output inside this directory. The module output is separeted in individual user-named sub-directories, each containing files produced by the module. The file format is either yaml for basic configuration and summary data, or Apache parquet for any other (binary, compressed)
+data. Parquet is optimal for columnar/tabular data as it is produced by CORSIKA 8.
+
+One advantage of this format is that with normal filesystem utilties users can manage the libraries. On all systems there are tools available to
+directly read/process yaml as well as parquet files. If you, for example, don't need the particle data for space reasons, this is very simple to remove from a library. Individual
+output stream (modules) can be easily separated with no extra effort.
+
+For example, the output of the "vertical_EAS" example program looks like this:
+
+::
+
+ vertical_EAS_outputs/
+ config.yaml
+ summary.yaml
+ particles/
+ config.yaml
+ summary.yaml
+ particles.parquet
+
+The "vertical_EAS_outputs" and the "particles" are user-defined names and can be arranged/changed. But the type
+of data is well defined, e.g. in "particles" the data from an ObservationPlane object is stored. This is relevant,
+since it allows python to access this data in a controlled way.
+
+The top level "config.yaml" contains top-level library information:
+
+.. code-block:: YAML
+
+ name: vertical_EAS_outputs
+ creator: CORSIKA8
+ version: 8.0.0-prealpha
+
+and the "summary.yaml" is written in the very end (thus, the presence of the summary also indicates that a run is finished):
+
+.. code-block:: YAML
+
+ showers: 2
+ start time: 06/02/2021 23:46:18 HST
+ end time: 06/02/2021 23:46:42 HST
+ runtime: 00:00:24.260
+
+Each module has its own "config.yaml" and "summary.yaml" file, too.
+To handle thus output for analysis any tool of your preference is feasible. We recommend python. There is a python library accompanied with
+CORSIKA~8 to facilitate analysis and output handling (>>> is python prompt):
+
+.. code-block:: python
+
+ >>> import corsika
+ >>> lib = corsika.Library("vertical_EAS_outputs")
+ >>> lib.config # this gets the library configuration as a Python dictionary
+ {'name': 'vertical_EAS_outputs',
+ 'creator': 'CORSIKA8',
+ 'version': '8.0.0-prealpha'}
+ >>> lib.names # get a list of all registered processes in the library
+ ['particles']
+ >>> lib.summary # you can also load the summary information
+ {'showers': 1,
+ 'start time': '06/02/2021 23:46:18 HST',
+ 'end time': '06/02/2021 23:46:30 HST',
+ 'runtime': 11.13}
+ >>> lib.get("particles") # you can then get the process by its registered name.
+ ObservationPlane('particles')
+ >>> lib.get("particles").config # and you can also get its config as well
+ {'type': 'ObservationPlane',
+ 'plane': {'center': [0, 0, 6371000],
+ 'center.units': 'm',
+ 'normal': [0, 0, 1]},
+ 'x-axis': [1, 0, 0],
+ 'y-axis': [0, 1, 0],
+ 'delete_on_hit': True,
+ 'name': 'particles'}
+ >>> lib.get("particles").data # this returns the data as a Pandas data frame
+ shower pdg energy x y radius
+ 0 0 211 9.066702e+10 2.449931 -5.913341 7.093710
+ 1 0 22 2.403024e+11 -1.561504 -1.276160 2.024900
+ 2 0 211 1.306354e+11 -4.626045 -3.237780 6.009696
+ 3 0 211 1.773324e+11 -1.566567 4.172961 4.461556
+ 4 0 211 7.835374e+10 3.152863 -1.049201 3.330416
+ .. ... ... ... ... ... ...
+ >>> lib.get("particles").astype("arrow") # you can also request the data in a different format
+ pyarrow.Table
+ shower: int32 not null
+ pdg: int32 not null
+ energy: double not null
+ x: double not null
+ y: double not null
+ radius: double not null
+ >>>lib.get("particles").astype("pandas") # or astype("arrow"), or astype("pandas").to_numpy()
+
+You can locally install the corsika python analysis library from within your corsika
+source code directory by `pip3 install --user -e python pyarrow==0.17.0`. Note, the pyarrow
+version fix has shown to be needed on some older systems. You may not need this, or you may
+need additional packages, too.
diff --git a/documentation/particle_stack.rst b/documentation/particle_stack.rst
new file mode 100644
index 0000000000000000000000000000000000000000..33f22f996ceebbd735665fc2d8d088e7e0b7dc7b
--- /dev/null
+++ b/documentation/particle_stack.rst
@@ -0,0 +1,12 @@
+Particle Properties
+===================
+
+.. toctree::
+ particle_classes
+
+.. doxygengroup:: Particles
+ :project: CORSIKA8
+ :members:
+
+
+
diff --git a/documentation/particles.rst b/documentation/particles.rst
index 33f22f996ceebbd735665fc2d8d088e7e0b7dc7b..ace733680c3af9b061288897f24b87febe70ba57 100644
--- a/documentation/particles.rst
+++ b/documentation/particles.rst
@@ -1,12 +1,33 @@
-Particle Properties
-===================
+Particle Storage and Stack
+==========================
-.. toctree::
- particle_classes
+Particles in memory are stored in a Stack. The Stack handles the memory management and access to particle data
+properties. The stack can be extended with additional fields that then become part of the particle object.
-.. doxygengroup:: Particles
- :project: CORSIKA8
- :members:
+The standard Stack object in CORSIKA 8 is the NuclearStackExtension, further extended with internal
+fields like `geometry node` and `weight` or `history`. But the important part is the `NuclearStackExtension`.
+There are several ways to create particles:
+* novel particles on the stack:
+ * stack::addParticle(particle_data_type const&)
+ * stack::addParticle(nuclear_particle_data_type const&)
+ * stack::addParticle(particle_data_momentum_type const&)
+ * stack::addParticle(nuclear_particle_data_momentum_type const&)
+* secondary particles:
+ * stack::addSecondary(parent const&, particle_data_type const&)
+ * stack::addSecondary(parent const&, nuclear_particle_data_type const&)
+ * stack::addSecondary(parent const&, particle_data_momentum_type const&)
+ * stack::addSecondary(parent const&, nuclear_particle_data_momentum_type const&)
+ or directly:
+ * particle::addSecondary(particle_data_type const&)
+ * particle::addSecondary(nuclear_particle_data_type const&)
+ * particle::addSecondary(particle_data_momentum_type const&)
+ * particle::addSecondary(nuclear_particle_data_momentum_type const&)
+
+The content of these method parameters are:
+
+* particle_data_type: {PID [corsika::Code], kinetic energy [HEPEnergyType], direction [DirectionVector], position [Point], time[TimeType]}
+* particle_data_momentum_type: {PID [corsika::Code], momentum [MomentumVector], position [Point], time[TimeType]}
+* nuclear_particle_data_type: {PID [corsika::Code], kinetic energy [HEPEnergyType], direction [DirectionVector], position [Point], time[TimeType], A [unsigned int], Z [unsigned int]}
+* nuclear_particle_data_momentum_type: {PID [corsika::Code], momentum [MomentumVector], position [Point], time[TimeType], A [unsigned int], Z [unsigned int]}
-
diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index 7d39abccb0db3e54506b24b34525a1f8e599d42a..31e96c319741975b523bfb43482ce1de4d82b195 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -101,6 +101,10 @@ TEST_CASE("Geometry CoordinateSystems") {
// vector norm invariant under rotation
CHECK(v1.getComponents(rotatedCS).getNorm().magnitude() ==
Approx(v1.getComponents(rootCS).getNorm().magnitude()));
+
+ // this is not possible
+ QuantityVector<length_d> const axis_invalid{0_m, 0_m, 0_km};
+ CHECK_THROWS(make_rotation(rootCS, axis_invalid, angle));
}
SECTION("multiple rotations") {
@@ -201,6 +205,10 @@ TEST_CASE("Geometry CoordinateSystem-hirarchy") {
CoordinateSystemPtr root = get_root_CoordinateSystem();
Point const p1(root, {0_m, 0_m, 0_m}); // the origin of the root CS
+ CHECK(p1.getX(root) == 0_m);
+ CHECK(p1.getY(root) == 0_m);
+ CHECK(p1.getZ(root) == 0_m);
+
// root -> cs2
CoordinateSystemPtr cs2 = make_translation(root, {0_m, 0_m, 1_m});
Point const p2(cs2, {0_m, 0_m, -1_m});
@@ -280,6 +288,8 @@ TEST_CASE("Geometry Trajectories") {
.getNorm()
.magnitude() == Approx(0).margin(absMargin));
+ CHECK((line.getTimeFromArclength(10_m) / (10_m / v0.getNorm()) == Approx(1)));
+
auto const t = 1_s;
StraightTrajectory base(line, t);
CHECK(line.getPosition(t).getCoordinates() == base.getPosition(1.).getCoordinates());
@@ -348,8 +358,6 @@ TEST_CASE("Distance between points") {
CHECK(distance(p5, p6) / 1_m == Approx(1));
}
-TEST_CASE("Geometry Tree") {}
-
TEST_CASE("Path") {
// define a known CS
CoordinateSystemPtr root = get_root_CoordinateSystem();
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index bbd94cb347e35adfec4b2eada8199258912b6d43..243cbffa7df9bc47b2486b048be356381587d6db 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -186,6 +186,7 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
particle.setNode(nextVol);
particle.setPosition(traj.getPosition(1));
particle.setMomentum(traj.getDirection(1) * particle.getMomentum().getNorm());
+ SpeedType const speed_0 = particle.getVelocity().getNorm();
if (outer) {
// now we know we are in target volume, depending on "outer"
CHECK(traj.getLength(1) / 1_m == Approx(0).margin(1e-3));
@@ -211,6 +212,7 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
particle.getMomentum().getNorm(),
particle.getVelocity().getNorm(), traj2.getLength(1),
traj2.getLength(1) / particle.getVelocity().getNorm());
+ CHECK(speed_0 / traj2.getVelocity(1).getNorm() == Approx(1));
}
CHECK_FALSE(hit_2nd_behind); // this can never happen
// the next line is maybe an actual BUG: this should be investigated and eventually