Units of azimuth and zenith in event header are wrong, should be rad

[VictorBarbosaMartins]

After discussions with @orelgueta , we have decided to open this issue due to 2 reasons:

1 - It seems that the units for azimuth and zenith are not correct here (they are actually given in radians):

pycorsikaio/corsikaio/subblocks/event_header.py

Lines 17 to 18 in ebaefd5

	Field(11, "zenith", unit="deg"),
	Field(12, "azimuth", unit="deg"),

I came across this by testing the implementation of the corsika output module in gammasim-tools. I suggest to change the two units to "rad".

2 - It would be very useful for us in the gammasim-tools team to have the units of the Fields in the event_header and run_header as astropy.units. Given that some of the "units" are non-trivial (e.g. GeV/c), I had to implement a quick solution that is working for us in gammasim-tools and that I copy it here as a suggestion or starting point:

def parse_astropy_unit(unit_in_string):
    """
    Convert units passed in form of a string to astropy.units.Unity.

    Parameters
    ----------
    unit_in_string: str
        Unit to be parsed, e.g. GeV/c
    """

    try:
        astropy_unit = u.Unit(unit_in_string)
    except ValueError:
        # Define the possible math operatos present in the string
        ops_symbol = {
            "+": operator.add,
            "-": operator.sub,
            "/": operator.truediv,
            "*": operator.mul,
            "^": operator.pow,
        }
        # Define the initial unit as dimensionless
        astropy_unit = u.dimensionless_unscaled
        for symbol in ops_symbol:
            # Split the given string according to the math operators
            splitted_unit_string = unit_in_string.split(symbol)
            # If it was indeed split into more than 1 piece, we handle it below
            if len(splitted_unit_string) > 1:
                for unit in splitted_unit_string:
                    # For now, `c` is the only constant that needed this manual adjustment.
                    if unit == "c":
                        unit = const.c
                    # For each unit element we use the current operator to bring it to the final
                    # unit.
                    astropy_unit = ops_symbol[symbol](astropy_unit, u.Unit(unit))
    return astropy_unit

Thank you.

[HealthyPear]

Hi @VictorBarbosaMartins,

I was also interested in this as I am working with CORSIKA and astropy in SWGO.

I developed and CorsikaReader which gives back an astropy table and I started #28 as a result.

It is not at all concluded, just the version I called "stable" in SWGO, but I was hoping that could go on if also you and the members of the gamma-sim group are interested in this.

[maxnoe]

2 - It would be very useful for us in the gammasim-tools team to have the units of the Fields in the event_header and run_header as astropy.units. Given that some of the "units" are non-trivial (e.g. GeV/c), I had to implement a quick solution that is working for us in gammasim-tools and that I copy it here as a suggestion or starting point:

There is a much simpler solution for this part, c as constant speed of line unit is defined as part of the cds subpackage in astropy (it's a bit arbitrary, but works):

In [1]: from astropy.units import cds

In [2]: import astropy.units as u

In [3]: with cds.enable():
   ...:     unit = u.Unit("GeV/c")
   ...: 

In [4]: unit
Out[4]: Unit("GeV / c")

Alternatively, you can also define c as a custom unit yourself, without having to enable the full cds units package:

In [1]: from astropy.constants import c as speed_of_light

In [2]: import astropy.units as u

In [3]: c_unit = u.def_unit("c", speed_of_light)

In [4]: u.add_enabled_units(c_unit)
Out[4]: <astropy.units.core._UnitContext at 0x7f76d46dc730>

In [5]: u.Unit("GeV / c")
Out[5]: Unit("GeV / c")

In [6]: (1 * u.Unit("GeV / c")).to(u.N * u.s)
Out[6]: <Quantity 5.34428599e-19 N s>

[maxnoe]

For the unit issue, could you just open a PR fixing the units? I.e. changing "deg" to "rad", I just checked with the user guide and you are right, the units are wrong.

[maxnoe]

I opened a PR

[VictorBarbosaMartins]

Hi @HealthyPear, thanks a lot for pointing me to this. It is interesting, and I am gonna have a look at your implementation. At the moment, in gammasim-tools we are only looking into the IACT output file but it might be interesting in the future to read and get astropy tables from the other output files as well. I will keep an eye on this PR.

[VictorBarbosaMartins]

2 - It would be very useful for us in the gammasim-tools team to have the units of the Fields in the event_header and run_header as astropy.units. Given that some of the "units" are non-trivial (e.g. GeV/c), I had to implement a quick solution that is working for us in gammasim-tools and that I copy it here as a suggestion or starting point:

There is a much simpler solution for this part, c as constant speed of line unit is defined as part of the cds subpackage in astropy (it's a bit arbitrary, but works):

In [1]: from astropy.units import cds

In [2]: import astropy.units as u

In [3]: with cds.enable():
   ...:     unit = u.Unit("GeV/c")
   ...: 

In [4]: unit
Out[4]: Unit("GeV / c")

Alternatively, you can also define c as a custom unit yourself, without having to enable the full cds units package:

In [1]: from astropy.constants import c as speed_of_light

In [2]: import astropy.units as u

In [3]: c_unit = u.def_unit("c", speed_of_light)

In [4]: u.add_enabled_units(c_unit)
Out[4]: <astropy.units.core._UnitContext at 0x7f76d46dc730>

In [5]: u.Unit("GeV / c")
Out[5]: Unit("GeV / c")

In [6]: (1 * u.Unit("GeV / c")).to(u.N * u.s)
Out[6]: <Quantity 5.34428599e-19 N s>

Hi Max, thanks a lot. That is indeed much nicer and simpler. I will go for your first suggestion and update my solution in gammasim-tools accordingly.

[VictorBarbosaMartins]

I opened a PR

Sorry for the late reply and thanks for updating it!