add sophia files

corsika/detail/modules/sophia/ParticleConversion.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+//#include <sibyll2.3d.hpp>
+
+namespace corsika::sophia {
+
+  inline HEPMassType getSophiaMass(Code const pCode) {
+    if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
+    auto sCode = convertToSophiaRaw(pCode);
+    if (sCode == 0)
+      throw std::runtime_error("getSophiaMass: unknown particle!");
+    else
+      return sqrt(get_sophia_mass2(sCode)) * 1_GeV;
+  }
+} // namespace corsika::sophia
\ No newline at end of file

corsika/modules/Sophia.hpp

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/modules/sophia/ParticleConversion.hpp>
+
+#include <corsika/modules/sophia/InteractionModel.hpp>
+
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+/**
+ * @file Sophia.hpp
+ *
+ * Includes all the parts of the Sophia model. Defines the InteractionProcess<TModel>
+ * classes needed for the ProcessSequence.
+ */
+
+namespace corsika::sophia {
+  /**
+   * @brief sophia::Interaction is the process for ProcessSequence.
+   *
+   * The sophia::InteractionModel is wrapped as an InteractionProcess here in order
+   * to provide all the functions for ProcessSequence.
+   */
+  struct Interaction : public InteractionModel, public InteractionProcess<Interaction> {
+    template <typename TEnvironment>
+    Interaction(TEnvironment const& env)
+        : InteractionModel{env} {}
+  };
+
+} // namespace corsika::sophia

corsika/modules/sophia/ParticleConversion.hpp

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+//#include <sophia.hpp>
+
+#include <string>
+
+namespace corsika::sophia {
+
+  enum class SophiaCode : int8_t;
+  using SophiaCodeIntType = std::underlying_type<SophiaCode>::type;
+
+  // /**
+  //    These are the possible projectile for which Sibyll knows the cross section
+  //  */
+  // enum class SibyllXSClass : int8_t {
+  //   CannotInteract = 0,
+  //   Baryon = 1,
+  //   Pion = 2,
+  //   Kaon = 3,
+  // };
+  // using SibyllXSClassIntType = std::underlying_type<SibyllXSClass>::type;
+
+#include <corsika/modules/sophia/Generated.inc>
+
+  SophiaCode constexpr convertToSophia(Code const pCode) {
+    return corsika2sophia[static_cast<CodeIntType>(pCode)];
+  }
+
+  Code constexpr convertFromSophia(SophiaCode const pCode) {
+    auto const s = static_cast<SophiaCodeIntType>(pCode);
+    auto const corsikaCode = sophia2corsika[s - minSophia];
+    if (corsikaCode == Code::Unknown) {
+      throw std::runtime_error(std::string("SOPHIA/CORSIKA conversion of ")
+                                   .append(std::to_string(s))
+                                   .append(" impossible"));
+    }
+    return corsikaCode;
+  }
+
+  int constexpr convertToSophiaRaw(Code const code) {
+    return static_cast<int>(convertToSophia(code));
+  }
+
+  // int constexpr getSibyllXSCode(Code const code) {
+  //   if (is_nucleus(code))
+  //     return static_cast<SibyllXSClassIntType>(SibyllXSClass::CannotInteract);
+  //   return static_cast<SibyllXSClassIntType>(
+  //       corsika2sibyllXStype[static_cast<CodeIntType>(code)]);
+  // }
+
+  bool constexpr canInteract(Code const pCode) {
+    return (pCode == Code::Photon ? True : False);
+  }
+
+  HEPMassType getSophiaMass(Code const);
+
+} // namespace corsika::sophia
+
+#include <corsika/detail/modules/sophia/ParticleConversion.inl>

modules/sophia/CMakeLists.txt

+set (
+  MODEL_SOURCES
+#  sibyll2.3d.cpp
+  eventgen.f
+  inpoutput.f
+  jetset74dp.f
+  sampling.f
+  SOPHIA20.f
+  )
+
+#set (
+#  MODEL_HEADERS
+#  sibyll2.3d.hpp
+#  nuclib.hpp
+#  )
+
+enable_language (Fortran)
+add_library (Sophia_static STATIC ${MODEL_SOURCES})
+add_library (Sophia SHARED ${MODEL_SOURCES})
+
+set_target_properties (
+  Sophia_static
+  PROPERTIES
+  POSITION_INDEPENDENT_CODE 1
+  )
+
+target_include_directories (
+  Sophia_static
+  PUBLIC
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+  $<INSTALL_INTERFACE:include/corsika_modules/sophia>
+  )
+  target_include_directories (
+    Sophia
+    PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+    $<INSTALL_INTERFACE:include/corsika_modules/sophia>
+    )
+
+target_link_libraries (
+  Sophia_static
+  PUBLIC
+  gfortran
+  )
+target_link_libraries (
+  Sophia
+  PUBLIC
+  gfortran
+  )
+  
+install (
+  FILES
+  ${MODEL_HEADERS}
+  DESTINATION include/corsika_modules/sophia
+  )
+
+install (
+  TARGETS Sophia_static Sophia
+  EXPORT CORSIKA8PublicTargets
+  ARCHIVE DESTINATION lib/corsika
+  LIBRARY DESTINATION lib/corsika # just for cmake 3.10.x (ubuntu 18)
+  )
+
+# add sophia to corsika8 build
+add_dependencies (CORSIKA8 Sophia_static)
+target_link_libraries (CORSIKA8 INTERFACE Sophia_static)

modules/sophia/README.TXT

+
+
+                           SOPHIA 2.0
+                           ==========
+
+
+
+A.M"ucke, Ralph Engel, J.P.Rachen, R.J.Protheroe, and Todor Stanev
+
+(astro-ph/9903478, to appear in Comp.Phys.Commun.)
+
+
+Files:
+-----
+
+
+SOPHIA20.f           (main program)
+
+eventgen.f           (hadronic event generator)
+
+sampling.f           (sampling of photon energies)
+
+inpoutput.f          (histogramming, input/output)
+
+jetset74dp.f         (Lund fragmentation program, changed to double prec.)
+
+rndm.f               (random number generator)
+
+testPL.inp            (sample input for straight power law calculation)
+
+testBP.inp            (sample input for broken power law calculation)
+
+testBB.inp            (sample input for black body radiation calculation)
+
+README.TXT           (this file)
+
+
+
+
+
+Particle record:
+---------------
+
+
+After event generation the final state particles are found in 
+
+       COMMON /S_PLIST/ P(2000,5), LLIST(2000), NP, Ideb
+
+NP is the number of particles produced in the last event.
+The particle momenta of particle I in x, y, z directions are 
+P(I,1), P(I,2), P(I,3) and P(I,4) is the energy (in GeV). 
+The particle mass is stored in P(I,5).
+
+The particle identity is given by LLIST with the following numbering
+scheme:
+
+ -----------------------
+ code particle     mass
+ -----------------------
+   1  gam          0.0000
+   2  e+           0.0005
+   3  e-           0.0005
+   4  mu+          0.1057
+   5  mu-          0.1057
+   6  pi0          0.1350
+   7  pi+          0.1396
+   8  pi-          0.1396
+   9  k+           0.4937
+  10  k-           0.4937
+  11  k0l          0.4977
+  12  k0s          0.4977
+  13  p            0.9383
+  14  n            0.9396
+  15  nue          0.0000
+  16  nueb         0.0000
+  17  num          0.0000
+  18  numb         0.0000
+  21  k0           0.4977
+  22  k0b          0.4977
+  23  eta          0.5488
+  24  etap         0.9576
+  25  rho+         0.7714
+  26  rho-         0.7714
+  27  rho0         0.7717
+  28  k*+          0.8921
+  29  k*-          0.8921
+  30  k*0          0.8965
+  31  k*0b         0.8965
+  32  omeg         0.7826
+  33  phi          1.1020
+  34  SIG+         1.1894
+  35  SIG0         1.1925
+  36  SIG-         1.1973
+  37  XI0          1.3149
+  38  XI-          1.3213
+  39  LAM          1.1156
+  40  DELT++       1.2300
+  41  DELT+        1.2310
+  42  DELT0        1.2320
+  43  DELT-        1.2330
+  44  SIG*+        1.3828
+  45  SIG*0        1.3837
+  46  SIG*-        1.3872
+  47  XI*0         1.5318
+  48  XI*-         1.5350
+  49  OME*-        1.6724
+
+Antibaryons have negative ID numbers correspondingly.
+Decayed particles are marked by adding 10000 to their ID code.
+
+The decay of a particle with code ID can be turned off via
+
+      IDB(ID) = -ABS(IDB(ID))
+
+and turned on via
+
+      IDB(ID) = ABS(IDB(ID))
+
+This has to be done only once at the beginning of event generation or
+might be changed on an event-by-event basis.
+
+
+
+The particles  produced in the last EVENTGEN call can be listed 
+conveniently by 
+      CALL print_event(1)
+
+
+
+
+Related publications:
+--------------------
+
+- M"ucke A., et al. 1999, astro-ph/9903478, to appear in Comp.Phys.Commun.
+
+- M"ucke A., et al. 1999, astro-ph/9808279, to appear in PASA.
+
+- M"ucke A., et al. 1999, to appear in: Proc. of the
+  19th Texas Symposium on Relativistic Astrophysics, Paris,
+  France, Dec. 1998. Eds.: J.~Paul, T.~Montmerle \& E.~Aubourg 
+  (CEA Saclay)
+
+- M"ucke A., et al. 1999, astro-ph/9905153, to appear in: Proc. of
+  19th Texas Symposium on Relativistic Astrophysics, Paris, France,
+  Dec. 1998. Eds.: J.~Paul, T.~Montmerle \& E.~Aubourg 
+  (CEA Saclay)
+
+- M"ucke A., et al 1999, to appear in: Proc. of 26th Int.Cosmic
+  Ray Conf. (Salt Lake City, Utah) 

modules/sophia/SOPHIA20.f

+c*****************************************************************************
+c**!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!***
+c**!!              IF YOU USE THIS PROGRAM, PLEASE CITE:                 !!***
+c**!! A.M"ucke, Ralph Engel, J.P.Rachen, R.J.Protheroe and Todor Stanev, !!***
+c**!!  1999, astro-ph/9903478, to appear in Comp.Phys.Commun.            !!***
+c**!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!***
+c*****************************************************************************
+c** Further SOPHIA related papers:                                         ***
+c** (1) M"ucke A., et al 1999, astro-ph/9808279, to appear in PASA.        ***
+c** (2) M"ucke A., et al 1999, to appear in: Proc. of the                  ***
+c**      19th Texas Symposium on Relativistic Astrophysics, Paris, France, ***
+c**      Dec. 1998. Eds.: J.~Paul, T.~Montmerle \& E.~Aubourg (CEA Saclay) ***
+c** (3) M"ucke A., et al 1999, astro-ph/9905153, to appear in: Proc. of    ***
+c**      19th Texas Symposium on Relativistic Astrophysics, Paris, France, ***
+c**      Dec. 1998. Eds.: J.~Paul, T.~Montmerle \& E.~Aubourg (CEA Saclay) ***
+c** (4) M"ucke A., et al 1999, to appear in: Proc. of 26th Int.Cosmic Ray  ***
+c**      Conf. (Salt Lake City, Utah)                                      ***
+c*****************************************************************************
+
+
+
+       program SOPHIA20
+
+c************************************************************************
+c*** Simulation Of PhotoHadronic Interactions in Astrophysics ***********
+c***           VERSION 2.0                                    ***********
+c************************************************************************
+
+c************************************************************************
+c** Main program for photopion production of relativistic nucleons     **
+c**  in a radiation field (blackbody or power law)                     **
+c************************************************************************
+c** Date: 20/01/98          **
+c** correct.:19/02/98       **
+c** first release to        ** 
+c**  collab.:      25/05/98 **
+c** Version 1.3:   31/08/98 **
+c** Version 1.4:   12/10/98 **
+c** change to DP:  16/11/98 **
+c** last corr.:       07/99 **
+c** authors: A.G.F. Muecke  **
+c**          R.R. Engel     **
+c**   in collaboration with **
+c**         R.J. Protheroe  **
+c**         J.P. Rachen     **
+c**         T.S. Stanev     **
+c*****************************
+
+c****** INPUT ***********************************************************************
+c E0 = energy of incident proton (in lab frame) [in GeV]
+c L0 = code number of the incident nucleon (L0=13: proton, L0=14: neutron)
+c soft photon field:   for thermal spectrum: temperature tbb [in K]
+c                      for PL: (set tbb = 0) alpha = PL index (S \sim nu^{-alpha}), 
+c                                            epsmin = low energy cut off  [in eV]
+c                                            epsmax = high energy cut off [in eV]
+c****** OUTPUT **********************************************************************
+c** energy distribution P(x) (logarithmic scale, x=E_particle/E0) for ***
+c** photons, protons, neutrons, e-neutrinos, nu-neutrinos             ***
+c************************************************************************
+       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+       IMPLICIT INTEGER (I-N)
+       SAVE
+
+       COMMON/input/ tbb,E0,alpha1,alpha2,
+     &           epsm1,epsm2,epsb,L0
+       COMMON /S_PLIST/ P(2000,5), LLIST(2000), NP, Ideb
+       COMMON /S_MASS1/ AM(49), AM2(49)
+       COMMON /S_CHP/  S_LIFE(49), ICHP(49), ISTR(49), IBAR(49)
+       COMMON /S_CSYDEC/ CBR(102), IDB(49), KDEC(612), LBARP(49)
+
+      CHARACTER*6 NAMPRES
+      COMMON /RES_PROP/ AMRES(9), SIG0(9),WIDTH(9), 
+     +                    NAMPRES(0:9)
+
+      CHARACTER*6 NAMPRESp
+      COMMON /RES_PROPp/ AMRESp(9), BGAMMAp(9),WIDTHp(9),  
+     +                    RATIOJp(9),NAMPRESp(0:9)
+
+      CHARACTER*6 NAMPRESn
+      COMMON /RES_PROPn/ AMRESn(9), BGAMMAn(9),WIDTHn(9),  
+     +                    RATIOJn(9),NAMPRESn(0:9)
+
+
+       DIMENSION DNg(201),DNnum(201),DNnue(201),DNp(201),DNn(201)
+       DIMENSION DNnuma(201),DNnuea(201),E0_arr(101)
+       DIMENSION Dg(101,201),Dnum(101,201),Dnue(101,201)
+       DIMENSION Dp(101,201),Dn(101,201),Dnuma(101,201),Dnuea(101,201)
+       DIMENSION Dem(101,201),Dep(101,201),DNem(201),DNep(201)
+
+       dimension c_feps(10000),eps_arr(10000)
+
+       character*6 nameinc
+       character*1 ans
+
+       external sample_eps,sample_s,eventgen,
+     &  listdistr,output,initial,prob_epskt
+
+       DATA pi /3.141593D0/ 
+
+        do i=1,201
+          DNg(i) = 0.D0
+          DNnum(i) = 0.D0
+          DNnue(i) = 0.D0
+          DNp(i) = 0.D0
+          DNn(i) = 0.D0
+          DNem(i) = 0.D0
+          DNep(i) = 0.D0
+          DNnuma(i) = 0.D0
+          DNnuea(i) = 0.D0
+
+          do jm=1,101
+            Dg(jm,i) = 0.D0
+            Dnum(jm,i) = 0.D0
+            Dnue(jm,i) = 0.D0
+            Dp(jm,i) = 0.D0
+            Dn(jm,i) = 0.D0
+            Dem(jm,i) = 0.D0
+            Dep(jm,i) = 0.D0
+            Dnuma(jm,i) = 0.D0
+            Dnuea(jm,i) = 0.D0
+         enddo
+
+       enddo
+
+
+c****** INPUT **************************************************
+       print*
+       print*,'Give in code number of incident nucleon: '
+       print*,' (13 = proton, 14 = neutron)'
+       read(*,*) L0
+
+       delE = 0.D0
+       print*
+       print*,'Incident nucleon spectrum:'
+       print*,'energy grid [s] or single nucleon [n] ? '
+       read(*,'(A1)') ans
+
+       if (ans.eq.'n') then
+        jm = 1
+        ninc = 1
+ 4     print*,'Give energy [in GeV] of incident nucleon: '
+       read(*,*) E0
+
+       if (E0.lt.AM(L0)) then
+        print*,'No valid input !'
+        goto 4
+       endif
+       Emin = log10(E0)
+       Emax = Emin
+       endif
+       if (ans.eq.'s') then
+  5    print*,'Give in low-energy cut off (in GeV) of nucleon 
+     &  energy grid:'
+       read(*,*) Emin
+
+       if (Emin.lt.AM(L0)) then
+        print*,'No valid input !'
+        goto 5
+       endif
+       Emin = log10(Emin)
+       print*,'Give in high-energy cut off (in GeV) of nucleon 
+     &  energy grid:'
+       read(*,*) Emax
+
+       if (Emax.lt.Emin) then
+        print*,'No valid input !'
+        goto 5
+       endif
+       Emax = log10(Emax)
+       print*,'Give number of bins (< 101):'
+       read(*,*) ninc
+       delE = (Emax-Emin)/ninc
+       ninc = ninc+1
+       endif
+
+       print*
+  2    print*,'Give soft photon spectrum: '
+       print*,' blackbody spectrum ? (y/n): '
+       read(*,'(A1)') ans
+
+       if (ans.eq.'y') then
+         rad_den = 1.D0
+         print*,'Give temperature [in K]: '
+         read(*,*) tbb
+
+       else if (ans.eq.'n') then
+       tbb = 0.D0
+       print*,' for power law spectrum n(eps) ~ eps^-alpha:'
+       print*,'  (alpha = 0 ... 3)'
+  8    print*,'straight power law ? (y/n): '
+
+       read(*,'(A1)') ans
+       if (ans.eq.'y') then
+  3    print*,'   spectral index alpha = '
+       read(*,*) alpha2
+
+        if (alpha2.lt.0.D0.or.alpha2.gt.3.D0) then
+         print*,'no valid input !'
+         goto 3
+        endif
+       alpha1 = alpha2
+       else if (ans.eq.'n') then
+        print*,'broken power law:'
+        print*,' spectral index alpha1, alpha2 = '
+        read(*,*) alpha1,alpha2
+        if (alpha1.lt.0.D0.or.alpha1.gt.3.D0.or.
+     &  alpha2.lt.0.D0.or.alpha2.gt.3.D0) then
+         print*,'no valid input !'
+         goto 8
+        endif
+        print*,'break energy [eV] = '
+        read(*,*) epsb        
+       else
+         print*,'no valid input !'
+         goto 8
+       endif 
+       print*,'low-energy cut off [eV] = '
+       read(*,*) epsmin
+       print*,'high-energy cut off [eV] = '
+       read(*,*) epsmax
+
+       if (ans.eq.'y') epsb = epsmin/100.D0
+       else
+        print*,'no valid input !'
+        goto 2
+       endif
+      
+       print*
+       print*,'Give number of trials: '
+       read(*,*) ntrial
+
+       print*
+       print*,
+     &'Give number of bins (<201) for output particle spectra: '
+       print*,
+     &'(spectra in logarithmic equal bins, with stepsize Delta x'
+       print*,'  with x=E_particle/E_initial nucleon )'
+       read(*,*) nbins
+
+       print*
+       print*,
+     &'Give stepsize Delta x for output particle spectra: '
+       read(*,*) delx
+
+       print*
+       print*,'Give filename (< 7 CH) for output: '
+       read(*,'(A6)') nameinc
+       print*
+
+       call initial(L0)
+
+c... nucleon energy loop:
+       do jm=1,ninc
+         Elog = Emin+(jm-1)*delE
+         E0 = 10.D0**Elog
+         E0_arr(jm) = E0
+
+c... trial loop:
+       do nt=1,ntrial 
+
+c*******************************************************
+c sample epsilon = energy of photon in lab frame
+c*******************************************************
+       pm = AM(L0)
+  6    call sample_eps(epseV,epsmin,epsmax)
+       if (epseV.le.0.) goto 7
+c *** eps in GeV:
+       eps = epseV/1.D9
+       Etot = E0+eps
+
+c*******************************************************
+c sample s = total mass of center of momentum frame
+c*******************************************************
+ 18    call sample_s(s,eps)
+       gammap = E0/pm
+       betap = sqrt(1.D0-1.D0/gammap/gammap)
+       theta = ((pm*pm-s)/2.D0/E0/eps+1.D0)/betap
+       if (abs(theta).gt.1.D0) STOP
+       theta = acos(theta)*180.D0/pi 
+
+c***********************************************************
+c*** CALL PHOTOPION EVENT GENERATOR 
+c***********************************************************
+
+        call eventgen(L0,E0,eps,theta,Imode)
+
+c*********************************************************
+c store decayed particles and their energy distribution: *
+c*********************************************************
+c... calculate particle distribution dN/dlog(f) with f=E/E0:
+c    [note: E*dN/dE = dN/dlog(f)]
+
+        call listdistr(E0,DNg,DNnum,DNnuma,DNnue,DNnuea,
+     &   DNp,DNn,DNem,DNep,nbins,delx)
+
+        xini = -nbins*delx
+        do 55 i=1,nbins
+         Dg(jm,i) = Dg(jm,i)+DNg(i)/delx/ntrial
+         Dnum(jm,i) = Dnum(jm,i)+DNnum(i)/delx/ntrial
+         Dnue(jm,i) = Dnue(jm,i)+DNnue(i)/delx/ntrial
+         Dnuma(jm,i) = Dnuma(jm,i)+DNnuma(i)/delx/ntrial
+         Dnuea(jm,i) = Dnuea(jm,i)+DNnuea(i)/delx/ntrial
+         Dp(jm,i) = Dp(jm,i)+DNp(i)/delx/ntrial
+         Dn(jm,i) = Dn(jm,i)+DNn(i)/delx/ntrial
+         Dem(jm,i) = Dem(jm,i)+DNem(i)/delx/ntrial
+         Dep(jm,i) = Dep(jm,i)+DNep(i)/delx/ntrial
+  55    continue
+      
+c... trial loop:
+       enddo
+  7    continue
+
+c... nucleon energy loop:
+       enddo
+
+       call output(Dg,Dnum,Dnuma,Dnue,Dnuea,Dp,Dn,Dem,Dep,
+     &  nbins,ninc,nameinc,delx,Emin,Emax,E0_arr,epsmin,epsmax)
+
+       STOP
+       
+
+       END

src/modules/sophia/CMakeLists.txt

+set (input_dir ${PROJECT_SOURCE_DIR}/src/modules/sophia)
+set (output_dir ${PROJECT_BINARY_DIR}/corsika/modules/sophia)
+
+file (MAKE_DIRECTORY ${output_dir})
+
+add_custom_command (
+  OUTPUT  ${output_dir}/Generated.inc
+  COMMAND ${input_dir}/code_generator.py 
+          ${PROJECT_BINARY_DIR}/corsika/framework/core/particle_db.pkl
+          ${input_dir}/sophia_codes.dat
+  DEPENDS ${input_dir}/code_generator.py
+          ${input_dir}/sophia_codes.dat
+          ${PROJECT_BINARY_DIR}/corsika/framework/core/particle_db.pkl
+          GenParticlesHeaders
+  WORKING_DIRECTORY
+          ${output_dir}/
+  COMMENT "Generate conversion tables for particle codes sophia <-> CORSIKA"
+  VERBATIM
+  )
+
+set_source_files_properties (
+  ${output_dir}/Generated.inc
+  PROPERTIES GENERATED TRUE
+)
+
+add_custom_target (SourceDirLinkSoph DEPENDS ${output_dir}/Generated.inc)
+add_dependencies (CORSIKA8 SourceDirLinkSoph)
+
+install (
+  FILES ${output_dir}/Generated.inc
+  DESTINATION include/corsika/modules/sophia
+  )

src/modules/sophia/code_generator.py

+#!/usr/bin/env python3
+
+# (c) Copyright 2018-2019 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.
+
+
+import pickle, sys, itertools
+
+
+
+def load_particledb(filename):
+    '''
+    loads the pickled particle_db (which is an OrderedDict)
+    '''
+    with open(filename, "rb") as f:
+        particle_db = pickle.load(f)
+    return particle_db
+
+
+
+def read_sophia_codes(filename, particle_db):
+    '''
+    reads to sophia codes data file
+
+    For particls known to sophia, add 'sophia_code' and 'sophia_xsType' to particle_db
+    '''
+    with open(filename) as f:
+        for line in f:
+            line = line.strip()
+            if len(line)==0 or line[0] == '#':
+                continue            
+            identifier, sib_code, canInteractFlag, xsType = line.split()
+            try:
+                particle_db[identifier]["sophia_code"] = int(sib_code)
+                particle_db[identifier]["sophia_xsType"] = xsType
+            except KeyError as e:
+                raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
+
+
+            
+
+def generate_sophia_enum(particle_db):
+    '''
+     generates the enum to access sophia particles by readable names
+    '''
+    output = "enum class sophiaCode : int8_t {\n"
+    for identifier, pData in particle_db.items():
+        if 'sophia_code' in pData:
+            output += "  {:s} = {:d},\n".format(identifier, pData['sophia_code'])
+    output += "};\n"
+    return output
+
+
+
+def generate_corsika2sophia(particle_db):    
+    '''
+    generates the look-up table to convert corsika codes to sophia codes
+    '''
+    string = "std::array<sophiaCode, {:d}> constexpr corsika2sophia = {{\n".format(len(particle_db))
+    for identifier, pData in particle_db.items():
+        if pData['isNucleus']: continue
+        if 'sophia_code' in pData:
+            string += "  sophiaCode::{:s}, \n".format(identifier)
+        else:
+            string += "  sophiaCode::Unknown, // {:s}\n".format(identifier + ' not implemented in sophia')
+    string += "};\n"
+    return string
+    
+
+
+def generate_corsika2sophia_xsType(particle_db):    
+    '''
+    generates the look-up table to convert corsika codes to sophia codes
+    '''
+    string = "std::array<sophiaXSClass, {:d}> constexpr corsika2sophiaXStype = {{\n".format(len(particle_db))
+    for identifier, pData in particle_db.items():
+        if pData['isNucleus']: continue
+        if 'sophia_xsType' in pData:
+            string += "  sophiaXSClass::{:s}, // {:s}\n".format(pData['sophia_xsType'], identifier)
+        else:
+            string += "  sophiaXSClass::CannotInteract, // {:s}\n".format(identifier + ' not implemented in sophia')
+    string += "};\n"
+    return string
+
+
+def generate_sophia2corsika(particle_db) :
+    '''
+    generates the look-up table to convert sophia codes to corsika codes    
+    '''
+    string = ""
+    
+    minID = 0
+    for identifier, pData in particle_db.items() :
+        if 'sophia_code' in pData:
+            minID = min(minID, pData['sophia_code'])
+
+    string += "sophiaCodeIntType constexpr minsophia = {:d};\n\n".format(minID)
+
+    pDict = {}
+    for identifier, pData in particle_db.items() :
+        if 'sophia_code' in pData:
+            sib_code = pData['sophia_code'] - minID
+            pDict[sib_code] = identifier
+    
+    nPart = max(pDict.keys()) - min(pDict.keys()) + 1
+    string += "std::array<corsika::Code, {:d}> constexpr sophia2corsika = {{\n".format(nPart)
+    
+    for iPart in range(nPart) :
+        if iPart in pDict:
+            identifier = pDict[iPart]
+        else:
+            identifier = "Unknown"
+        string += "  corsika::Code::{:s}, \n".format(identifier)
+    
+    string += "};\n"
+    return string
+
+if __name__ == "__main__":
+    if len(sys.argv) != 3:
+        print("usage: {:s} <particle_db.pkl> <sophia_codes.dat>".format(sys.argv[0]), file=sys.stderr)
+        sys.exit(1)
+        
+    print("code_generator.py for sophia")
+    
+    particle_db = load_particledb(sys.argv[1])
+    read_sophia_codes(sys.argv[2], particle_db)
+    
+    with open("Generated.inc", "w") as f:
+        print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
+        print(generate_sophia_enum(particle_db), file=f)
+        print(generate_corsika2sophia(particle_db), file=f)
+        print(generate_sophia2corsika(particle_db), file=f)
+        print(generate_corsika2sophia_xsType(particle_db), file=f)

src/modules/sophia/sophia_codes.dat

+# input file for particle conversion to/from sophia
+# the format of this file is: "corsika-identifier" "sophia-id" "can-interact-in-sophia" "cross-section-type"
+
+# The unknown particle is to handle all particles that are not known to sophia.
+# It is important that sophia-id does not overlap with any existing sophia particle!
+# Be careful
+Unknown                  0       0       CannotInteract
+
+# Here is the list of particles known to sophia 
+Electron                 3       0       CannotInteract
+Positron                 2       0       CannotInteract
+NuE                      15      0       CannotInteract
+NuEBar                   16      0       CannotInteract
+MuMinus                  5       0       CannotInteract
+MuPlus                   4       0       CannotInteract
+NuMu                     17      0       CannotInteract
+NuMuBar                  18      0       CannotInteract
+Photon                   1       1       Photon
+Pi0                      6       0       CannotInteract
+#                        rho0    could interact but sophia has no cross section/interaction length. was used for gamma had int
+Rho0                     27      0       CannotInteract
+K0Long                   11      0       CannotInteract
+K0                       21      0       CannotInteract
+K0Bar                    22      0       CannotInteract
+PiPlus                   7       0       CannotInteract
+PiMinus                  8       0       CannotInteract
+RhoPlus                  25      0       CannotInteract
+RhoMinus                 26      0       CannotInteract
+Eta                      23      0       CannotInteract
+EtaPrime                 24      0       CannotInteract
+Omega                    32      0       CannotInteract
+K0Short                  12      0       CannotInteract
+KStar0                   30      0       CannotInteract
+KStar0Bar                31      0       CannotInteract
+KPlus                    9       0       CannotInteract
+KMinus                   10      0       CannotInteract
+KStarPlus                28      0       CannotInteract
+KStarMinus               29      0       CannotInteract
+Neutron                  14      0       CannotInteract
+AntiNeutron              -14     0       CannotInteract
+Delta0                   42      0       CannotInteract
+Delta0Bar                -42     0       CannotInteract
+DeltaMinus               43      0       CannotInteract
+DeltaPlusBar             -43     0       CannotInteract
+Proton                   13      0       CannotInteract
+AntiProton               -13     0       CannotInteract
+DeltaPlus                41      0       CannotInteract
+DeltaMinusBar            -41     0       CannotInteract
+DeltaPlusPlus            40      0       CannotInteract
+DeltaMinusMinusBar       -40     0       CannotInteract
+SigmaMinus               36      0       CannotInteract
+SigmaPlusBar             -36     0       CannotInteract
+SigmaStarMinus           46      0       CannotInteract
+SigmaStarPlusBar         -46     0       CannotInteract
+SigmaStarPlus            44      0       CannotInteract
+SigmaStarMinusBar        -44     0       CannotInteract
+SigmaStar0               45      0       CannotInteract
+SigmaStar0Bar            -45     0       CannotInteract
+Lambda0                  39      0       CannotInteract
+Lambda0Bar               -39     0       CannotInteract
+Sigma0                   35      0       CannotInteract
+Sigma0Bar                -35     0       CannotInteract
+SigmaPlus                34      0       CannotInteract
+SigmaMinusBar            -34     0       CannotInteract
+XiMinus                  38      0       CannotInteract
+XiPlusBar                -38     0       CannotInteract
+Xi0                      37      0       CannotInteract
+Xi0Bar                   -37     0       CannotInteract
+XiStarMinus              48      0       CannotInteract
+XiStarPlusBar            -48     0       CannotInteract  
+XiStar0                  47      0       CannotInteract
+XiStar0Bar               -47     0       CannotInteract
+OmegaMinus               49      0       CannotInteract
+OmegaPlusBar             -49     0       CannotInteract
+Phi                      33      0       CannotInteract
\ No newline at end of file