Documentation

[HealthyPear]

The only docs we have its the README and features like that introduced by #21 are not easy to understand if you are not an experienced developer.

What if we strat proper documentation like in other CTA projects like ctapipe/protopipe?

[mathewpotts]

I have a question about using this. I am using the CorsikaParticleFile and I don't understand the particle_description parameter. They don't seem to correspond with particle IDs I know. For reference, I want to look at all the muons from an event.

[maxnoe]

Hi @mathewpotts

The particle_description is what is written by CORSIKA, you can find it's definition in the CORSIKA manual, Table 10 on page 135 (page number on the page) / 147 (page number of the PDF).
https://www.iap.kit.edu/corsika/downloads/CORSIKA_GUIDE7.7550.pdf

It's a composite value of the particle id and other information:

particle description encoded as:
part. id×1000 + hadr. generation × 10 + no. of obs. level
[for additional muon information or id 95/96:
part. id×1000 + hadr. generation

So to get the particle id, you need to compute:

particle_description // 1000

[maxnoe]

The scikit-hep particle package can help you with this a bit:

from particle import Particle, Corsika7ID

In [8]: Corsika7ID.from_pdgid(Particle.from_name("mu-").pdgid)
Out[8]: <Corsika7ID: 6>

In [9]: Corsika7ID.from_pdgid(Particle.from_name("mu+").pdgid)
Out[9]: <Corsika7ID: 5>

https://github.com/scikit-hep/particle

Small example:

import numpy as np
from corsikaio import CorsikaParticleFile
from particle import Corsika7ID, Particle


MUON_IDS = [
    # get corsika 7 id from the particle name
    Corsika7ID.from_pdgid(Particle.from_name("mu-").pdgid), 
    Corsika7ID.from_pdgid(Particle.from_name("mu+").pdgid), 
]

with CorsikaParticleFile("./DAT000006") as f:
    for e in f:
        particle_ids = e.particles["particle_description"].astype(int) // 1000

        # boolean mask for all muons
        muons = np.isin(particle_ids, MUON_IDS)

        # count them
        ids, counts = np.unique(particle_ids[muons], return_counts=True)

        for cid, count in zip(ids, counts):
            # create particle from the corsika 7 id
            p = Particle.from_pdgid(Corsika7ID(cid).to_pdgid())
            print(p, count)