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_nominal,
 )
 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_nominal",
     "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_nominal",
 ]
 
 
@@ -85,3 +87,69 @@
     )
 
     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.array([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)]).T,
+            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_nominal(
+    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_nominal(energy_bins, fov_longitude_bins, fov_latitude_bins)
+
+    hist_selected, _ = np.histogramdd(
+        np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_lon'].value, selected_events['true_source_fov_lat'].value]).T,
+        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,81 @@
         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.ones(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_nominal(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 @@
     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