Effective Area 3D

docs/changes/281.feature.rst

@@ -0,0 +1,2 @@
+
+Add 3D effective area functions for lon/lat and theta/phi coordinates and some necessary utiliy functions.

pyirf/coordinates.py

@@ -0,0 +1,72 @@
+import astropy.units as u
+from astropy.coordinates import SkyCoord, SkyOffsetFrame, angular_separation, position_angle
+
+__all__ = [
+    'gadf_fov_coords_lon_lat',
+    'gadf_fov_coords_theta_phi',
+]
+
+
+def gadf_fov_coords_lon_lat(lon, lat, pointing_lon, pointing_lat):
+    """Transform sky coordinates to field-of-view longitude-latitude coordinates in accordance with
+    the definitions laid out by the Gamma Astro Data Format.
+
+    GADF documentation here:
+    https://gamma-astro-data-formats.readthedocs.io/en/latest/general/coordinates.html
+
+    Parameters
+    ----------
+    lon, lat : `~astropy.units.Quantity`
+        Sky coordinate to be transformed.
+    pointing_lon, pointing_lat : `~astropy.units.Quantity`
+        Coordinate specifying the pointing position.
+        (i.e. the center of the field of view.)
+
+    Returns
+    -------
+    lon_t, lat_t : `~astropy.units.Quantity`
+        Transformed field-of-view coordinate.
+    """
+    # Create a frame that is centered on the pointing position
+    center = SkyCoord(pointing_lon, pointing_lat)
+    fov_frame = SkyOffsetFrame(origin=center)
+
+    # Define coordinate to be transformed.
+    target_sky = SkyCoord(lon, lat)
+
+    # Transform into FoV-system
+    target_fov = target_sky.transform_to(fov_frame)
+
+    # Switch sign of longitude angle since this axis is
+    # reversed in our definition of the FoV-system
+    return -target_fov.lon, target_fov.lat
+
+
+def gadf_fov_coords_theta_phi(lon, lat, pointing_lon, pointing_lat):
+    """Transform sky coordinates to field-of-view theta-phi coordinates in accordance with
+    the definitions laid out by the Gamma Astro Data Format.
+
+    GADF documentation here:
+    https://gamma-astro-data-formats.readthedocs.io/en/latest/general/coordinates.html
+
+    Parameters
+    ----------
+    lon, lat : `~astropy.units.Quantity`
+        Sky coordinate to be transformed.
+    pointing_lon, pointing_lat : `~astropy.units.Quantity`
+        Coordinate specifying the pointing position.
+        (i.e. the center of the field of view.)
+
+    Returns
+    -------
+    theta, phi : `~astropy.units.Quantity`
+        Transformed field-of-view coordinate.
+    """
+
+    theta = angular_separation(pointing_lon, pointing_lat, lon, lat)
+
+    # astropy defines the position angle as increasing towards east of north
+    phi = position_angle(pointing_lon, pointing_lat, lon, lat)
+
+    # GADF defines FOV PHI opposite to the position angle way so the sign is switched
+    return theta.to(u.deg), (-phi).wrap_at(360 * u.deg).to(u.deg)

pyirf/irf/__init__.py

@@ -2,6 +2,8 @@
     effective_area,
     effective_area_per_energy_and_fov,
     effective_area_per_energy,
+    effective_area_3d_polar,
+    effective_area_3d_lonlat,
 )
 from .energy_dispersion import energy_dispersion
 from .psf import psf_table
@@ -11,6 +13,8 @@
     "effective_area",
     "effective_area_per_energy",
     "effective_area_per_energy_and_fov",
+    "effective_area_3d_polar",
+    "effective_area_3d_lonlat",
     "energy_dispersion",
     "psf_table",
     "background_2d",

pyirf/irf/effective_area.py

@@ -7,6 +7,8 @@
     "effective_area",
     "effective_area_per_energy",
     "effective_area_per_energy_and_fov",
+    "effective_area_3d_polar",
+    "effective_area_3d_lonlat",
 ]
 
 
@@ -85,3 +87,104 @@ def effective_area_per_energy_and_fov(
     )
 
     return effective_area(hist_selected, hist_simulated, area)
+
+
+def effective_area_3d_polar(
+    selected_events,
+    simulation_info,
+    true_energy_bins,
+    fov_offset_bins,
+    fov_position_angle_bins,
+):
+    """
+    Calculate effective area in bins of true energy, field of view offset, and field of view position angle.
+
+    Parameters
+    ----------
+    selected_events: astropy.table.QTable
+        DL2 events table, required columns for this function:
+        - `true_energy`
+        - `true_source_fov_offset`
+        - `true_source_fov_position_angle`
+    simulation_info: pyirf.simulations.SimulatedEventsInfo
+        The overall statistics of the simulated events
+    true_energy_bins: astropy.units.Quantity[energy]
+        The true energy bin edges in which to calculate effective area.
+    fov_offset_bins: astropy.units.Quantity[angle]
+        The field of view radial bin edges in which to calculate effective area.
+    fov_position_angle_bins: astropy.units.Quantity[radian]
+        The field of view azimuthal bin edges in which to calculate effective area.
+    """
+    area = np.pi * simulation_info.max_impact**2
+
+    hist_simulated = simulation_info.calculate_n_showers_3d_polar(
+        true_energy_bins, fov_offset_bins, fov_position_angle_bins
+    )
+
+    hist_selected, _ = np.histogramdd(
+        np.column_stack(
+            [
+                selected_events["true_energy"].to_value(u.TeV),
+                selected_events["true_source_fov_offset"].to_value(u.deg),
+                selected_events["true_source_fov_position_angle"].to_value(u.rad),
+            ]
+        ),
+        bins=(
+            true_energy_bins.to_value(u.TeV),
+            fov_offset_bins.to_value(u.deg),
+            fov_position_angle_bins.to_value(u.rad),
+        ),
+    )
+
+    return effective_area(hist_selected, hist_simulated, area)
+
+
+def effective_area_3d_lonlat(
+    selected_events,
+    simulation_info,
+    true_energy_bins,
+    fov_longitude_bins,
+    fov_latitude_bins,
+    subpixels=20,
+):
+    """
+    Calculate effective area in bins of true energy, field of view longitude, and field of view latitude.
+
+    Parameters
+    ----------
+    selected_events: astropy.table.QTable
+        DL2 events table, required columns for this function:
+        - `true_energy`
+        - `true_source_fov_lon`
+        - `true_source_fov_lat`
+    simulation_info: pyirf.simulations.SimulatedEventsInfo
+        The overall statistics of the simulated events
+    true_energy_bins: astropy.units.Quantity[energy]
+        The true energy bin edges in which to calculate effective area.
+    fov_longitude_bins: astropy.units.Quantity[angle]
+        The field of view longitude bin edges in which to calculate effective area.
+    fov_latitude_bins: astropy.units.Quantity[angle]
+        The field of view latitude bin edges in which to calculate effective area.
+    """
+    area = np.pi * simulation_info.max_impact**2
+
+    hist_simulated = simulation_info.calculate_n_showers_3d_lonlat(
+        true_energy_bins, fov_longitude_bins, fov_latitude_bins, subpixels=subpixels
+    )
+
+    selected_columns = np.column_stack(
+            [
+                selected_events["true_energy"].to_value(u.TeV),
+                selected_events["true_source_fov_lon"].to_value(u.deg),
+                selected_events["true_source_fov_lat"].to_value(u.deg),
+            ]
+        )
+    bins = (
+            true_energy_bins.to_value(u.TeV),
+            fov_longitude_bins.to_value(u.deg),
+            fov_latitude_bins.to_value(u.deg),
+    )
+
+    hist_selected, _ = np.histogramdd(selected_columns, bins=bins)
+
+    return effective_area(hist_selected, hist_simulated, area)

pyirf/irf/tests/test_effective_area.py

@@ -93,3 +93,148 @@ def test_effective_area_energy_fov():
     assert area.unit == u.m ** 2
     assert u.allclose(area[:, 0], [200, 20] * u.m ** 2)
     assert u.allclose(area[:, 1], [100, 10] * u.m ** 2)
+
+
+def test_effective_area_3d_polar():
+    from pyirf.irf import effective_area_3d_polar
+    from pyirf.simulations import SimulatedEventsInfo
+
+    true_energy_bins = [0.1, 1.0, 10.0] * u.TeV
+    # choose edges so that half are in each bin in fov
+    fov_offset_bins = [0, np.arccos(0.98), np.arccos(0.96)] * u.rad
+    # choose edges so that they are equal size and cover the whole circle
+    fov_pa_bins = [0, np.pi, 2*np.pi] * u.rad
+    center_1_o, center_2_o = 0.5 * (fov_offset_bins[:-1] + fov_offset_bins[1:]).to_value(
+        u.deg
+    )
+    center_1_pa, center_2_pa = 0.5 * (fov_pa_bins[:-1] + fov_pa_bins[1:]).to_value(
+        u.deg
+    )
+
+    selected_events = QTable(
+        {
+            "true_energy": np.concatenate(
+                [
+                    np.full(1000, 0.5),
+                    np.full(10, 5),
+                    np.full(500, 0.5),
+                    np.full(5, 5),
+                    np.full(1000, 0.5),
+                    np.full(10, 5),
+                    np.full(500, 0.5),
+                    np.full(5, 5),
+                ]
+            )
+            * u.TeV,
+            "true_source_fov_offset": np.concatenate(
+                [
+                    np.full(1010, center_1_o),
+                    np.full(505, center_2_o),
+                    np.full(1010, center_1_o),
+                    np.full(505, center_2_o),
+                ]
+            )
+            * u.deg,
+            "true_source_fov_position_angle": np.append(
+                np.full(1515, center_1_pa), np.full(1515, center_2_pa)
+            )
+            * u.deg,
+        }
+    )
+
+    # this should give 100000 events in the first bin and 10000 in the second
+    simulation_info = SimulatedEventsInfo(
+        n_showers=110000,
+        energy_min=true_energy_bins[0],
+        energy_max=true_energy_bins[-1],
+        max_impact=100 / np.sqrt(np.pi) * u.m,  # this should give a nice round area
+        spectral_index=-2,
+        viewcone_min=0 * u.deg,
+        viewcone_max=fov_offset_bins[-1],
+    )
+
+    area = effective_area_3d_polar(
+        selected_events, simulation_info, true_energy_bins, fov_offset_bins, fov_pa_bins
+    )
+
+    assert area.shape == (
+        len(true_energy_bins) - 1, len(fov_offset_bins) - 1, len(fov_pa_bins) - 1
+    )
+    assert area.unit == u.m ** 2
+    assert u.allclose(area[:, 0, :], [[400, 400],[40,40]] * u.m ** 2)
+    assert u.allclose(area[:, 1, :], [[200, 200],[20,20]] * u.m ** 2)
+
+
+def test_effective_area_3d_lonlat():
+    from pyirf.irf import effective_area_3d_lonlat
+    from pyirf.simulations import SimulatedEventsInfo
+
+    true_energy_bins = [0.1, 1.0, 10.0] * u.TeV
+    # choose edges so that a quarter are in each bin in fov
+    fov_lon_bins = [-1.0, 0, 1.0] * u.deg
+    fov_lat_bins = [-1.0, 0, 1.0] * u.deg
+    center_1_lon, center_2_lon = 0.5 * (fov_lon_bins[:-1] + fov_lon_bins[1:]).to_value(
+        u.deg
+    )
+    center_1_lat, center_2_lat = 0.5 * (fov_lat_bins[:-1] + fov_lat_bins[1:]).to_value(
+        u.deg
+    )
+
+    selected_events = QTable(
+        {
+            "true_energy": np.concatenate(
+                [
+                    np.full(1000, 0.5),
+                    np.full(10, 5),
+                    np.full(500, 0.5),
+                    np.full(5, 5),
+                    np.full(1000, 0.5),
+                    np.full(10, 5),
+                    np.full(500, 0.5),
+                    np.full(5, 5),
+                ]
+            )
+            * u.TeV,
+            "true_source_fov_lon": np.concatenate(
+                [
+                np.full(1010, center_1_lon),
+                np.full(505, center_2_lon),
+                np.full(1010, center_1_lon),
+                np.full(505, center_2_lon),
+                ]
+            )
+            * u.deg,
+            "true_source_fov_lat": np.append(
+                np.full(1515, center_1_lat), np.full(1515, center_2_lat)
+            )
+            * u.deg,
+        }
+    )
+
+    # this should give 100000 events in the first bin and 10000 in the second
+    simulation_info = SimulatedEventsInfo(
+        n_showers=110000,
+        energy_min=true_energy_bins[0],
+        energy_max=true_energy_bins[-1],
+        max_impact=100 / np.sqrt(np.pi) * u.m,  # this should give a nice round area
+        spectral_index=-2,
+        viewcone_min=0 * u.deg,
+        viewcone_max=fov_lon_bins[-1],
+    )
+
+    area = effective_area_3d_lonlat(
+        selected_events,
+        simulation_info,
+        true_energy_bins,
+        fov_lon_bins,
+        fov_lat_bins,
+        subpixels=20,
+    )
+
+    assert area.shape == (
+        len(true_energy_bins) - 1, len(fov_lon_bins) - 1, len(fov_lat_bins) - 1
+    )
+    assert area.unit == u.m ** 2
+    # due to inexact approximation of the circular FOV area in the lon-lat bins tolerance of 1%
+    assert u.allclose(area[:, 0, :], [[400, 400],[40,40]] * u.m ** 2,rtol=1e-2)
+    assert u.allclose(area[:, 1, :], [[200, 200],[20,20]] * u.m ** 2,rtol=1e-2)