Merge pull request #1270 from cta-observatory/fix_interp_aeff

Fix the shape of the interpolated Effective Area

lstchain/high_level/interpolate.py

@@ -422,7 +422,7 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
     main_headers = fits.open(irfs[0])[1].header
 
     point_like = main_headers["HDUCLAS3"] == "POINT-LIKE"
-    
+
     # Update headers to be added to the final IRFs
     extra_headers = dict()
 
@@ -471,19 +471,19 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
     }
 
     for hdu_name in hdu_names:
-        
+
         if hdu_name == "BACKGROUND":
             log.warning("The interpolation of BACKGROUND is not yet supported")
             continue
-            
+
         if interp_col_keys[hdu_name] == 'cut':
             gadf_irf = False
         else:
-            gadf_irf = True    
-            
+            gadf_irf = True
+
         irf_list = load_irf_grid(
-            irfs, 
-            extname=hdu_name, 
+            irfs,
+            extname=hdu_name,
             interp_col=interp_col_keys[hdu_name],
             gadf_irf=gadf_irf
         )
@@ -492,17 +492,17 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
         if hdu_name in ['EFFECTIVE AREA', 'ENERGY DISPERSION', 'PSF']:
             e_true = np.append(temp_irf["ENERG_LO"][0], temp_irf["ENERG_HI"][0][-1])
             fov_off = np.append(temp_irf["THETA_LO"][0], temp_irf["THETA_HI"][0][-1])
-            
+
             if hdu_name == "EFFECTIVE AREA":
                 aeff_estimator = EffectiveAreaEstimator(
                     grid_points=irf_pars_sel,
                     effective_area=irf_list,
                     interpolator_kwargs={"method": interp_method},
                 )
                 aeff_interp = aeff_estimator(interp_pars_sel)
-                
+
                 aeff_hdu_interp = create_aeff2d_hdu(
-                    effective_area=aeff_interp.T[0],
+                    effective_area=aeff_interp[0],
                     true_energy_bins=e_true,
                     fov_offset_bins=fov_off,
                     point_like=point_like,
@@ -519,7 +519,7 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
                     energy_dispersion=irf_list,
                 )
                 edisp_interp = edisp_estimator(interp_pars_sel)
-                
+
                 edisp_hdu_interp = create_energy_dispersion_hdu(
                     energy_dispersion=edisp_interp[0],
                     true_energy_bins=e_true,
@@ -556,7 +556,7 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
             cut_header = fits.Header()
             cut_header["CREATOR"] = f"lstchain v{__version__}"
             cut_header["DATE"] = Time.now().utc.iso
-            
+
             for k, v in extra_headers.items():
                 cut_header[k] = v
 

lstchain/high_level/tests/test_interpolate.py

@@ -257,17 +257,17 @@ def test_interp_irf(
                 psf_1_meta = fits.open(irf)["PSF"].header
                 psf_2 = psf_1.copy()
                 psf_2_meta = psf_1_meta.copy()
-                
+
                 psf_1["RPSF"][0] *= factor_czd
                 psf_2["RPSF"][0] *= factor_czd * factor_sd
-                
+
                 psf_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
                 psf_1_meta["AZ_PNT"] = (az_1 * 180 / np.pi, "deg")
                 psf_1_meta["B_DELTA"] = (del_1 * 180 / np.pi, "deg")
                 psf_2_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
                 psf_2_meta["AZ_PNT"] = (az_2 * 180 / np.pi, "deg")
                 psf_2_meta["B_DELTA"] = (del_2 * 180 / np.pi, "deg")
-                
+
                 psf_hdu_1 = fits.BinTableHDU(
                     psf_1, header=psf_1_meta, name="PSF"
                 )
@@ -333,10 +333,14 @@ def test_interp_irf(
     hdu_en_srcdep = interpolate_irf(irfs_en_srcdep, data_pars)
     hdu_en_srcdep.writeto(irf_file_en_srcdep_final, overwrite=True)
 
+    aeff_shape_final = Table.read(irf_file_g_final, hdu=1)["EFFAREA"].shape
+    aeff_shape_2 = Table.read(irf_file_g_2, hdu=1)["EFFAREA"].shape
+
     assert hdu_g[1].header["ZEN_PNT"] == zen_t
     assert irf_file_g_2.exists()
     assert irf_file_g_3.exists()
     assert irf_file_g_final.exists()
+    assert aeff_shape_final == aeff_shape_2
 
     assert hdu_en[1].header["ZEN_PNT"] == zen_t
     assert irf_file_en_2.exists()