diff --git a/docs/changes/281.feature.rst b/docs/changes/281.feature.rst new file mode 100644 index 000000000..1ad9368c6 --- /dev/null +++ b/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. diff --git a/pyirf/coordinates.py b/pyirf/coordinates.py new file mode 100644 index 000000000..d8e2eca66 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/pyirf/irf/__init__.py b/pyirf/irf/__init__.py index 3a7224684..f6eb1608d 100644 --- a/pyirf/irf/__init__.py +++ b/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", diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index daa13e767..a204cd556 100644 --- a/pyirf/irf/effective_area.py +++ b/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) diff --git a/pyirf/irf/tests/test_effective_area.py b/pyirf/irf/tests/test_effective_area.py index 802fdf32d..abfdf2881 100644 --- a/pyirf/irf/tests/test_effective_area.py +++ b/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) \ No newline at end of file diff --git a/pyirf/simulations.py b/pyirf/simulations.py index 2bedf4c1b..d7f6c3be9 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -1,5 +1,9 @@ import astropy.units as u import numpy as np +from astropy.coordinates import angular_separation +from .utils import cone_solid_angle +from .utils import rectangle_solid_angle +from .binning import bin_center __all__ = [ 'SimulatedEventsInfo', @@ -168,6 +172,107 @@ 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 + ) + viewcone_integral = self.calculate_n_showers_per_fov( + [self.viewcone_min, self.viewcone_max] * u.deg + ) + + n_bins_pa = len(fov_position_angle_bins) - 1 + position_angle_integral = np.full(n_bins_pa, viewcone_integral / n_bins_pa) + + total_integral = e_fov_offset_integral[:, :, np.newaxis] * position_angle_integral + + return total_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, subpixels=20 + ): + """ + 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_overlap = _fov_lonlat_grid_overlap( + self, fov_longitude_bins, fov_latitude_bins, subpixels=subpixels + ) + + bin_grid_lon, bin_grid_lat = np.meshgrid(fov_longitude_bins,fov_latitude_bins) + bin_area = rectangle_solid_angle( + bin_grid_lon[:-1,:-1], + bin_grid_lon[1:,1:], + bin_grid_lat[:-1,:-1], + bin_grid_lat[1:,1:], + ) + viewcone_area = cone_solid_angle(self.viewcone_max) - cone_solid_angle(self.viewcone_min) + + shower_density = self.n_showers / viewcone_area + + fov_integral = shower_density * bin_area + e_integral = self.calculate_n_showers_per_energy(energy_bins) + + fov_integral = fov_overlap * fov_integral + + return (e_integral[:, np.newaxis, np.newaxis] * fov_integral) / self.n_showers + def __repr__(self): return ( f"{self.__class__.__name__}(" @@ -224,3 +329,116 @@ def _viewcone_pdf_integral(viewcone_min, viewcone_max, fov_low, fov_high): if scalar: return np.squeeze(integral) return integral + +def _fov_lonlat_grid_overlap(self, bin_edges_lon, bin_edges_lat, subpixels=20): + """ + When binning in longitude and latitude there might be bins that are fully or + partially outside of the simulated viewcone window. + + Subdivides each bin on the edge of the viewcone and calculates the fraction of + 'subpixels' whose centers are inside the viewcone, which approximates the area of the + bin inside the viewcone. + """ + # treat edge cases where all bins are either fully inside or outside of the viewcone + bin_grid_lon, bin_grid_lat = np.meshgrid(bin_edges_lon, bin_edges_lat) + + bin_dist = angular_separation( + bin_grid_lon, + bin_grid_lat, + 0, + 0, + ).to_value(u.deg) + + all_outside = np.logical_or( + np.all(bin_dist <= self.viewcone_min.to_value(u.deg)), + np.all(bin_dist >= self.viewcone_max.to_value(u.deg)), + ) + all_inside = np.logical_and( + np.all(bin_dist <= self.viewcone_max.to_value(u.deg)), + np.all(bin_dist >= self.viewcone_min.to_value(u.deg)), + ) + + if all_outside: + return 0 + elif all_inside: + return 1 + + + # define grid of bin centers + fov_bin_centers_lon = bin_center(bin_edges_lon) + fov_bin_centers_lat = bin_center(bin_edges_lat) + + bin_centers_grid_lon, bin_centers_grid_lat = np.meshgrid( + fov_bin_centers_lon, fov_bin_centers_lat, + ) + + # calculate angular separation of bin centers to FOV center + radius_bin_center = angular_separation( + bin_centers_grid_lon, + bin_centers_grid_lat, + 0, + 0, + ) + + # simple area mask with all bin centers outside FOV = 0 + mask_simple = np.logical_or( + radius_bin_center > self.viewcone_max, + radius_bin_center < self.viewcone_min, + ) + + area = np.ones(mask_simple.shape) + area[mask_simple] = 0 + + # select only bins partially covered by the FOV + bin_width_lon = bin_edges_lon[1] - bin_edges_lon[0] + bin_width_lat = bin_edges_lat[1] - bin_edges_lat[0] + bin_max_radius = np.sqrt(bin_width_lon ** 2 + bin_width_lat ** 2) * 0.5 + + fov_outer_edge_mask = np.logical_and( + radius_bin_center < self.viewcone_max + bin_max_radius, + radius_bin_center > self.viewcone_max - bin_max_radius, + ) + fov_inner_edge_mask = np.logical_and( + radius_bin_center < self.viewcone_min + bin_max_radius, + radius_bin_center > self.viewcone_min - bin_max_radius, + ) + + fov_edge_mask = np.logical_or(fov_inner_edge_mask, fov_outer_edge_mask) + + # get indices of relevant bin corners + corner_idx = np.nonzero(fov_edge_mask) + + # define start and endpoints for subpixels + edges_lon = np.array( + [bin_grid_lon[corner_idx], bin_grid_lon[corner_idx] + bin_width_lon] + ) + edges_lat = np.array( + [bin_grid_lat[corner_idx], bin_grid_lat[corner_idx] + bin_width_lat] + ) + + n_edge_pixels = len(fov_edge_mask[fov_edge_mask==True]) + + # define subpixel bin centers and grid + subpixel_lon = bin_center(np.linspace(*edges_lon, subpixels + 1) * u.deg) + subpixel_lat = bin_center(np.linspace(*edges_lat, subpixels + 1) * u.deg) + + # shape: (n_edge_pixels, 2, subpixels, subpixels) + subpixel_grid = np.array( + [np.meshgrid(subpixel_lon[:,i], subpixel_lat[:,i]) for i in range(n_edge_pixels)] + ) + subpixel_grid_lon = subpixel_grid[:,0] * u.deg + subpixel_grid_lat = subpixel_grid[:,1] * u.deg + + # make mask with subpixels inside the FOV + radius_subpixel = angular_separation(subpixel_grid_lon, subpixel_grid_lat, 0, 0) + mask = np.logical_and( + radius_subpixel <= self.viewcone_max, + radius_subpixel >= self.viewcone_min, + ) + + # calculates the fraction of subpixel centers within the FOV + FOV_covered_area = mask.sum(axis=(1,2)) / (subpixels ** 2) + + area[corner_idx] = FOV_covered_area + + return area \ No newline at end of file diff --git a/pyirf/tests/test_coordinates.py b/pyirf/tests/test_coordinates.py new file mode 100644 index 000000000..e96a639dd --- /dev/null +++ b/pyirf/tests/test_coordinates.py @@ -0,0 +1,60 @@ +from numpy.testing import assert_allclose +import astropy.units as u + +def test_gadf_fov_coords_lon_lat(): + from pyirf.coordinates import gadf_fov_coords_lon_lat + # test some simple cases + lon, lat = gadf_fov_coords_lon_lat(1 * u.deg, 1 * u.deg, 0 * u.deg, 0 * u.deg) + assert_allclose(lon.value, -1) + assert_allclose(lat.value, 1) + + lon, lat = gadf_fov_coords_lon_lat(269 * u.deg, 0 * u.deg, 270 * u.deg, 0 * u.deg) + assert_allclose(lon.value, 1) + assert_allclose(lat.value, 0, atol=1e-7) + + lon, lat = gadf_fov_coords_lon_lat(1 * u.deg, 60 * u.deg, 0 * u.deg, 60 * u.deg) + assert_allclose(lon.value, -0.5, rtol=1e-3) + assert_allclose(lat.value, 0.003779, rtol=1e-3) + + # these are cross-checked with the + # transformation as implemented in H.E.S.S. + az = [51.320575, 50.899125, 52.154053, 48.233023] + alt = [49.505451, 50.030165, 51.811739, 54.700102] + az_pointing = [52.42056255, 52.24706061, 52.06655505, 51.86795724] + alt_pointing = [51.11908203, 51.23454751, 51.35376141, 51.48385814] + lon, lat = gadf_fov_coords_lon_lat( + az * u.deg, alt * u.deg, az_pointing * u.deg, alt_pointing * u.deg + ) + assert_allclose( + lon.value, [0.7145614, 0.86603433, -0.05409698, 2.10295248], rtol=1e-5 + ) + assert_allclose( + lat.value, [-1.60829115, -1.19643974, 0.45800984, 3.26844192], rtol=1e-5 + ) + +def test_gadf_fov_coords_theta_phi(): + from pyirf.coordinates import gadf_fov_coords_theta_phi + + theta, phi = gadf_fov_coords_theta_phi( + lat=1 * u.deg, lon=0 * u.deg, pointing_lat=0 * u.deg, pointing_lon=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 0 * u.deg) + + theta, phi = gadf_fov_coords_theta_phi( + lat=-1 * u.deg, lon=0 * u.deg, pointing_lat=0 * u.deg, pointing_lon=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 180 * u.deg) + + theta, phi = gadf_fov_coords_theta_phi( + lat=0 * u.deg, lon=-1 * u.deg, pointing_lat=0 * u.deg, pointing_lon=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 90 * u.deg) + + theta, phi = gadf_fov_coords_theta_phi( + lat=0 * u.deg, lon=1 * u.deg, pointing_lat=0 * u.deg, pointing_lon=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 270 * u.deg) \ No newline at end of file diff --git a/pyirf/tests/test_simulations.py b/pyirf/tests/test_simulations.py index fbb5ac335..5ccabf827 100644 --- a/pyirf/tests/test_simulations.py +++ b/pyirf/tests/test_simulations.py @@ -106,6 +106,94 @@ def test_integrate_energy_fov_pointlike(): info.calculate_n_showers_per_energy_and_fov(energy_bins, fov_bins) +def test_integrate_3d_polar(): + from pyirf.simulations import SimulatedEventsInfo + + # simple case, max viewcone on bin edge + info = SimulatedEventsInfo( + n_showers=int(1e6), + energy_min=100 * u.GeV, + energy_max=10 * u.TeV, + max_impact=500 * u.m, + spectral_index=-2, + viewcone_min=0 * u.deg, + viewcone_max=10 * u.deg, + ) + + fov_bins_theta = [0, 10, 20] * u.deg + fov_bins_phi = [0, 180, 360] * u.deg + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) + + n_events = info.calculate_n_showers_3d_polar(energy_bins, fov_bins_theta, fov_bins_phi) + + assert np.all(n_events[:, 1:] == 0) + assert np.isclose(np.sum(n_events[...,0]), int(0.5e6)) + assert np.isclose(np.sum(n_events[...,1]), int(0.5e6)) + assert np.isclose(np.sum(n_events), int(1e6)) + + # viewcone inside of offset bin + info = SimulatedEventsInfo( + n_showers=int(1e6), + energy_min=100 * u.GeV, + energy_max=10 * u.TeV, + max_impact=500 * u.m, + spectral_index=-2, + viewcone_min=0 * u.deg, + viewcone_max=10 * u.deg, + ) + + fov_bins_theta = [0, 9, 11, 20] * u.deg + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) + n_events = info.calculate_n_showers_3d_polar(energy_bins, fov_bins_theta, fov_bins_phi) + + assert np.all(n_events[:, 1:2] > 0) + np.testing.assert_equal(n_events[:, 2:], 0) + assert np.isclose(np.sum(n_events), int(1e6)) + + +def test_integrate_3d_lonlat(): + from pyirf.simulations import SimulatedEventsInfo + + # simple case, max viewcone on bin edge + info = SimulatedEventsInfo( + n_showers=int(1e6), + energy_min=100 * u.GeV, + energy_max=10 * u.TeV, + max_impact=500 * u.m, + spectral_index=-2, + viewcone_min=0 * u.deg, + viewcone_max=10 * u.deg, + ) + + fov_bins = [-20, -10, 0, 10, 20] * u.deg + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) + + n_events = info.calculate_n_showers_3d_lonlat(energy_bins, fov_bins, fov_bins) + + assert np.all(n_events[:, 0, :] == 0) + assert np.all(n_events[:, :, 0] == 0) + assert np.allclose(np.sum(n_events, axis = 0)[1:3,1:3], int(0.25e6), rtol=1e-2) + assert np.isclose(np.sum(n_events), int(1e6), rtol=1e-2) + + # viewcone inside of offset bin + info = SimulatedEventsInfo( + n_showers=int(1e6), + energy_min=100 * u.GeV, + energy_max=10 * u.TeV, + max_impact=500 * u.m, + spectral_index=-2, + viewcone_min=0 * u.deg, + viewcone_max=10 * u.deg, + ) + + fov_bins = [-15, -5, 5, 15] * u.deg + energy_bins = np.geomspace(info.energy_min, info.energy_max, 20) + n_events = info.calculate_n_showers_3d_lonlat(energy_bins, fov_bins, fov_bins) + + assert np.all(n_events > 0) + assert np.isclose(np.sum(n_events), int(1e6), rtol=1e-2) + + def test_viewcone_integral(): from pyirf.simulations import _viewcone_pdf_integral diff --git a/pyirf/tests/test_utils.py b/pyirf/tests/test_utils.py index 83480ab2d..38186e521 100644 --- a/pyirf/tests/test_utils.py +++ b/pyirf/tests/test_utils.py @@ -30,6 +30,28 @@ def test_cone_solid_angle(): assert u.isclose(cone_solid_angle(0 * u.deg), 0 * u.sr) +def test_rectangle_solid_angle(): + from pyirf.utils import rectangle_solid_angle + + # whole sphere + assert u.isclose( + rectangle_solid_angle(0 * u.deg, 360 * u.deg, -90 * u.deg, 90 * u.deg), + 4 * np.pi * u.sr, + ) + + # half the sphere + assert u.isclose( + rectangle_solid_angle(0 * u.deg, 180 * u.deg, -90 * u.deg, 90 * u.deg), + 2 * np.pi * u.sr, + ) + + # zero + assert u.isclose( + rectangle_solid_angle(0 * u.deg, 0 * u.deg, 0* u.deg, 0 * u.deg), + 0 * u.sr, + ) + + def test_check_table(): from pyirf.exceptions import MissingColumns, WrongColumnUnit from pyirf.utils import check_table diff --git a/pyirf/utils.py b/pyirf/utils.py index d5ed730e5..4129fa0c3 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -1,6 +1,7 @@ import numpy as np import astropy.units as u from astropy.coordinates import angular_separation +from .coordinates import gadf_fov_coords_theta_phi, gadf_fov_coords_lon_lat 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,62 @@ 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 = gadf_fov_coords_theta_phi( + events[f"{prefix}_az"], + events[f"{prefix}_alt"], + events["pointing_az"], + events["pointing_alt"], + ) + + return phi.to(u.deg) + + +def calculate_source_fov_lonlat(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. + """ + lon, lat = gadf_fov_coords_lon_lat( + events[f"{prefix}_az"], + events[f"{prefix}_alt"], + events["pointing_az"], + events["pointing_alt"], + ) + + return lon.to(u.deg), lat.to(u.deg) + + def check_histograms(hist1, hist2, key="reco_energy"): """ Check if two histogram tables have the same binning @@ -128,6 +186,32 @@ def cone_solid_angle(angle): return solid_angle +def rectangle_solid_angle(lon_low, lon_high, lat_low, lat_high): + """Calculate the solid angle of a latitude-longitude rectangle + + Parameters + ---------- + lon_low: astropy.units.Quantity[angle] + Lower longitude coordinate of the rectangle corner + lat_low: astropy.units.Quantity[angle] + Lower latitude coordinate of the rectangle corner + lon_high: astropy.units.Quantity[angle] + Higher longitude coordinate of the rectangle corner + lat_high: astropy.units.Quantity[angle] + Higher Latitude coordinates of the rectangle corner + + Returns + ------- + solid angle: astropy.units.Quantity[solid angle] + + """ + diff_lon = (lon_high - lon_low).to_value(u.rad) + diff_lat = np.sin(lat_high.to_value(u.rad)) - np.sin(lat_low.to_value(u.rad)) + + solid_angle = diff_lon * diff_lat * u.sr + return solid_angle + + def check_table(table, required_columns=None, required_units=None): """Check a table for required columns and units.