Improving OGIP compatibility of EDispKernel.from_hdulist()

[k-egg]

Hello, I was asked to open an issue to discuss this topic further. Some discussion on this has already taken place in the gammapyXray and gammapy dev meeting.

The EDispKernel.from_hdulist() (also called by EDispKernel.read()) function was created to read OGIP format RMF files, the current version does, however, not accommodate all OGIP compliant formats. OGIP format specifications are defined in this memo.
Fermi LAT DRM files face additional issues, these go beyond the scope of the OGIP format.

To ensure that all OGIP compliant RMF files can be seamlessly read by EDispKernel.from_hdulist() two things need to be changed in the code:

1. Use the MATRIX extension TLMIN4 keyword
This keyword, taking a value of 0 or 1, signifies whether the indexing of the matrix values conducted with the F_CHAN column (usually the 4th column in the table) starts at 0 or 1. This is necessary to ensure proper indexing, the matrix might otherwise be shifted by one energy bin or, as in the case of the eROSITA RMF, fail to load at all.

2. Reinstate scalar format option for N_CHAN and F_CHAN keywords
The version of this function in Luca's and Mireia's repository includes the scalar option, the current 1.2 release, however does not. If N_CHAN and F_CHAN are given as scalars (e.g. by XMM MOS1/2), the function tries to index them and fails.

Suggested changes:
Try reading TLMIN* keyword, use value as index offset for indexing the matrix elements, should no TLMIN* keyword be present default to 0. Reinstate distinction between scalar and array N_CHAN/F_CHAN entries.

This version of the function can be found on github here.

TLMIN is found by:

         #check for TLMIN keyword to determine if indexing starts at 0 or 1:
        ind_offset=0
        try:
            ind_offset=int(matrix_hdu.header["TLMIN*"][0])

Enumerating over the data is then implemented this way:

        for i, l in enumerate(data):
            if l.field("N_GRP"):
                m_start = 0
                for k in range(l.field("N_GRP")):
                    if np.isscalar(l.field("N_CHAN")):
                        chan_min = l.field("F_CHAN")-ind_offset
                        chan_max = l.field("F_CHAN") + l.field("N_CHAN")-ind_offset
                        n_chan=l.field("N_CHAN")
                    else:
                        chan_min = l.field("F_CHAN")[k]-ind_offset
                        chan_max = l.field("F_CHAN")[k] + l.field("N_CHAN")[k]-ind_offset
                        n_chan=l.field("N_CHAN")[k]

                    pdf_matrix[i, chan_min:chan_max] = l.field("MATRIX")[
                        m_start : m_start + n_chan  # noqa: E203
                    ]
                    m_start += n_chan

I think it would be great to include something like this in gammapy to make it more compatible with OGIP data and consequently all kinds of X-ray missions. Please let me know what you think, I am very open to feedback!

I believe Mireia and Fabio also have further notes on this, as they have implemented, e.g., how to avoid energy bins with a lower threshold of zero for XMM RMFs, and more. Maybe we can try to find a good combination of all these aspects with maximum compatibility for OGIP data!

Let's discuss how this should be best implemented in gammapy!