interpolate_irfs reduces dimentions of EFFAREA IRF

[IevgenVovk]

lstchain.high_level.interpolate_irfs() reduces the dimensions of the returned effective area array, resulting in inconsistency between the specified energy-offset binning and array shape of the written IRF. To fix this, apparently, it is sufficient to replace

    aeff_hdu_interp = create_aeff2d_hdu(
            effective_area=aeff_interp.T[0],
            ...

with

    aeff_hdu_interp = create_aeff2d_hdu(
            effective_area=aeff_interp.T,
            ...

in L521 of "interpolate.py".

The issue manifests itself in v0.10.11 and is not fixed in master so far.

[chaimain]

Can you show an example of the difference? I cannot reproduce this difference and get the same array shape of the EFFAREA column - (1, 1, 25) for IRFs and DL3 files produced using the current main branch of lstchain.

[IevgenVovk]

Creation of a DL3 file:

lstchain_create_dl3_file \
            --input-dl2 /fefs/aswg/workspace/ievgen.vovk/lst/crab/data/dl2/dl2_LST-1.Run15864.h5 \
            --output-dl3-path "/fefs/aswg/workspace/ievgen.vovk/lst/crab/data/dl3/" \
            --input-irf-path "/fefs/aswg/workspace/ievgen.vovk/lst/allsky-mc-irf/irf/all/" \
            --irf-file-pattern "*irf*.fits" \
            --source-name "crab" \
            --config "/fefs/aswg/workspace/ievgen.vovk/lst/allsky-mc-irf/irf/irf-config.json" \
            --overwrite \
            --gzip

Aeff array dimensions check:

from astropy.io import fits

hdus = fits.open('/fefs/aswg/workspace/ievgen.vovk/lst/crab/data/dl3/dl3_LST-1.Run15864.fits.gz')
hdus[4].data

->

FITS_rec([([5.00000000e-03, 1.58113883e-02, 5.00000000e-02, 1.58113883e-01, 5.00000000e-01, 1.58113883e+00, 5.00000000e+00, 1.58113883e+01, 5.00000000e+01, 1.58113883e+02], [1.58113883e-02, 5.00000000e-02, 1.58113883e-01, 5.00000000e-01, 1.58113883e+00, 5.00000000e+00, 1.58113883e+01, 5.00000000e+01, 1.58113883e+02, 5.00000000e+02], [0. , 0.5, 1. , 1.5], [0.5, 1. , 1.5, 2. ], [[7.46435097e+00, 1.98198758e+03, 4.70403938e+04, 1.89589412e+05, 2.56331808e+05, 2.11730543e+05, 1.58451248e+05, 9.11920427e+04, 2.70711109e+02, 1.00000000e+00]])],
         dtype=(numpy.record, [('ENERG_LO', '>f8', (10,)), ('ENERG_HI', '>f8', (10,)), ('THETA_LO', '>f8', (4,)), ('THETA_HI', '>f8', (4,)), ('EFFAREA', '>f8', (1, 10))]))

The resulting dimensions are (1, 10) instead of (4, 10)

[IevgenVovk]

Checking this more closely, the issue might be in the order of original dimensions of the interpolated collection area. I've added a few lines before L521 like this:

        aeff_interp = aeff_estimator(interp_pars_sel)
        print(aeff_interp.shape)
        print(aeff_interp.T.shape)

this results in

(1, 10, 4)
(4, 10, 1)

This way to keep the useful dimensions one would need to use effective_area=aeff_interp[0].T instead of effective_area=aeff_interp.T[0]. My original suggestion thus looks wrong, while the issue exists.

[chaimain]

Thanks for the example and suggestion! Indeed, I probably only tested with point-like IRFs, and hence made that error.

[IevgenVovk]

Trying to load the resulting files with gammapy, it further seems the transpose operation is not need (gammapy complains about not matching - effectively swapped - axes).

[IevgenVovk]

I.e. that line should read effective_area=aeff_interp[0]

[chaimain]

Indeed. I checked with the test data in lstchain as we have diffuse gamma MC as test data. I will make a quick PR on this fix.

[IevgenVovk]

Thank you, @chaimain