Modifications for CORSIKA, the IACT interface, and sim_telarray to enable tracing back Cherenkov photons that reached the camera to their emission points, and to more easily correlate them with particles that reach the ground. Designed for use with ctapipe.

Requirements:

• CORSIKA 7.7500.
• IACT/ATMO (bernlohr) 1.67.
• sim_telarray 2024-02-28.
• pyeventio at least v1.15.0 (2025-03-06).
• ctapipe at least 0.23.1.

• Instead of writing MASS and CHARGE in the emitter table, write PARTICLE_ID and GENERATION. These can be used for matches to the particle table.
• Instead of writing X and Y in the emitter table, write XEM and YEM. This is the emission position in the shower frame.
• Write an additional MC_PHOTONS block out from sim_telarray that records the Cherenkov photons that survived raytracing to the focal plane (plane of the camera). This block includes the emitter information of each photon if available.

A reminder of how exactly data gets from CORSIKA to a simtel file:

CORSIKA -> IACT interface -> eventio-format binary stream -> multipipe_corsika ---> write to file
                                                                                \-> sim_telarray with telescope class A -> eventio-format to file
                                                                                \-> sim_telarray with telescope class B -> eventio-format to file
                                                                                \-> ...

The "IACT interface" is a set of C routines called from Fortran, which convert data into eventio format and pass it on to (generally) multipipe_corsika.

multipipe_corsika takes the eventio-format binary stream from the IACT interface and splits it between multiple output locations, including potentially saving directly to disk and being passed to sim_telarray to simulate different types of telescopes or NSB levels from the same air shower.

sim_telarray writes data in eventio format to simtel files, but it doesn't necessarily pass through all of the eventio data that it originally received: some of it is reformatted into different data structures, and some is dropped.

There are two "stages" of configuration of CORSIKA/IACT relevant to these tools: CORSIKA compile flags and CORSIKA steering commands (keywords) issued in the input card.

CORSIKA compile flags can be managed by the build_all script distributed with sim_telarray. This is recommended as the easiest way (see step 2, below).

Useful references:

• Chapters 3-4 in the CORSIKA manual.
• Chapter 3 in the sim_telarray manual.

Follow these steps:

1. Patch: From this project directory, run:

   patch --backup /path/to/install/dir/corsika-77500/src/corsika.F patches/corsika.F.patch
   patch --backup /path/to/install/dir/corsika-77500/bernlohr/iact.c patches/iact.c.patch
   patch --backup /path/to/install/dir/sim_telarray/common/sim_telarray.c patches/sim_telarray.c.patch

2. Compile: From the install directory of CORSIKA/sim_telarray, run:

   TAR_KEEP_NEWER=1 ./build_all prod6 qgs2 iactext muprod store-emitter

   (assuming you want CTAO Prod6 and QGSJET-II).

   The effect of each of these compilation flags is as follows:

   • IACTEXT: Pass photon "emitter" information from CORSIKA into the IACT interface.
   • MUPROD: Include "fated" muons (that don't survive to the observation level) in the particle table passed to the IACT interface.
   • STORE_EMITTER: Set the default value of the IACT STORE-EMITTER keyword to true.

   Fated muons are called "decaying" muons in official documentation. We use this more general term because the muons could have been removed from the particle stack for other reasons besides decay (recorded by the fate index).

3. Simulate: Add your selection of the following commands to your CORSIKA input card:

   • MUADDI YES: Include "birth" muon entries (point at which a muon was created) in the particle table passed to the IACT interface.
   • IACT STORE-EMITTER YES: Pass photon "emitter" information from the IACT interface into the eventio output stream. Defaults to YES if sim_telarray was compiled with the STORE_EMITTER flag.
   • IACT STORE-PARTICLES YES: Pass the particle table from the IACT interface into the eventio output stream.

   Another relevant command is OBSLEV which defines a layer (observation level) at which to record particle information. There can be up to 10 observation levels, listed in any order. The lowest observation level is taken to be the ground level that the telescopes are on. The counting of the level numbers (only relevant to the particle table) is such that level 1 is the highest level and level N is the lowest level.

   You can use the demo notebook to test your created simtel file.

4. Analyse: The Python module cherentrace provides two functions for accessing the photon and particle tables in a convenient format. These functions take care of transforming from CORSIKA units (cm) to ctapipe default units (m) and deriving other useful quantities.

   • import cherentrace
     df_photons = cherentrace.get_photons(source, event, tel_id, to_telescope_frame = True)

     Get a Pandas DataFrame of the Cherenkov photons that reached the focal plane, including emitter information of each photon if available.

   • import cherentrace
     df_particles = cherentrace.get_particles(source, event)

     Get a Pandas DataFrame of the particles that reached each observation level, including additional information about muons if available.

Note that the positions and vector components (cx/cy/cz) are reported in the CORSIKA coordinate system: x is positive towards magnetic north, y is positive towards west, and z is positive upwards. Azimuth is measured anti-clockwise from the x axis (north). The exception to this statement is the az column included in the photon table when to_telescope_frame=True: this is reported in astropy convention.

Astropy uses an azimuth convention of clockwise from north, so you must invert CORSIKA azimuths any time you pass them through astropy tools.

References:

• Table 4 in the CORSIKA manual
• The MUPROD option manual

The emitter information in the photon table only identifies muons as ID 5 or 6 (μ+ and μ-, respectively). In the particle table, there may be additional muon information with other IDs.

The following table lists the meaning of each muon-related particle ID in the particle table.

The following table lists which types of muon entries are written to the output for each combination of compilation and steering options.

In other words, MUPROD enables information about muons that don't survive to the next observation level. MUADDI YES enables information about muon birth points.

Reference:

• Section 3.1 in the MUPROD option manual

The generation counter tracks the interaction history of each particle by adding specific units together to represent different types of interactions. In the photon table we get the full generation counter, however in the particle table the generation counter is only the last two or three digits of the generation counter (two digits if an observation level is given, three digits if the observation level is -1).

For a reference on the full generation counter, refer to Section 3.5.15 in the CORSIKA manual. Although this section is describing the EHISTORY option, this is still the relevant description of the full generation counter that we get in the photon table.

For a reference on how hadronic interactions are handled in the generation counter, refer to Section 4 in the MUPROD option manual.

Note that the last three digits of the generation counter are for the hadronic interactions, so in the particle table (which only lists the last two or three digits), only hadronic interactions are seen.