Effective Area 3D

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,99 @@ def effective_area_per_energy_and_fov(
     )
 
     return effective_area(hist_selected, hist_simulated, area)
+
+
+def effective_area_3d_polar(
+    selected_events,
+    simulation_info,
+    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(
+        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=(
+            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, energy_bins, fov_longitude_bins, fov_latitude_bins
+):
+    """
+    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(
+        energy_bins, fov_longitude_bins, fov_latitude_bins
+    )
+
+    hist_selected, _ = np.histogramdd(
+        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=(
+            energy_bins.to_value(u.TeV),
+            fov_longitude_bins.to_value(u.deg),
+            fov_latitude_bins.to_value(u.deg),
+        ),
+    )
+
+    return effective_area(hist_selected, hist_simulated, area)

pyirf/simulations.py

@@ -1,5 +1,6 @@
 import astropy.units as u
 import numpy as np
+from .utils import cone_solid_angle
 
 __all__ = [
     'SimulatedEventsInfo',
@@ -168,6 +169,115 @@ def calculate_n_showers_per_energy_and_fov(self, energy_bins, fov_bins):
         fov_integral = self.calculate_n_showers_per_fov(fov_bins)
         return e_integral[:, np.newaxis] * fov_integral / self.n_showers
 
+    @u.quantity_input(
+        energy_bins=u.TeV, fov_offset_bins=u.deg, fov_position_angle_bins=u.rad
+    )
+    def calculate_n_showers_3d_polar(
+        self, energy_bins, fov_offset_bins, fov_position_angle_bins
+    ):
+        """
+        Calculate number of showers that were simulated in the given
+        energy and 2D fov bins in polar coordinates.
+
+        This assumes the events were generated uniformly distributed per solid angle,
+        and from a powerlaw in energy like CORSIKA simulates events.
+
+        Parameters
+        ----------
+        energy_bins: astropy.units.Quantity[energy]
+            The energy bin edges for which to calculate the number of simulated showers
+        fov_offset_bins: astropy.units.Quantity[angle]
+            The FOV radial bin edges for which to calculate the number of simulated showers
+        fov_position_angle_bins: astropy.units.Quantity[radian]
+            The FOV azimuthal bin edges for which to calculate the number of simulated showers
+
+        Returns
+        -------
+        n_showers: numpy.ndarray(ndim=3)
+            The expected number of events inside each of the
+            ``energy_bins``, ``fov_offset_bins`` and ``fov_position_angle_bins``.
+            Dimension (n_energy_bins, n_fov_offset_bins, n_fov_position_angle_bins)
+            This is a floating point number.
+            The actual numbers will follow a poissionian distribution around this
+            expected value.
+        """
+        e_fov_offset_integral = self.calculate_n_showers_per_energy_and_fov(
+            energy_bins, fov_offset_bins
+        )
+        position_angle_integral = np.full(
+            len(fov_position_angle_bins) - 1,
+            self.calculate_n_showers_per_fov(
+                np.linspace(self.viewcone_min, self.viewcone_max, 2)
+            )
+            / (len(fov_position_angle_bins) - 1),
+        )
+
+        return (
+            e_fov_offset_integral[:, :, np.newaxis]
+            * position_angle_integral
+            / self.n_showers
+        )
+
+    @u.quantity_input(
+        energy_bins=u.TeV, fov_longitude_bins=u.deg, fov_latitude_bins=u.rad
+    )
+    def calculate_n_showers_3d_lonlat(
+        self, energy_bins, fov_longitude_bins, fov_latitude_bins
+    ):
+        """
+        Calculate number of showers that were simulated in the given
+        energy and 2D fov bins in nominal coordinates.
+
+        This assumes the events were generated uniformly distributed per solid angle,
+        and from a powerlaw in energy like CORSIKA simulates events.
+
+        Parameters
+        ----------
+        energy_bins: astropy.units.Quantity[energy]
+            The energy bin edges for which to calculate the number of simulated showers
+        fov_longitude_bins: astropy.units.Quantity[angle]
+            The FOV longitude bin edges for which to calculate the number of simulated showers
+        fov_latitude_bins: astropy.units.Quantity[angle]
+            The FOV latitude bin edges for which to calculate the number of simulated showers
+
+        Returns
+        -------
+        n_showers: numpy.ndarray(ndim=3)
+            The expected number of events inside each of the
+            ``energy_bins``, ``fov_longitude_bins`` and ``fov_latitude_bins``.
+            Dimension (n_energy_bins, n_fov_longitude_bins, n_fov_latitude_bins)
+            This is a floating point number.
+            The actual numbers will follow a poissionian distribution around this
+            expected value.
+        """
+        fov_bin_midpoints_lon, fov_bin_midpoints_lat = np.meshgrid(
+            (fov_longitude_bins[:-1] + fov_longitude_bins[1:]) / 2,
+            (fov_latitude_bins[:-1] + fov_latitude_bins[1:]) / 2,
+        )
+        mask_outside_fov = np.logical_or(
+            np.sqrt(fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2)
+            > self.viewcone_max,
+            np.sqrt(fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2)
+            < self.viewcone_min,
+        )
+        A_bins = (
+            np.outer(
+                np.diff(np.sin(fov_latitude_bins.to_value(u.rad))),
+                np.diff(fov_longitude_bins.to_value(u.rad)),
+            )
+            * u.sr
+        )
+        shower_density = self.n_showers / (
+            cone_solid_angle(self.viewcone_max) - cone_solid_angle(self.viewcone_min)
+        )
+
+        fov_integral = shower_density * A_bins
+        e_integral = self.calculate_n_showers_per_energy(energy_bins)
+
+        fov_integral[mask_outside_fov] = 0
+
+        return (e_integral[:, np.newaxis, np.newaxis] * fov_integral) / self.n_showers
+
     def __repr__(self):
         return (
             f"{self.__class__.__name__}("

pyirf/utils.py

@@ -1,6 +1,7 @@
 import numpy as np
 import astropy.units as u
 from astropy.coordinates import angular_separation
+from astropy.coordinates import position_angle
 
 from .exceptions import MissingColumns, WrongColumnUnit
 
@@ -9,6 +10,7 @@
     "is_scalar",
     "calculate_theta",
     "calculate_source_fov_offset",
+    "calculate_source_fov_position_angle",
     "check_histograms",
     "cone_solid_angle",
 ]
@@ -88,6 +90,34 @@ def calculate_source_fov_offset(events, prefix="true"):
     return theta.to(u.deg)
 
 
+def calculate_source_fov_position_angle(events, prefix="true"):
+    """Calculate position_angle of true positions relative to pointing positions.
+
+    Parameters
+    ----------
+    events : astropy.QTable
+        Astropy Table object containing the reconstructed events information.
+
+    prefix: str
+        Column prefix for az / alt, can be used to calculate reco or true
+        source fov position_angle.
+
+    Returns
+    -------
+    phi: astropy.units.Quantity
+        Position angle of the true positions relative to the pointing positions
+        in the sky.
+    """
+    phi = position_angle(
+        events["pointing_az"],
+        events["pointing_alt"],
+        events[f"{prefix}_az"],
+        events[f"{prefix}_alt"],
+    )
+
+    return phi.to(u.deg)
+
+
 def check_histograms(hist1, hist2, key="reco_energy"):
     """
     Check if two histogram tables have the same binning