Framework/Cascade/CMakeLists.txt

-
-# namespace of library -> location of header files
-set (
-  CORSIKAcascade_NAMESPACE
-  corsika/cascade
-  )
-
-# header files of this library
-set (
-  CORSIKAcascade_HEADERS
-  Cascade.h
-  sibyll2.3c.h
-  SibStack.h
-  )
-
-#set (
-#  CORSIKAcascade_SOURCES
-#  TrackingStep.cc
-#  Cascade.cc
-#  )
-
-set (
-  CORSIKAsibyll_NAMESPACE
-  corsika/cascade
-  )
-
-set (
-  CORSIKAsibyll_SOURCES
-  sibyll2.3c.f
-  gasdev.f
-  )
-add_library (CORSIKAsibyll STATIC ${CORSIKAsibyll_SOURCES})
-
-#add_library (CORSIKAcascade STATIC ${CORSIKAcascade_SOURCES})
-add_library (CORSIKAcascade INTERFACE)
-
-
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
-
-#target_link_libraries (
-#  CORSIKAcascade
-#  CORSIKAparticles
-#  CORSIKAunits
-#  CORSIKAthirdparty # for catch2
-#  )
-
-# include directive for upstream code
-target_include_directories (
-  CORSIKAcascade
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-# install library
-install (
-  FILES ${CORSIKAcascade_HEADERS}
-  DESTINATION include/${CORSIKAcascade_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-add_executable (
-  testCascade
-  testCascade.cc
-  )
-
-target_link_libraries (
-  testCascade
-  #  CORSIKAutls
-  CORSIKArandom
-  CORSIKAsibyll
-  CORSIKAcascade
-  ProcessStackInspector
-  CORSIKAstackinterface
-  CORSIKAprocesses
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAprocesssequence
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (
-  NAME testCascade
-  COMMAND testCascade
-  )

Framework/Cascade/Cascade.cc

-
-/**
- * (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/cascade/Cascade.h>
-
-using namespace corsika::cascade;
-
-template <typename ProcessList, typename Stack> //, typename Trajectory>
-void Cascade::Init() {
-  fStack.Init();
-  fProcesseList.Init();
-}
-
-template <typename ProcessList, typename Stack> //, typename Trajectory>
-void Cascade::Run() {
-  if (!fStack.IsEmpty()) {
-    if (!fStack.IsEmpty()) {
-      Particle& p = fStack.GetNextParticle();
-      Step(p);
-    }
-    // do cascade equations, which can put new particles on Stack,
-    // thus, the double loop
-    // DoCascadeEquations(); //
-  }
-}
-
-template <typename ProcessList, typename Stack> //, typename Trajectory>
-void Cascade::Step(Particle& particle) {
-  double nextStep = fProcesseList.MinStepLength(particle);
-  corsika::geometry::LineTrajectory trajectory =
-      fProcesseList.Transport(particle, nextStep);
-  sequence.DoContinuous(particle, trajectory);
-  // whats going on here? Everywhere else DoDiscrete is passed a Stack reference as well
-  sequence.DoDiscrete(particle);
-}

Framework/Cascade/Cascade.h

-
-/**
- * (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_cascade_Cascade_h_
-#define _include_corsika_cascade_Cascade_h_
-
-#include <corsika/process/ProcessReturn.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <type_traits>
-
-#include <corsika/setup/SetupTrajectory.h>
-
-using namespace corsika::units::si;
-
-namespace corsika::cascade {
-
-  template <typename Tracking, typename ProcessList, typename Stack>
-  class Cascade {
-
-    typedef typename Stack::ParticleType Particle;
-
-    Cascade() = delete;
-
-  public:
-    Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
-        : fTracking(tr)
-        , fProcesseList(pl)
-        , fStack(stack) {
-      // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoDiscrete)>::value,
-      //"ProcessList has not function DoDiscrete.");
-      // static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoContinuous)>::value,
-      //	    "ProcessList has not function DoContinuous.");
-    }
-
-    void Init() {
-      fTracking.Init();
-      fProcesseList.Init();
-      fStack.Init();
-    }
-
-    void Run() {
-      while (!fStack.IsEmpty()) {
-        while (!fStack.IsEmpty()) {
-          Particle& pNext = *fStack.GetNextParticle();
-          Step(pNext);
-        }
-        // do cascade equations, which can put new particles on Stack,
-        // thus, the double loop
-        // DoCascadeEquations(); //
-      }
-    }
-
-    void Step(Particle& particle) {
-      corsika::setup::Trajectory step = fTracking.GetTrack(particle);
-      fProcesseList.MinStepLength(particle, step);
-
-      /// here the particle is actually moved along the trajectory to new position:
-      std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
-
-      corsika::process::EProcessReturn status =
-	fProcesseList.DoContinuous(particle, step, fStack);
-      if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
-        fStack.Delete(particle); // TODO: check if this is really needed
-      } else {
-        fProcesseList.DoDiscrete(particle, fStack);
-      }
-    }
-    
-  private:
-    Tracking& fTracking;
-    ProcessList& fProcesseList;
-    Stack& fStack;
-  };
-
-} // namespace corsika::cascade
-
-#endif

Framework/Cascade/SibStack.h

-#ifndef _include_sibstack_h_
-#define _include_sibstack_h_
-
-#include <vector>
-#include <string>
-
-#include <corsika/stack/Stack.h>
-#include <corsika/cascade/sibyll2.3c.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
-
-using namespace std;
-using namespace corsika::stack;
-
-class SibStackData {
-  
- public:
-  void Init();
-  
-  void Clear() { s_plist_.np = 0; }
-  
-  int GetSize() const { return s_plist_.np;  }
-  int GetCapacity() const { return 8000; }
-
-  
-  void SetId(const int i, const int v) { s_plist_.llist[i]=v; }
-  void SetEnergy(const int i, const double v) { s_plist_.p[3][i]=v;  }
-
-  int GetId(const int i) const { return s_plist_.llist[i]; }
-  double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
-  
-  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 double v) { GetStackData().SetEnergy(GetIndex(), v); }
-  double 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())); }
-  
-};
-
-
-typedef Stack<SibStackData, ParticleInterface> SibStack;
-
-#endif

Framework/Cascade/Step.cc

-
-/**
- * (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.
- */
-
-namespace cascade;
-
-void Cascade::Step(auto& sequence, Particle& particle) {
-  double nextStep = sequence.MinStepLength(particle);
-  Trajectory trajectory = sequence.Transport(particle, nextStep);
-  sequence.DoContinuous(particle, trajectory);
-  sequence.DoDiscrete(particle);
-}

Framework/Cascade/rndm_dbl.f

-C***********************************************************************
-C
-C    interface to PHOJET double precision random number generator 
-C    for SIBYLL \FR'14
-C
-C***********************************************************************
-      DOUBLE PRECISION FUNCTION S_RNDM(IDUMMY)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      DUMMY = dble(IDUMMY)
-      S_RNDM= PHO_RNDM(DUMMY)
-      END
-
-C***********************************************************************
-C
-C    initialization routine for double precision random number generator
-C    calls PHO_RNDIN \FR'14
-C
-C***********************************************************************
-      SUBROUTINE RND_INI
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      COMMON /RNDMGAS/ ISET
-      ISET = 0
-      CALL PHO_RNDIN(12,34,56,78)
-      END
-
-
-      DOUBLE PRECISION FUNCTION GASDEV(Idum)
-C***********************************************************************
-C     Gaussian deviation
-C***********************************************************************
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      IMPLICIT INTEGER(I-N)
-      COMMON /RNDMGAS/ ISET
-      SAVE
-      DATA ISET/0/      
-      gasdev=idum
-      IF (ISET.EQ.0) THEN
-1       V1=2.D0*S_RNDM(0)-1.D0
-        V2=2.D0*S_RNDM(1)-1.D0
-        R=V1**2+V2**2
-        IF(R.GE.1.D0)GO TO 1
-        FAC=SQRT(-2.D0*LOG(R)/R)
-        GSET=V1*FAC
-        GASDEV=V2*FAC
-        ISET=1
-      ELSE
-        GASDEV=GSET
-        ISET=0
-      ENDIF
-      RETURN
-      END
-C***********************************************************************
-      
-
-      DOUBLE PRECISION FUNCTION PHO_RNDM(DUMMY)
-C***********************************************************************
-C
-C    random number generator
-C
-C    initialization by call to PHO_RNDIN needed!
-C     
-C    the algorithm is taken from
-C      G.Marsaglia, A.Zaman: 'Toward a unversal random number generator'
-C      Florida State Univ. preprint FSU-SCRI-87-70
-C
-C    implementation by K. Hahn (Dec. 88), changed to include possibility
-C    of saving / reading generator registers to / from file (R.E. 10/98)
-C
-C    generator should not depend on the hardware (if a real has
-C    at least 24 significant bits in internal representation),
-C    the period is about 2**144,
-C
-C    internal registers:
-C       U(97),C,CD,CM,I,J  - seed values as initialized in PHO_RNDIN
-C
-C
-C***********************************************************************
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      SAVE
-
-      COMMON /PORAND/ U(97),C,CD,CM,I,J
-
- 100  CONTINUE
-      RNDMI = DUMMY
-      RNDMI = U(I)-U(J)
-      IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
-      U(I) = RNDMI
-      I    = I-1
-      IF ( I.EQ.0 ) I = 97
-      J    = J-1
-      IF ( J.EQ.0 ) J = 97
-      C    = C-CD
-      IF ( C.LT.0.D0 ) C = C+CM
-      RNDMI = RNDMI-C
-      IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
-
-      IF((ABS(RNDMI).LT.0.D0).OR.(ABS(RNDMI-1.D0).LT.1.D-10)) GOTO 100
-      PHO_RNDM = RNDMI
-
-      END
-
-
-CDECK  ID>, PHO_RNDIN
-      SUBROUTINE PHO_RNDIN(NA1,NA2,NA3,NB1)
-C***********************************************************************
-C
-C     initialization of PHO_RNDM, has to be called before using PHO_RNDM
-C
-C     input:
-C       NA1,NA2,NA3,NB1  - values for initializing the generator
-C                          NA? must be in 1..178 and not all 1;
-C                          12,34,56  are the standard values
-C                          NB1 must be in 1..168;
-C                          78  is the standard value
-C
-C***********************************************************************
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      SAVE
-
-      COMMON /PORAND/ U(97),C,CD,CM,I,J
-      MA1 = NA1
-      MA2 = NA2
-      MA3 = NA3
-      MB1 = NB1
-      I   = 97
-      J   = 33
-      DO 20 II2 = 1,97
-        S = 0.D0
-        T = 0.5D0
-        DO 10 II1 = 1,24
-          MAT  = MOD(MOD(MA1*MA2,179)*MA3,179)
-          MA1  = MA2
-          MA2  = MA3
-          MA3  = MAT
-          MB1  = MOD(53*MB1+1,169)
-          IF ( MOD(MB1*MAT,64).GE.32 ) S = S+T
-          T    = 0.5D0*T
- 10     CONTINUE
-        U(II2) = S
- 20   CONTINUE
-      C  =   362436.D0/16777216.D0
-      CD =  7654321.D0/16777216.D0
-      CM = 16777213.D0/16777216.D0
-
-      END
-
-
-CDECK  ID>, PHO_RNDSI
-      SUBROUTINE PHO_RNDSI(UIN,CIN,CDIN,CMIN,IIN,JIN)
-C***********************************************************************
-C
-C     updates internal random number generator registers using
-C     registers given as arguments
-C
-C***********************************************************************
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      SAVE
-
-      DIMENSION UIN(97)
-      COMMON /PORAND/ U(97),C,CD,CM,I,J
-      DO 10 KKK = 1,97
-        U(KKK) = UIN(KKK)
- 10   CONTINUE
-      C  = CIN
-      CD = CDIN
-      CM = CMIN
-      I  = IIN
-      J  = JIN
-
-      END
-
-
-CDECK  ID>, PHO_RNDSO
-      SUBROUTINE PHO_RNDSO(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
-C***********************************************************************
-C
-C     copies internal registers from randon number generator
-C     to arguments
-C
-C***********************************************************************
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      SAVE
-
-      DIMENSION UOUT(97)
-      COMMON /PORAND/ U(97),C,CD,CM,I,J
-      DO 10 KKK = 1,97
-        UOUT(KKK) = U(KKK)
- 10   CONTINUE
-      COUT  = C
-      CDOUT = CD
-      CMOUT = CM
-      IOUT  = I
-      JOUT  = J
-
-      END
-
-
-CDECK  ID>, PHO_RNDTE
-      SUBROUTINE PHO_RNDTE(IO)
-C***********************************************************************
-C
-C     test of random number generator PHO_RNDM
-C
-C     input:
-C       IO defines output
-C           0  output only if an error is detected
-C           1  output independend on an error
-C
-C     uses PHO_RNDSI and PHO_RNDSO to bring the random number generator
-C     to same status as it had before the test run
-C
-C***********************************************************************
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-      SAVE
-
-
-C  input/output channels
-      INTEGER LI,LO
-      COMMON /POINOU/ LI,LO
-
-      DIMENSION UU(97)
-      DIMENSION U(6),X(6),D(6)
-      DATA U / 6533892.D0 , 14220222.D0 ,  7275067.D0 ,
-     &         6172232.D0 ,  8354498.D0 , 10633180.D0 /
-
-      CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
-
-      CALL PHO_RNDIN(12,34,56,78)
-      DO 10 II1 = 1,20000
-        XX      = PHO_RNDM(SD)
- 10   CONTINUE
-
-      SD        = 0.D0
-      DO 20 II2 = 1,6
-        X(II2)  = 4096.D0*(4096.D0*PHO_RNDM(XX))
-        D(II2)  = X(II2)-U(II2)
-        SD      = SD+ABS(D(II2))
- 20   CONTINUE
-
-      CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
-
-      IF ((IO.EQ.1).OR.(ABS(SD).GT.0.D-10)) THEN
-        WRITE(LO,50) (U(I),X(I),D(I),I=1,6)
-      ENDIF
-
- 50   FORMAT(/,' PHO_RNDTE: test of the random number generator:',/,
-     &  '    expected value    calculated value     difference',/,
-     &  6(F17.1,F20.1,F15.3,/),
-     &  ' generator has the same status as before calling PHO_RNDTE',/)
-
-      END
-
-
-CDECK  ID>, PHO_RNDST
-      SUBROUTINE PHO_RNDST(MODE,FILENA)
-C***********************************************************************
-C
-C     read / write random number generator status from / to file
-C
-C     input:    MODE        1   read registers from file
-C                           2   dump registers to file
-C
-C               FILENA      file name
-C
-C***********************************************************************
-
-      IMPLICIT NONE
-
-
-
-      SAVE
-
-      INTEGER       MODE
-      CHARACTER*(*) FILENA
-
-
-C  input/output channels
-      INTEGER LI,LO
-      COMMON /POINOU/ LI,LO
-
-
-      DOUBLE PRECISION UU,CC,CCD,CCM
-      DIMENSION UU(97)
-
-      INTEGER I,II,JJ
-
-      CHARACTER*80 CH_DUMMY
-
-      IF(MODE.EQ.1) THEN
-
-        WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
-     &    'reading random number registers from file ',FILENA
-
-        OPEN(12,FILE=FILENA,ERR=1010,STATUS='OLD')
-        READ(12,*,ERR=1010) CH_DUMMY
-        DO I=1,97
-          READ(12,*,ERR=1010) UU(I)
-        ENDDO
-        READ(12,*,ERR=1010) CC
-        READ(12,*,ERR=1010) CCD
-        READ(12,*,ERR=1010) CCM
-        READ(12,*,ERR=1010) II,JJ
-        CLOSE(12)
-        CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
-
-      ELSE IF(MODE.EQ.2) THEN
-
-        WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
-     &    'dumping random number registers to file ',FILENA
-
-        OPEN(12,FILE=FILENA,ERR=1010,STATUS='UNKNOWN')
-        CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
-        WRITE(12,'(1X,A)',ERR=1020) 'random number status registers:'
-        DO I=1,97
-          WRITE(12,'(1X,1P,E28.20)',ERR=1020) UU(I)
-        ENDDO
-        WRITE(12,'(1X,1P,E28.20)',ERR=1020) CC
-        WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCD
-        WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCM
-        WRITE(12,'(1X,2I4)',ERR=1020) II,JJ
-        CLOSE(12)
-
-      ELSE
-
-        WRITE(LO,'(/,1X,2A,I6,/)') 'PHO_RNDST: ',
-     &    'called with invalid mode, nothing done (mode)',MODE
-
-      ENDIF
-
-      RETURN
-
- 1010 CONTINUE
-      WRITE(LO,'(1X,2A,A,/)') 'PHO_RNDST: ',
-     &  'cannot open or read file ',FILENA
-      RETURN
-
- 1020 CONTINUE
-      WRITE(LO,'(1X,A,A,/)') 'PHO_RNDST: ',
-     &  'cannot open or write file ',FILENA
-      RETURN
-    
-      END
-
-C----------------------------------------
-C standard generator
-C----------------------------------------
-      REAL FUNCTION S_RNDM_std(IDUMMY)
-C...Generator  from the LUND montecarlo
-C...Purpose: to generate random numbers uniformly distributed between
-C...0 and 1, excluding the endpoints.
-      COMMON/LUDATR/MRLU(6),RRLU(100)
-      SAVE /LUDATR/
-      EQUIVALENCE (MRLU1,MRLU(1)),(MRLU2,MRLU(2)),(MRLU3,MRLU(3)),
-     &(MRLU4,MRLU(4)),(MRLU5,MRLU(5)),(MRLU6,MRLU(6)),
-     &(RRLU98,RRLU(98)),(RRLU99,RRLU(99)),(RRLU00,RRLU(100))
- 
-C...  Initialize generation from given seed.
-      S_RNDM_std = real(idummy)
-      IF(MRLU2.EQ.0) THEN
-        IF (MRLU1 .EQ. 0)  MRLU1 = 19780503    ! initial seed
-        IJ=MOD(MRLU1/30082,31329)
-        KL=MOD(MRLU1,30082)
-        I=MOD(IJ/177,177)+2
-        J=MOD(IJ,177)+2
-        K=MOD(KL/169,178)+1
-        L=MOD(KL,169)
-        DO 110 II=1,97
-        S=0.
-        T=0.5
-        DO 100 JJ=1,24
-        M=MOD(MOD(I*J,179)*K,179)
-        I=J
-        J=K
-        K=M
-        L=MOD(53*L+1,169)
-        IF(MOD(L*M,64).GE.32) S=S+T
-        T=0.5*T
-  100   CONTINUE
-        RRLU(II)=S
-  110   CONTINUE
-        TWOM24=1.
-        DO 120 I24=1,24
-        TWOM24=0.5*TWOM24
-  120   CONTINUE
-        RRLU98=362436.*TWOM24
-        RRLU99=7654321.*TWOM24
-        RRLU00=16777213.*TWOM24
-        MRLU2=1
-        MRLU3=0
-        MRLU4=97
-        MRLU5=33
-      ENDIF
- 
-C...Generate next random number.
-  130 RUNI=RRLU(MRLU4)-RRLU(MRLU5)
-      IF(RUNI.LT.0.) RUNI=RUNI+1.
-      RRLU(MRLU4)=RUNI
-      MRLU4=MRLU4-1
-      IF(MRLU4.EQ.0) MRLU4=97
-      MRLU5=MRLU5-1
-      IF(MRLU5.EQ.0) MRLU5=97
-      RRLU98=RRLU98-RRLU99
-      IF(RRLU98.LT.0.) RRLU98=RRLU98+RRLU00
-      RUNI=RUNI-RRLU98
-      IF(RUNI.LT.0.) RUNI=RUNI+1.
-      IF(RUNI.LE.0.OR.RUNI.GE.1.) GOTO 130
- 
-C...Update counters. Random number to output.
-      MRLU3=MRLU3+1
-      IF(MRLU3.EQ.1000000000) THEN
-        MRLU2=MRLU2+1
-        MRLU3=0
-      ENDIF
-      S_RNDM_std=RUNI
-
-      END

Framework/Cascade/sibyll2.3c.h

-#ifndef _include_sib23c_interface_h_
-#define _include_sib23c_interface_h_
-//----------------------------------------------
-//  C++ interface for the SIBYLL event generator
-//----------------------------------------------
-//wrapper
-
-
-
-
-extern"C"{
-
-  typedef char s_name[6];
-
-  // SIBYLL particle stack (FORTRAN COMMON)
-  // variables are: np : numer of particles on stack
-  //                 p : 4momentum + mass of particles on stack
-  //             llist : id of particles on stack
-  extern struct { 
-    double p[5][8000]; 
-    int llist[8000]; 
-    int np; 
-  }s_plist_;
-
-
-  extern struct {
-    double  cbr[223+16+12+8];
-    int    kdec[1338+6*(16+12+8)];
-    int    lbarp[99];
-    int     idb[99];
-  }s_csydec_;
-
-  
-  // additional particle stack for the mother particles of unstable particles
-  // stable particles have entry zero
-  extern struct { 
-    int llist1[8000]; 
-  }s_plist1_;
-
-
-  // tables with particle properties
-  // charge, strangeness and baryon number
-  extern struct {
-    int ichp[99]; 
-    int istr[99]; 
-    int ibar[99]; 
-  }s_chp_;
-
-  // tables with particle properties
-  // mass and mass squared
-  extern struct {
-    double am[99]; 
-    double am2[99]; 
-  }s_mass1_;
-
-  // table with particle names
-  extern struct {char namp[6][99];}s_cnam_;
-
-  // debug info
-  extern struct {int ncall; int ndebug; int lun; }s_debug_;
-
-  // lund random generator setup
-  //extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
-
-
-  // sibyll main subroutine
-  void sibyll_(int&,int&,double&);
-    
-  // subroutine to initiate sibyll
-  void sibyll_ini_();
-
-  // subroutine to SET DECAYS
-  void dec_ini_();
-
-  // subroutine to initiate random number generator
-  //void rnd_ini_();
-  
-  // print event
-  void sib_list_(int&);
-
-  // decay routine
-  void decsib_();
-
-  // interaction length
-  //double fpni_(double&, int&);
-
-  void sib_sigma_hnuc_(int&,int&,double&,double&,double&);
-  void sib_sigma_hp_(int&,double&,double&,double&,double&,double*,double&,double&);
-  
-  double s_rndm_(int&);
-
-  // phojet random generator setup
-  void pho_rndin_(int&,int&,int&,int&);
-}
-#endif
-

Framework/Cascade/testCascade.cc

-
-/**
- * (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/cascade/Cascade.h>
-
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
-
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-using corsika::setup::Trajectory;
-
-#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;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::geometry;
-
-#include <iostream>
-using namespace std;
-
-static int fCount = 0;
-
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
-public:
-  ProcessSplit() {}
-
-  template <typename Particle>
-  void MinStepLength(Particle&, setup::Trajectory&) const {}
-
-  template <typename Particle, typename Stack>
-  EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
-    return EProcessReturn::eOk;
-  }
-
-  template <typename Particle, typename Stack>
-  void DoDiscrete(Particle& p, Stack& s) const {
-    EnergyType E = p.GetEnergy();
-    if (E < 85_MeV) {
-      p.Delete();
-      fCount++;
-    } else {
-      p.SetEnergy(E / 2);
-      auto pnew = s.NewParticle();
-      // s.Copy(p, pnew);
-      pnew.SetEnergy(E / 2);
-      pnew.SetPosition(p.GetPosition());
-      pnew.SetMomentum(p.GetMomentum());
-    }
-  }
-
-  void Init() { fCount = 0; }
-
-  int GetCount() const { return fCount; }
-
-private:
-};
-
-template <typename Stack>
-class NewtonTracking { // naja.. not yet
-  typedef typename Stack::ParticleType Particle;
-
-public:
-  void Init() {}
-  corsika::setup::Trajectory GetTrack(Particle& p) {
-    corsika::geometry::Vector<SpeedType::dimension_type> v =
-        p.GetDirection(); 
-    corsika::geometry::Line traj(p.GetPosition(), v);
-    {
-      CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-      cout << v.GetComponents(rootCS) << endl;
-    }
-    return corsika::geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
-  }
-};
-
-TEST_CASE("Cascade", "[Cascade]") {
-
-  NewtonTracking<setup::Stack> tracking;
-
-  stack_inspector::StackInspector<setup::Stack> p0(true);
-  ProcessSplit p1;
-  const auto sequence = p0 + p1;
-  setup::Stack stack;
-
-  corsika::cascade::Cascade EAS(tracking, sequence, stack);
-
-  CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
-  stack.Clear();
-  auto particle = stack.NewParticle();
-  EnergyType E0 = 100_GeV;
-  particle.SetEnergy(E0);
-  particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
-  particle.SetMomentum(
-      corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_MeV}));
-  EAS.Init();
-  EAS.Run();
-
-  SECTION("sectionTwo") {
-    for (int i = 0; i < 0; ++i) {
-      stack.Clear();
-      auto particle = stack.NewParticle();
-      EnergyType E0 = 100_GeV * pow(10, i);
-      particle.SetEnergy(E0);
-      EAS.Init();
-      EAS.Run();
-
-      // cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
-    }
-  }
-}

Framework/Geometry/BaseTrajectory.h

-
-/**
- * (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_BASETRAJECTORY_H
-#define _include_BASETRAJECTORY_H
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <string>
-
-namespace corsika::geometry {
-
-  /*!
-   * Interface / base class for trajectories.
-   */
-  class BaseTrajectory {
-
-    BaseTrajectory() = delete;
-
-  public:
-    BaseTrajectory(corsika::units::si::TimeType start, corsika::units::si::TimeType end)
-        : fTStart(start)
-        , fTEnd(end) {}
-
-    //!< for \f$ t = 0 \f$, the starting Point shall be returned.
-    virtual Point GetPosition(corsika::units::si::TimeType) const = 0;
-
-    //!< the Point is return from u=0 (start) to u=1 (end)
-    virtual Point GetPosition(double u) const = 0;
-
-    /*!
-     * returns the length between two points of the trajectory
-     * parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
-     */
-
-    virtual LengthType GetDistance(corsika::units::si::TimeType t1,
-                                   corsika::units::si::TimeType t2) const = 0;
-
-    virtual corsika::units::si::TimeType GetDuration(
-        corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
-      return t2 - t1;
-    }
-
-    virtual Point GetEndpoint() const { return GetPosition(fTEnd); }
-    virtual Point GetStartpoint() const { return GetPosition(fTStart); }
-
-  protected:
-    corsika::units::si::TimeType const fTStart, fTEnd;
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/BaseVector.h

-
-/**
- * (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_BASEVECTOR_H_
-#define _include_BASEVECTOR_H_
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/geometry/QuantityVector.h>
-
-namespace corsika::geometry {
-
-  /*!
-   * Common base class for Vector and Point. Currently it does basically nothing.
-   */
-
-  template <typename dim>
-  class BaseVector {
-  protected:
-    QuantityVector<dim> qVector;
-    CoordinateSystem const* cs;
-
-  public:
-    BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
-        : qVector(pQVector)
-        , cs(&pCS) {}
-
-
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/CMakeLists.txt

-
-set (
-  GEOMETRY_SOURCES
-  CoordinateSystem.cc
-  )
-
-set (
-  GEOMETRY_HEADERS
-  Vector.h
-  Point.h
-  Line.h
-  Sphere.h
-  CoordinateSystem.h
-  RootCoordinateSystem.h
-  Helix.h
-  BaseVector.h
-  QuantityVector.h
-  Trajectory.h
-  # BaseTrajectory.h
-  )
-
-set (
-  GEOMETRY_NAMESPACE
-  corsika/geometry
-  )
-
-add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAgeometry ${GEOMETRY_NAMESPACE} ${GEOMETRY_HEADERS})
-
-set_target_properties (
-  CORSIKAgeometry
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${GEOMETRY_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAgeometry
-  CORSIKAunits
-  )
-
-target_include_directories (
-  CORSIKAgeometry
-  PUBLIC    ${EIGEN3_INCLUDE_DIR}
-  INTERFACE ${EIGEN3_INCLUDE_DIR}
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS CORSIKAgeometry
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${GEOMETRY_NAMESPACE}
-  )
-
-
-# --------------------
-# code unit testing
-add_executable (testGeometry testGeometry.cc)
-
-target_link_libraries (
-  testGeometry
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (NAME testGeometry COMMAND testGeometry)
-

Framework/Geometry/CoordinateSystem.cc

-
-/**
- * (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/geometry/CoordinateSystem.h>
-
-using namespace corsika::geometry;
-
-/**
- * returns the transformation matrix necessary to transform primitives with coordinates
- * in \a pFrom to \a pTo, e.g.
- * \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
- * (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
- * the indicated CoordinateSystem).
- */
-EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom,
-                                                   CoordinateSystem const& pTo) {
-  CoordinateSystem const* a{&pFrom};
-  CoordinateSystem const* b{&pTo};
-  CoordinateSystem const* commonBase{nullptr};
-
-  while (a != b && b != nullptr) {
-    a = &pFrom;
-
-    while (a != b && a != nullptr) { a = a->GetReference(); }
-
-    if (a == b) break;
-
-    b = b->GetReference();
-  }
-
-  if (a == b && a != nullptr) {
-    commonBase = a;
-
-  } else {
-    throw std::string("no connection between coordinate systems found!");
-  }
-
-  EigenTransform t = EigenTransform::Identity();
-  auto* p = &pFrom;
-
-  while (p != commonBase) {
-    t = p->GetTransform() * t;
-    p = p->GetReference();
-  }
-
-  p = &pTo;
-
-  while (p != commonBase) {
-    t = t * p->GetTransform().inverse(Eigen::TransformTraits::Isometry);
-    p = p->GetReference();
-  }
-
-  return t;
-}

Framework/Geometry/CoordinateSystem.h

-
-/**
- * (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_COORDINATESYSTEM_H_
-#define _include_COORDINATESYSTEM_H_
-
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <Eigen/Dense>
-
-typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
-typedef Eigen::Translation<double, 3> EigenTranslation;
-
-namespace corsika::geometry {
-
-  class RootCoordinateSystem;
-  
-  using corsika::units::si::length_d;
-
-  class CoordinateSystem {
-    CoordinateSystem const* reference = nullptr;
-    EigenTransform transf;
-
-    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
-        : reference(&reference)
-        , transf(transf) {}
-        
-    CoordinateSystem()
-        : // for creating the root CS
-        transf(EigenTransform::Identity()) {}
-
-  protected:
-	static auto CreateCS() { return CoordinateSystem(); }
-	friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can creat ONE unique root CS
-
-  public:
-  
-    static EigenTransform GetTransformation(CoordinateSystem const& c1,
-                                            CoordinateSystem const& c2);
-
-    auto& operator=(const CoordinateSystem& pCS) {
-      reference = pCS.reference;
-      transf = pCS.transf;
-      return *this;
-    }
-
-    auto translate(QuantityVector<length_d> vector) const {
-      EigenTransform const translation{EigenTranslation(vector.eVector)};
-
-      return CoordinateSystem(*this, translation);
-    }
-
-    auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
-      if (axis.eVector.isZero()) {
-          throw std::string("null-vector given as axis parameter");
-      }
-      
-      EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
-
-      return CoordinateSystem(*this, rotation);
-    }
-
-    auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
-                            QuantityVector<phys::units::length_d> axis, double angle) {
-      if (axis.eVector.isZero()) {
-          throw std::string("null-vector given as axis parameter");
-      }
-                                
-      EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
-                                  EigenTranslation(translation.eVector)};
-
-      return CoordinateSystem(*this, transf);
-    }
-
-    auto const* GetReference() const { return reference; }
-
-    auto const& GetTransform() const { return transf; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Helix.h

-
-/**
- * (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_HELIX_H_
-#define _include_HELIX_H_
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <cmath>
-
-namespace corsika::geometry {
-  /*!
-   * A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
-   * Point r0 and
-   * the velocity vectors \f$ \vec{v}_{\parallel} \f$ and \f$ \vec{v}_{\perp} \f$
-   * denoting the projections of the initial velocity \f$ \vec{v}_0 \f$ parallel
-   * and perpendicular to the axis \f$ \vec{B} \f$, respectively, i.e.
-   * \f{align*}{
-        \vec{v}_{\parallel} &= \frac{\vec{v}_0 \cdot \vec{B}}{\vec{B}^2} \vec{B} \\
-        \vec{v}_{\perp} &= \vec{v}_0 - \vec{v}_{\parallel}
-     \f}
-   */
-
-  class Helix {
-
-    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
-
-    Point const r0;
-    corsika::units::si::FrequencyType const omegaC;
-    VelocityVec const vPar;
-    VelocityVec const vPerp, uPerp;
-
-    corsika::units::si::LengthType const radius;
-
-  public:
-    Helix(Point const& pR0, corsika::units::si::FrequencyType pOmegaC,
-          VelocityVec const& pvPar, VelocityVec const& pvPerp)
-        : r0(pR0)
-        , omegaC(pOmegaC)
-        , vPar(pvPar)
-        , vPerp(pvPerp)
-        , uPerp(vPerp.cross(vPar.normalized()))
-        , radius(pvPar.norm() / abs(pOmegaC)) {}
-
-    Point GetPosition(corsika::units::si::TimeType t) const {
-      return r0 + vPar * t +
-             (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
-    }
-
-    auto GetRadius() const { return radius; }
-
-    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
-                                  corsika::units::si::TimeType t2) const {
-      return (vPar + vPerp).norm() * (t2 - t1);
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Line.h

-
-/**
- * (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_LINETRAJECTORY_H
-#define _include_LINETRAJECTORY_H
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  class Line {
-    
-    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
-
-    Point const r0;
-    VelocityVec const v0;
-
-  public:
-    Line(Point const& pR0, VelocityVec const& pV0)
-        : r0(pR0)
-        , v0(pV0) {}
-
-    Point GetPosition(corsika::units::si::TimeType t) const {
-      return r0 + v0 * t;
-    }
-    
-    LengthType GetDistanceBetween(corsika::units::si::TimeType t1,
-				  corsika::units::si::TimeType t2) const {
-      // assert(t2 >= t1);
-      return v0.norm() * (t2 - t1);
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Point.h

-
-/**
- * (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_POINT_H_
-#define _include_POINT_H_
-
-#include <corsika/geometry/BaseVector.h>
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  using corsika::units::si::length_d;
-  using corsika::units::si::LengthType;
-
-  /*!
-   * A Point represents a point in position space. It is defined by its
-   * coordinates with respect to some CoordinateSystem.
-   */
-  class Point : public BaseVector<phys::units::length_d> {
-  public:
-    Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector)
-        : BaseVector<phys::units::length_d>(pCS, pQVector) {}
-
-    Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
-        : BaseVector<phys::units::length_d>(cs, {x, y, z}) {}
-
-    // TODO: this should be private or protected, we don NOT want to expose numbers without reference to outside:
-    auto GetCoordinates() const { return BaseVector<phys::units::length_d>::qVector; }
-
-    /// this always returns a QuantityVector as triple
-    auto GetCoordinates(CoordinateSystem const& pCS) const {
-      if (&pCS == BaseVector<phys::units::length_d>::cs) {
-        return BaseVector<phys::units::length_d>::qVector;
-      } else {
-        return QuantityVector<phys::units::length_d>(
-            CoordinateSystem::GetTransformation(*BaseVector<phys::units::length_d>::cs,
-                                                pCS) *
-            BaseVector<phys::units::length_d>::qVector.eVector);
-      }
-    }
-
-    /*!
-     * transforms the Point into another CoordinateSystem by changing its
-     * coordinates interally
-     */
-    void rebase(CoordinateSystem const& pCS) {
-      BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS);
-      BaseVector<phys::units::length_d>::cs = &pCS;
-    }
-
-    Point operator+(Vector<phys::units::length_d> const& pVec) const {
-      return Point(
-          *BaseVector<phys::units::length_d>::cs,
-          GetCoordinates() + pVec.GetComponents(*BaseVector<phys::units::length_d>::cs));
-    }
-
-    /*!
-     * returns the distance Vector between two points
-     */
-    Vector<phys::units::length_d> operator-(Point const& pB) const {
-      auto& cs = *BaseVector<phys::units::length_d>::cs;
-      return Vector<phys::units::length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/QuantityVector.h

-
-/**
- * (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_QUANTITYVECTOR_H_
-#define _include_QUANTITYVECTOR_H_
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <Eigen/Dense>
-
-#include <iostream>
-#include <utility>
-
-namespace corsika::geometry {
-
-  /*!
-   * A QuantityVector is a three-component container based on Eigen::Vector3d
-   * with a phys::units::si::dimension. Arithmethic operators are defined that
-   * propagate the dimensions by dimensional analysis.
-   */
-
-  template <typename dim>
-  class QuantityVector {
-  protected:
-    // todo: check if we need to move "quantity" into namespace corsika::units
-    using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity
-                                                         // corresponding to the dimension
-
-  public:
-    Eigen::Vector3d eVector; //!< the actual container where the raw numbers are stored
-
-    typedef dim dimension; //!< should be a phys::units::dimension
-
-    QuantityVector(Quantity a, Quantity b, Quantity c)
-        : eVector{a.magnitude(), b.magnitude(), c.magnitude()} {}
-
-    QuantityVector(Eigen::Vector3d pBareVector)
-        : eVector(pBareVector) {}
-
-    auto operator[](size_t index) const {
-      return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
-    }
-
-    auto norm() const {
-      return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
-    }
-
-    auto squaredNorm() const {
-      using QuantitySquared =
-          decltype(std::declval<Quantity>() * std::declval<Quantity>());
-      return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
-    }
-
-    auto operator+(QuantityVector<dim> const& pQVec) const {
-      return QuantityVector<dim>(eVector + pQVec.eVector);
-    }
-
-    auto operator-(QuantityVector<dim> const& pQVec) const {
-      return QuantityVector<dim>(eVector - pQVec.eVector);
-    }
-
-    template <typename ScalarDim>
-    auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
-      using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
-
-      if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
-                                                              // a "Quantity" anymore
-      {
-        return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
-      } else {
-        return QuantityVector<typename ResQuantity::dimension_type>(eVector *
-                                                                    p.magnitude());
-      }
-    }
-
-    template <typename ScalarDim>
-    auto operator/(phys::units::quantity<ScalarDim, double> const p) const {
-      return (*this) * (1 / p);
-    }
-
-    auto operator*(double const p) const { return QuantityVector<dim>(eVector * p); }
-
-    auto operator/(double const p) const { return QuantityVector<dim>(eVector / p); }
-
-    auto& operator/=(double const p) {
-      eVector /= p;
-      return *this;
-    }
-
-    auto& operator*=(double const p) {
-      eVector *= p;
-      return *this;
-    }
-
-    auto& operator+=(QuantityVector<dim> const& pQVec) {
-      eVector += pQVec.eVector;
-      return *this;
-    }
-
-    auto& operator-=(QuantityVector<dim> const& pQVec) {
-      eVector -= pQVec.eVector;
-      return *this;
-    }
-
-    auto& operator-() const { return QuantityVector<dim>(-eVector); }
-
-    auto normalized() const { return (*this) * (1 / norm()); }
-
-    auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
-  };
-
-} // namespace corsika::geometry
-
-template <typename dim>
-auto& operator<<(std::ostream& os, corsika::geometry::QuantityVector<dim> qv) {
-  using Quantity = phys::units::quantity<dim, double>;
-
-  os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
-     << phys::units::to_unit_symbol<dim, double>(
-            Quantity(phys::units::detail::magnitude_tag, 1));
-  return os;
-}
-
-#endif

Framework/Geometry/RootCoordinateSystem.h

-#ifndef _include_corsika_geometry_rootcoordinatesystem_h_
-#define _include_corsika_geometry_rootcoordinatesystem_h_
-
-#include <corsika/utl/Singleton.h>
-
-#include <corsika/geometry/CoordinateSystem.h>
-
-/*!
- * This is the only way to get a root-coordinate system, and it is a
- * singleton. All other CoordinateSystems must be relative to the
- * RootCoordinateSystem
- */
-
-namespace corsika::geometry {
-
-  class RootCoordinateSystem : public corsika::utl::Singleton<RootCoordinateSystem> {
-    
-    friend class corsika::utl::Singleton<RootCoordinateSystem>;
-    
-  protected:
-    
-    RootCoordinateSystem() {}
-    
-  public:
-
-    corsika::geometry::CoordinateSystem& GetRootCS() { return fRootCS; }
-    const corsika::geometry::CoordinateSystem& GetRootCS() const { return fRootCS; }
-
-  private:
-    corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT
-    
-  };
-  
-}
-
-#endif

Framework/Geometry/Sphere.h

-
-/**
- * (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_SPHERE_H_
-#define _include_SPHERE_H_
-
-#include <corsika/geometry/Point.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  class Sphere {
-    Point center;
-    LengthType const radius;
-
-  public:
-    Sphere(Point const& pCenter, LengthType const pRadius)
-        : center(pCenter)
-        , radius(pRadius) {}
-
-    //! returns true if the Point \a p is within the sphere
-    auto Contains(Point const& p) const {
-      return radius * radius > (center - p).squaredNorm();
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Trajectory.h

-
-/**
- * (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_TRAJECTORY_H
-#define _include_TRAJECTORY_H
-
-#include <corsika/units/PhysicalUnits.h>
-
-using corsika::units::si::LengthType;
-using corsika::units::si::TimeType;
-
-namespace corsika::geometry {
-
-  template <typename T>
-  class Trajectory : public T {
-
-    corsika::units::si::TimeType fTimeLength; 
-
-  public:
-    using T::GetPosition;
-    using T::GetDistanceBetween;
-
-    Trajectory(T const& theT, 
-               corsika::units::si::TimeType timeLength)
-        : T(theT)
-        , fTimeLength(timeLength) {}
-
-    /*Point GetPosition(corsika::units::si::TimeType t) const {
-      return fTraj.GetPosition(t + fTStart);
-      }*/
-    
-    Point GetPosition(const double u) const {
-      return T::GetPosition(fTimeLength * u);
-    }
-
-    TimeType GetDuration() const { return fTimeLength; }
-    
-    LengthType GetDistance(const corsika::units::si::TimeType t) const {
-      assert(t>fTimeLength);
-      assert(t>=0*corsika::units::si::second);
-      return T::DistanceBetween(0, t);
-    }
-
-  };
-
-} // namespace corsika::geometry
-
-#endif