diff --git a/magicctapipe/io/gadf.py b/magicctapipe/io/gadf.py index 98a12773..2871bcb4 100644 --- a/magicctapipe/io/gadf.py +++ b/magicctapipe/io/gadf.py @@ -31,7 +31,9 @@ @u.quantity_input(reco_energy_bins=u.TeV, fov_offset_bins=u.deg) -def create_gh_cuts_hdu(gh_cuts, reco_energy_bins, fov_offset_bins, **header_cards): +def create_gh_cuts_hdu( + gh_cuts, reco_energy_bins, fov_offset_bins, point_like, **header_cards +): """ Creates a fits binary table HDU for dynamic gammaness cuts. @@ -44,6 +46,8 @@ def create_gh_cuts_hdu(gh_cuts, reco_energy_bins, fov_offset_bins, **header_card Bin edges in the reconstructed energy fov_offset_bins : astropy.units.quantity.Quantity Bin edges in the field of view offset + point_like : bool + Meaning: true = IRFs are point-like, false = IRFS are full-enclosure **header_cards : dict Additional metadata to add to the header @@ -73,7 +77,7 @@ def create_gh_cuts_hdu(gh_cuts, reco_energy_bins, fov_offset_bins, **header_card ("CREATOR", f"magicctapipe v{MCP_VERSION}"), ("HDUCLAS1", "RESPONSE"), ("HDUCLAS2", "GH_CUTS"), - ("HDUCLAS3", "POINT-LIKE"), + ("HDUCLAS3", "POINT-LIKE" if point_like else "FULL-ENCLOSURE"), ("HDUCLAS4", "GH_CUTS_2D"), ("DATE", Time.now().utc.iso), ] diff --git a/magicctapipe/io/tests/test_gadf.py b/magicctapipe/io/tests/test_gadf.py index c5b31edc..559a1c30 100644 --- a/magicctapipe/io/tests/test_gadf.py +++ b/magicctapipe/io/tests/test_gadf.py @@ -60,7 +60,9 @@ def event(): def test_create_gh_cuts_hdu(gammaness_cut, energy_bins, fov_bins, header): - g_cuts_fits = create_gh_cuts_hdu(gammaness_cut, energy_bins, fov_bins, **header) + g_cuts_fits = create_gh_cuts_hdu( + gammaness_cut, energy_bins, fov_bins, True, **header + ) g_cuts = g_cuts_fits.data assert np.allclose(g_cuts["ENERG_LO"][0], np.array([0.1, 1.0, 10.0])) assert np.allclose(g_cuts["ENERG_HI"][0], np.array([1.0, 10.0, 100.0])) diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py index e2727cee..9ea97784 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py @@ -366,10 +366,12 @@ def create_irf( # Calculate dynamic gammaness cuts gh_percentile = 100 * (1 - gh_efficiency) - + withinfov = ( + event_table_gamma["true_source_fov_offset"].to("deg") < fov_offset_bins[-1] + ) * (event_table_gamma["true_source_fov_offset"].to("deg") > fov_offset_bins[0]) cut_table_gh = calculate_percentile_cut( - values=event_table_gamma["gammaness"], - bin_values=event_table_gamma["reco_energy"], + values=event_table_gamma[withinfov]["gammaness"], + bin_values=event_table_gamma[withinfov]["reco_energy"], bins=energy_bins, fill_value=gh_cut_min, percentile=gh_percentile, @@ -423,6 +425,7 @@ def create_irf( gh_cuts=gh_cuts, reco_energy_bins=energy_bins, fov_offset_bins=fov_offset_bins, + point_like=is_point_like, **extra_header, ) @@ -472,9 +475,16 @@ def create_irf( # Calculate dynamic theta cuts theta_percentile = 100 * theta_efficiency + withinfov = ( + event_table_gamma["true_source_fov_offset"].to("deg") + < fov_offset_bins[-1] + ) * ( + event_table_gamma["true_source_fov_offset"].to("deg") + > fov_offset_bins[0] + ) cut_table_theta = calculate_percentile_cut( - values=event_table_gamma["theta"], - bin_values=event_table_gamma["reco_energy"], + values=event_table_gamma[withinfov]["theta"], + bin_values=event_table_gamma[withinfov]["reco_energy"], bins=energy_bins, fill_value=theta_cut_max, percentile=theta_percentile, diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py b/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py index 097f1443..7255ec34 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py @@ -206,11 +206,12 @@ def dl2_to_dl3(input_file_dl2, input_dir_irf, output_dir, config): aeff_interp = irf_data["effective_area"][0] print("skipping interpolation since only one point is given") + point_like = irf_data["effective_area"].shape[-1] == 1 aeff_hdu = create_aeff2d_hdu( effective_area=aeff_interp, true_energy_bins=irf_data["energy_bins"], fov_offset_bins=irf_data["fov_offset_bins"], - point_like=True, + point_like=point_like, extname="EFFECTIVE AREA", **extra_header, ) @@ -247,7 +248,7 @@ def dl2_to_dl3(input_file_dl2, input_dir_irf, output_dir, config): true_energy_bins=irf_data["energy_bins"], migration_bins=irf_data["migration_bins"], fov_offset_bins=irf_data["fov_offset_bins"], - point_like=True, + point_like=point_like, extname="ENERGY DISPERSION", ) @@ -328,6 +329,7 @@ def dl2_to_dl3(input_file_dl2, input_dir_irf, output_dir, config): gh_cuts=gh_cuts_interp, reco_energy_bins=irf_data["energy_bins"], fov_offset_bins=irf_data["fov_offset_bins"], + point_like=point_like, **extra_header, )