Merge branch 'docs_update' into 'master'

More documentation

See merge request !348

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: 

SCIENTIFIC_AUTHORS

-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

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)

documentation/howto_create_module.rst

+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.
+

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
    
 

documentation/output.rst

+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. 

documentation/particle_stack.rst

+Particle Properties
+===================
+
+.. toctree::
+   particle_classes
+
+.. doxygengroup:: Particles
+   :project: CORSIKA8
+   :members:
+
+
+      

documentation/particles.rst

-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]}
 
-      

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();

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