Specify IRF axis orders

[cdeil]

This pull request specifies the IRF axis orders.

It became apparent that this is needed because we accidentally chose an ENERGY, THETA order for BKG_2D, whereas the current AEFF and PSF formats use THETA, ENERGY (see gammapy/gammapy#461 (comment)).

I think this pull request is uncontroversial and I'll merge it later today.

Still, @jknodlseder, @mimayer, @JouvinLea, if you could have a look and specifically comment on these questions, that would be useful:

• OK to use THETA, ENERGY axis order everywhere? Or is ENERGY, THETA better? The main difference is memory access speed ... values that are typically accessed together should be consecutive in memory so that CPU cache performance is better. I haven't thought about / benchmarked what is better here.
• Does the axis order ENERGY, DETX, DETY correspond to the current Gammalib implementation? Or do you have ENERGY, DETY, DETX?
• It's a bit weird that we have energy first for 3D and last for 2D, but I guess that's OK?
• Any other thoughts?

[jknodlseder]

From the GammaLib perspective, I recently changed the code so that the order is not relevant anymore. I have however not tested this. And so far, I had used the ENERGY, THETA order for producing the CTA IRFs that are bundled with the software.

Personally I prefer the ENERGY, THETA order since (a) it's the order Fermi-LAT uses (hence our format would be inconsistent with their format), and (b) it orders things by importance of dependency (the THETA dependency for the PSF is smaller than the ENERGY dependency for example).

Concerning memory access, what is true for one application is false for another application. You cannot build a general data format on the assumed usage of that format. If a specific code needs things ordered in another way to be more efficient, the code has to deal with that reordering. In any case, GammaLib performs bilinear interpolations, hence needs to access both dimensions in any case.

The current GammaLib implementation is ENERGY, DETX and DETY. This is consistent with the ENERGY, THETA order, as DETX, DETY is just a 2D version of THETA. As argued above, I would always put ENERGY first.

[cdeil]

@jknodlseder - OK. ENERGY, THETA everywhere sounds good to me.
I can adapt the HAP-HD exporter and Gammapy today.

For now, I'm changing this spec update to say "recommended axis order: ENERGY, THETA".
@mimayer - OK?

A better description and example of how axis order is defined in this BINTABLE FITS format would be good here:
https://github.com/open-gamma-ray-astro/gamma-astro-data-formats/blob/master/source/info/fits-arrays.rst#bintable-hdu
@jknodlseder - you probably don't have time to help edit the spec in better shape, right?
(especially the event list spec is still in very bad shape, but really all the parts need to be improved)

[mimayer]

@mimayer - OK?

ok

[cdeil]

There is the case of the RPSF:
http://gamma-astro-data-formats.readthedocs.org/en/latest/irfs/psf/psf_table/index.html

For this the OGIP document recommends RAD, THETA, ENERGY order for better access performance:
http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_020/cal_gen_92_020.html

@jknodlseder - What's your preference here?

(to me, your argument that codes can re-order if they think it matters for performance makes sense, so I guess we should put ENERGY first there also?)

[jknodlseder]

I would try to follow OGIP as much as possible, so I would keep the established order.

Interesting to see that they were concerned by the access performance. They probably had specific applications in mind when writing the spec.

[cdeil]

OK, will keep RAD, THETA, ENERGY then.
I think their main argument is that the common case is to extract a PSF curve that depends on RAD, so that should be the first access for good memory / cache line performance.

In Gammapy we only compute a PSF kernel once at the start of the likelihood fit, so such concerns on IRF performance is irrelevant. It might matter for Gammalib, right?

[cdeil]

For EDISP2D is there any preference?

• ENERGY, THETA, MIGRA
• ENERGY, MIGRA, THETA

[jknodlseder]

Also here, energy axes first (ENERGY, MIGRA), offset angle follows (THETA).

[cdeil]

The recommended axis orders should now (a30965f ) be correct.

I've added this general statement:

The IRF format specs mention a recommended axis format and axis units.
But tools should not depend on this and instead:

• Use the axis order specified by the CREF header keyword (see :ref:fits-arrays-bintable-hdu)
• Use the axis unit specifiec by the CUNIT header keywords (see :ref:fits-arrays-bintable-hdu)

@jknodlseder - OK for the wording on CUNIT? You plan to change Gammalib to respect those units, right?

One last question: at the moment we have this general statement:

• For energy grids, see here <http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_003/cal_gen_92_003.html#tth_sEc7>__
  for basic recommendations. Column names should be ENERGY or ENERG_LO, ENERG_HI
  because that is used (consistently I think) for OGIP and Fermi-LAT.
  For separate HDUs, the extension names should be ENERGIES or EBOUNDS (used by Fermi-LAT consistently).

@jknodlseder - Does Gammalib currently support ENERGY instead of ENERG_LO, ENERG_HI if it's declared via CREF? Do you plan to support it?
I think we discussed this and agreed to only support ENERG_LO, ENERG_HI, right?

[cdeil]

Merging this now.

Further comments and suggestions welcome.
Feel free to make pull requests, or I can make edits in master or further pull requests.

FYI - For our Python prototype I'm still looking for feedback how to implement this well (abstract axis and array class, little code) based on astropy.io.fits: https://groups.google.com/forum/#!topic/astropy-dev/khqTG3741uk